Вычислите арктангенс фиксированной точки

Открыть скрипт

В этом примере показано, как использовать алгоритм CORDIC, полиномиальное приближение, и интерполяционная таблица приближается, чтобы вычислить фиксированную точку, четыре квадрантных обратных касательные. Эти реализации являются приближениями к встроенной функции MATLAB® atan2. Эффективный алгоритм арктангенса фиксированной точки, чтобы оценить угол очень важен для многих приложений, включая управление робототехники, отслеживание частоты в радиосвязях, и многое другое.

Вычисление atan2(y,x) Используя алгоритм CORDIC

Введение

cordicatan2 функция аппроксимирует MATLAB® atan2 функция, с помощью основанного на CORDIC алгоритма. CORDIC является акронимом для Координатного Компьютера Вращения. Основанный на вращении алгоритм CORDIC Givens (см. [1,2]) является одним из большей части оборудования эффективные алгоритмы, потому что это только требует итеративных операций shift-add. Алгоритм CORDIC избавляет от необходимости явные множители и подходит для вычисления множества функций, таков как синус, косинус, arcsine, arccosine, арктангенс, векторная величина, разделитесь, квадратный корень, гиперболические и логарифмические функции.

CORDIC векторизация режима расчета

CORDIC векторизация уравнений режима широко используются, чтобы вычислить atan(y/x). В векторизации режима вращающее устройство CORDIC вращает входной вектор к положительной Оси X, чтобы минимизировать$$ y $$ компонент вектора невязок. Для каждой итерации, если$$ y $$ координата вектора невязок положительна, вращающее устройство CORDIC вращается по часовой стрелке (использующий отрицательный угол); в противном случае это вращается против часовой стрелки (использующий положительный угол). Если угловой аккумулятор инициализируется к 0, в конце итераций, накопленный угол поворота является углом исходного входного вектора.

В векторизации режима уравнения CORDIC:

$$ x_{i+1} = x_{i} - y_{i}*d_{i}*2^{-i} $$

$$ y_{i+1} = y_{i} + x_{i}*d_{i}*2^{-i} $$

$$ z_{i+1} = z_{i} + d_{i}*atan(2^{-i}) $$ угловой аккумулятор

где$$ d_{i} = +1 $$, если$$ y_{i} < 0 $$, и$$ -1 $$ в противном случае;

$$ i = 0, 1, ..., N-1 $$, и$$ N $$ общее количество итераций.

Как$$ N $$ подходы$$ +\infty $$:

$$ x_{N} = A_{N}\sqrt{x_{0}^2+y_{0}^2} $$

$$ y_{N} = 0 $$

$$ z_{N} = z_{0} + atan(y_{0}/x_{0}) $$

$$ A_{N} =&#xA;1/(cos(atan(2^{0}))*cos(atan(2^{-1}))*...*cos(atan(2^{-(N-1)})))&#xA; = \prod_{i=0}^{N-1}{\sqrt{1+2^{-2i}}}&#xA; $$

Как объяснено выше, арктангенс может быть непосредственно вычислен с помощью вращающего устройства режима CORDIC векторизации с угловым аккумулятором, инициализированным, чтобы обнулить, т.е.$$ z_{0}=0, $$ и$$ z_{N} \approx atan(y_{0}/x_{0}) $$.

Понимание CORDICATAN2 Код

Введение

cordicatan2 функция вычисляет четыре квадрантных арктангенса элементов X и Y, где$$ -\pi \leq ATAN2(y,x) \leq +\pi $$. cordicatan2 вычисляет арктангенс с помощью алгоритма режима CORDIC векторизации, согласно вышеупомянутому уравнения CORDIC.

Инициализация

cordicatan2 функция выполняет выполняющие шаги инициализации:

  • $$ x_{0} $$ установлен в начальную букву X входных значений.

  • $$ y_{0} $$ установлен в начальную букву Y входное значение.

  • $$ z_{0} $$ обнуляется.

После$$ N $$ итераций эти начальные значения приводят к $$ z_{N} \approx atan(y_{0}/x_{0}) $$

Разделяемая фиксированная точка и код ядра CORDIC с плавающей точкой

Код MATLAB для алгоритма CORDIC (векторизовавший режим) фрагмент ядра можно следующим образом (для случая скалярного xY, и z). Этот тот же код используется и в фиксированной точке и в операциях с плавающей точкой:

function [x, y, z] = cordic_vectoring_kernel(x, y, z, inpLUT, n)
% Perform CORDIC vectoring kernel algorithm for N kernel iterations.
xtmp = x;
ytmp = y;
for idx = 1:n
    if y < 0
        x(:) = x - ytmp;
        y(:) = y + xtmp;
        z(:) = z - inpLUT(idx);
    else
        x(:) = x + ytmp;
        y(:) = y - xtmp;
        z(:) = z + inpLUT(idx);
    end
    xtmp = bitsra(x, idx); % bit-shift-right for multiply by 2^(-idx)
    ytmp = bitsra(y, idx); % bit-shift-right for multiply by 2^(-idx)
end

Визуализация режима векторизации итерации CORDIC

Алгоритм CORDIC обычно запускается через заданное (постоянное) количество итераций начиная с окончания итераций CORDIC, рано повредил бы конвейерный код, и усиление CORDIC$$ A_{n} $$ не будет постоянным, потому что$$ n $$ варьировался бы.

Для очень больших значений $$ n $$алгоритм CORDIC, как гарантируют, будет сходиться, но не всегда монотонно. Как будет показан в следующем примере, промежуточные итерации иногда вращают вектор ближе к положительной Оси X, чем следующая итерация. Можно обычно достигать большей точности путем увеличения общего числа итераций.

Пример

В следующем примере итерация 5 обеспечивает лучшую оценку угла, чем итерация 6, и алгоритм CORDIC сходится в более поздних итерациях.

Инициализируйте входной вектор угловыми$$ \theta = 43 $$ степенями, величина = 1

origFormat = get(0, 'format'); % store original format setting;
                               % restore this setting at the end.
format short
%
theta = 43*pi/180;  % input angle in radians
Niter = 10;         % number of iterations
inX   = cos(theta); % x coordinate of the input vector
inY   = sin(theta); % y coordinate of the input vector
%
% pre-allocate memories
zf = zeros(1, Niter);
xf = [inX, zeros(1, Niter)];
yf = [inY, zeros(1, Niter)];
angleLUT = atan(2.^-(0:Niter-1)); % pre-calculate the angle lookup table
%
% Call CORDIC vectoring kernel algorithm
for k = 1:Niter
   [xf(k+1), yf(k+1), zf(k)] = fixed.internal.cordic_vectoring_kernel_private(inX, inY, 0, angleLUT, k);
end

Следующий выход показывает угловое накопление CORDIC (в градусах) через 10 итераций. Обратите внимание на то, что 5-я итерация произвела меньше ошибки, чем 6-я итерация, и что расчетный угол быстро сходится к фактическому входному углу позже.

angleAccumulator = zf*180/pi; angleError = angleAccumulator - theta*180/pi;
fprintf('Iteration: %2d, Calculated angle: %7.3f, Error in degrees: %10g, Error in bits: %g\n',...
        [(1:Niter); angleAccumulator(:)'; angleError(:)';log2(abs(zf(:)'-theta))]);

Iteration:  1, Calculated angle:  45.000, Error in degrees:          2, Error in bits: -4.84036
Iteration:  2, Calculated angle:  18.435, Error in degrees:   -24.5651, Error in bits: -1.22182
Iteration:  3, Calculated angle:  32.471, Error in degrees:   -10.5288, Error in bits: -2.44409
Iteration:  4, Calculated angle:  39.596, Error in degrees:   -3.40379, Error in bits: -4.07321
Iteration:  5, Calculated angle:  43.173, Error in degrees:   0.172543, Error in bits: -8.37533
Iteration:  6, Calculated angle:  41.383, Error in degrees:   -1.61737, Error in bits: -5.14671
Iteration:  7, Calculated angle:  42.278, Error in degrees:  -0.722194, Error in bits: -6.3099
Iteration:  8, Calculated angle:  42.725, Error in degrees:   -0.27458, Error in bits: -7.70506
Iteration:  9, Calculated angle:  42.949, Error in degrees: -0.0507692, Error in bits: -10.1403
Iteration: 10, Calculated angle:  43.061, Error in degrees:  0.0611365, Error in bits: -9.87218

Как N подходы$$ +\infty $$, усиление вращающего устройства CORDIC$$ A_{N} $$ приближается 1.64676. В этом примере вход$$ (x_{0},y_{0}) $$ был на модульном круге, таким образом, начальная величина вращающего устройства равняется 1. Следующий выход показывает величину вращающего устройства через 10 итераций:

rotatorMagnitude = sqrt(xf.^2+yf.^2); % CORDIC rotator gain through iterations
fprintf('Iteration: %2d, Rotator magnitude: %g\n',...
    [(0:Niter); rotatorMagnitude(:)']);

Iteration:  0, Rotator magnitude: 1
Iteration:  1, Rotator magnitude: 1.41421
Iteration:  2, Rotator magnitude: 1.58114
Iteration:  3, Rotator magnitude: 1.6298
Iteration:  4, Rotator magnitude: 1.64248
Iteration:  5, Rotator magnitude: 1.64569
Iteration:  6, Rotator magnitude: 1.64649
Iteration:  7, Rotator magnitude: 1.64669
Iteration:  8, Rotator magnitude: 1.64674
Iteration:  9, Rotator magnitude: 1.64676
Iteration: 10, Rotator magnitude: 1.64676

Обратите внимание на то, что$y_{n}$ подходы 0 и$x_{n}$ подходы$$ A_{n} \sqrt{x_{0}^{2} + y_{0}^{2}} = A_{n}, $$, потому что$$ \sqrt{x_{0}^{2} + y_{0}^{2}} = 1 $$.

y_n = yf(end)

y_n =

   -0.0018


x_n = xf(end)

x_n =

    1.6468


figno = 1;
fidemo.fixpt_atan2_demo_plot(figno, xf, yf) %Vectoring Mode CORDIC Iterations


figno = figno + 1; %Cumulative Angle and Rotator Magnitude Through Iterations
fidemo.fixpt_atan2_demo_plot(figno,Niter, theta, angleAccumulator, rotatorMagnitude)

Выполнение полного анализа ошибок алгоритма CORDIC

Полная ошибка состоит из двух частей:

  1. Алгоритмическая ошибка, которая следует из угла поворота CORDIC, представляемого конечным числом основных углов.

  2. Квантование или погрешность округления, которая следует из конечного представления точности угловой интерполяционной таблицы, и от конечной арифметики точности, используемой в операциях фиксированной точки.

Вычислите алгоритмическую ошибку CORDIC

theta  = (-178:2:180)*pi/180; % angle in radians
inXflt = cos(theta); % generates input vector
inYflt = sin(theta);
Niter  = 12; % total number of iterations
zflt   = cordicatan2(inYflt, inXflt, Niter); % floating-point results

Вычислите максимальную величину алгоритмической ошибки CORDIC путем сравнения расчета CORDIC со встроенным atan2 функция.

format long
cordic_algErr_real_world_value = max(abs((atan2(inYflt, inXflt) - zflt)))

cordic_algErr_real_world_value =

     4.753112306290497e-04

Логарифмическая основа 2 ошибки связана с количеством итераций. В этом примере мы используем 12 итераций (т.е. с точностью до 11 двоичных цифр), таким образом, величина ошибки меньше $$ 2^{-11} $$

cordic_algErr_bits = log2(cordic_algErr_real_world_value)

cordic_algErr_bits =

 -11.038839889583048

Отношение между количеством итераций и точностью

Если ошибка квантования доминирует над полной ошибкой, т.е. ошибка квантования больше алгоритмической ошибки, увеличивание общего числа итераций не значительно уменьшит полную ошибку фиксированной точки алгоритм CORDIC. Необходимо выбрать дробные длины и общее количество итераций, чтобы гарантировать, что ошибка квантования меньше, чем алгоритмическая ошибка. В алгоритме CORDIC точность увеличивает на один бит каждую итерацию. Таким образом нет никакой причины выбрать много итераций, больше, чем точность входных данных.

Другой способ посмотреть на отношение между количеством итераций и точностью находится на шаге сдвига вправо алгоритма. Например, на против часовой стрелки вращении

x(:) = x0 - bitsra(y,i);
y(:) = y + bitsra(x0,i);

если i равна размеру слова y и x0, то bitsra(y,i) и bitsra(x0,i) переключите полностью, чтобы обнулить и ничего не внести в следующий шаг.

Чтобы измерить ошибку из алгоритма фиксированной точки, а не различия во входных значениях, вычисляют ссылку с плавающей точкой с теми же входными параметрами как фиксированная точка алгоритм CORDIC.

inXfix = sfi(inXflt, 16, 14);
inYfix = sfi(inYflt, 16, 14);
zref   = atan2(double(inYfix), double(inXfix));
zfix8  = cordicatan2(inYfix, inXfix, 8);
zfix10 = cordicatan2(inYfix, inXfix, 10);
zfix12 = cordicatan2(inYfix, inXfix, 12);
zfix14 = cordicatan2(inYfix, inXfix, 14);
zfix15 = cordicatan2(inYfix, inXfix, 15);
cordic_err = bsxfun(@minus,zref,double([zfix8;zfix10;zfix12;zfix14;zfix15]));

Ошибка зависит от количества итераций и точности входных данных. В вышеупомянутом примере входные данные находятся в области значений [-1, +1], и дробная длина равняется 14. Из следующих таблиц, показывающих максимальную погрешность в каждой итерации и фигуру, показывающую полную ошибку алгоритма CORDIC, вы видите, что ошибочные уменьшения приблизительно на 1 бит за итерацию до точности данных достигнуты.

iterations = [8, 10, 12, 14, 15];
max_cordicErr_real_world_value = max(abs(cordic_err'));
fprintf('Iterations: %2d, Max error in real-world-value: %g\n',...
    [iterations; max_cordicErr_real_world_value]);

Iterations:  8, Max error in real-world-value: 0.00773633
Iterations: 10, Max error in real-world-value: 0.00187695
Iterations: 12, Max error in real-world-value: 0.000501175
Iterations: 14, Max error in real-world-value: 0.000244621
Iterations: 15, Max error in real-world-value: 0.000244621

max_cordicErr_bits = log2(max_cordicErr_real_world_value);
fprintf('Iterations: %2d, Max error in bits: %g\n',[iterations; max_cordicErr_bits]);

Iterations:  8, Max error in bits: -7.01414
Iterations: 10, Max error in bits: -9.05739
Iterations: 12, Max error in bits: -10.9624
Iterations: 14, Max error in bits: -11.9972
Iterations: 15, Max error in bits: -11.9972

figno = figno + 1;
fidemo.fixpt_atan2_demo_plot(figno, theta, cordic_err)

Ускорение фиксированной точки CORDICATAN2 Алгоритм Используя FIACCEL

Можно сгенерировать MEX-функцию из кода MATLAB с помощью MATLAB® fiaccel команда. Как правило, выполнение сгенерированной MEX-функции может улучшить скорость симуляции, несмотря на то, что фактическое улучшение скорости зависит от используемой платформы симуляции. Следующий пример показывает, как ускорить фиксированную точку cordicatan2 алгоритм с помощью fiaccel.

fiaccel функционируйте компилирует код MATLAB в MEX-функцию. Этот шаг требует создания временной директории и полномочий записи в той директории.

tempdirObj = fidemo.fiTempdir('fixpt_atan2_demo');

Когда вы объявляете, что количество итераций константа (например, 12) использование coder.newtype('constant',12), скомпилированная угловая интерполяционная таблица также будет постоянной, и таким образом не будет вычислена в каждой итерации. Кроме того, когда вы вызываете скомпилированный файл MEX cordicatan2_mex, вы не должны будете давать ему входной параметр для количества итераций. Если вы передадите в количестве итераций, MEX-функция будет ошибка.

Тип данных входных параметров определяет ли cordicatan2 функция выполняет фиксированную точку или вычисления с плавающей точкой. Когда MATLAB генерирует код для этого файла, код только сгенерирован для определенного типа данных. Например, если входные параметры являются фиксированной точкой, только фиксированная точка сгенерирована.

inp = {inYfix, inXfix, coder.newtype('constant',12)}; % example inputs for the function
fiaccel('cordicatan2', '-o', 'cordicatan2_mex',  '-args', inp)

Во-первых, вычислите вектор 4 квадрантных atan2 путем вызова cordicatan2.

tstart = tic;
cordicatan2(inYfix,inXfix,Niter);
telapsed_Mcordicatan2 = toc(tstart);

Затем вычислите вектор 4 квадрантных atan2 путем вызова MEX-функции cordicatan2_mex.

cordicatan2_mex(inYfix,inXfix); % load the MEX file
tstart = tic;
cordicatan2_mex(inYfix,inXfix);
telapsed_MEXcordicatan2 = toc(tstart);

Теперь сравните скорость. Введите следующее в окне команды MATLAB, чтобы видеть улучшение скорости на вашей определенной платформе:

fiaccel_speedup = telapsed_Mcordicatan2/telapsed_MEXcordicatan2;

Чтобы очистить временную директорию, запустите следующие команды:

clear cordicatan2_mex;
status = tempdirObj.cleanUp;

Вычисление atan2(y,x) Используя приближение полинома Чебышева

Полиномиальное приближение является умножением - накапливают (MAC) центральный алгоритм. Это может быть хороший выбор для реализаций DSP нелинейных функций как atan(x).

Для данной степени полинома и заданной функции f(x) = atan(x) оцененный на интервале [-1, +1], теория полиномиального приближения пытается найти полином, который минимизирует максимальное значение$$ |P(x)-f(x)| $$, где P(x) полином аппроксимации. В общем случае можно получить полиномы очень близко к оптимальному путем аппроксимации заданной функции в терминах Полиномов Чебышева и отключения полинома в желаемой степени.

Приближение арктангенса на интервале [-1, +1] использование Полинома Чебышева первого вида получено в итоге в следующей формуле:

$$ atan(x) = 2\sum_{n=0}^{\infty} {(-1)^{n}q^{2n+1} \over (2n+1)}&#xA;T_{2n+1}(x) $$

где

$$ q = 1/(1+\sqrt{2}) $$

$$ x \in [-1, +1] $$

$$ T_{0}(x) = 1 $$

$$ T_{1}(x) = x $$

$$ T_{n+1}(x) = 2xT_{n}(x) - T_{n-1}(x). $$

Поэтому 3-й порядок приближение Полинома Чебышева

$$ atan(x) = 0.970562748477141*x - 0.189514164974601*x^{3}. $$

5-й порядок приближение Полинома Чебышева

$$ atan(x) = 0.994949366116654*x - 0.287060635532652*x^{3}&#xA; + 0.078037176446441*x^{5}. $$

7-й порядок приближение Полинома Чебышева

$$ \begin{array}{lllll}&#xA; atan(x) &#38; = &#38; 0.999133448222780*x &#38; - &#38; 0.320533292381664*x^{3} \\&#xA; &#38; + &#38; 0.144982490144465*x^{5} &#38; - &#38; 0.038254464970299*x^{7}.&#xA;\end{array} $$

Можно получить четыре квадрантных выхода посредством угловой коррекции на основе свойств функции арктангенса.

Сравнение алгоритмической ошибки CORDIC и алгоритмов полиномиального приближения

В общем случае более высокие степени полиномиального приближения приводят к более точным конечным результатам. Однако более высокие степени полиномиального приближения также увеличивают сложность алгоритма и требуют большего количества операций MAC и большей памяти. Быть сопоставимым с алгоритмом CORDIC и atan2 MATLAB функция, входные параметры состоят из обоих x и y координаты вместо отношения y/x.

Чтобы устранить ошибку квантования, реализации с плавающей точкой CORDIC и алгоритмов аппроксимации Полинома Чебышева используются в сравнении ниже. Алгоритмическое ошибочное сравнение показывает, что увеличение числа итераций CORDIC приводит к меньшей ошибке. Это также показывает, что алгоритм CORDIC с 12 итерациями обеспечивает немного лучшую угловую оценку, чем 5-й порядок приближение Полинома Чебышева. Ошибка приближения 3-го порядка Полином Чебышева приблизительно в 8 раз больше, чем тот из 5-го порядка Полином Чебышева. Необходимо выбрать порядок или степень полинома на основе требуемой точности угловой оценки и аппаратных ограничений.

Коэффициенты приближения Полинома Чебышева для atan(x) показаны в порядке возрастания x.

constA3 = [0.970562748477141, -0.189514164974601]; % 3rd order
constA5 = [0.994949366116654,-0.287060635532652,0.078037176446441]; % 5th order
constA7 = [0.999133448222780 -0.320533292381664 0.144982490144465...
          -0.038254464970299]; % 7th order

theta   = (-90:1:90)*pi/180; % angle in radians
inXflt  = cos(theta);
inYflt  = sin(theta);
zfltRef = atan2(inYflt, inXflt); % Ideal output from ATAN2 function
zfltp3  = fidemo.poly_atan2(inYflt,inXflt,3,constA3); % 3rd order polynomial
zfltp5  = fidemo.poly_atan2(inYflt,inXflt,5,constA5); % 5th order polynomial
zfltp7  = fidemo.poly_atan2(inYflt,inXflt,7,constA7); % 7th order polynomial
zflt8   = cordicatan2(inYflt, inXflt,  8); % CORDIC alg with 8 iterations
zflt12  = cordicatan2(inYflt, inXflt, 12); % CORDIC alg with 12 iterations

Максимальную алгоритмическую ошибочную величину (или норма по бесконечности алгоритмической ошибки) для алгоритма CORDIC с 8 и 12 итерациями показывают ниже:

cordic_algErr    = [zfltRef;zfltRef] - [zflt8;zflt12];
max_cordicAlgErr = max(abs(cordic_algErr'));
fprintf('Iterations: %2d, CORDIC algorithmic error in real-world-value: %g\n',...
    [[8,12]; max_cordicAlgErr(:)']);

Iterations:  8, CORDIC algorithmic error in real-world-value: 0.00772146
Iterations: 12, CORDIC algorithmic error in real-world-value: 0.000483258

Логарифмическая основа 2 ошибки показывает количество двоичных цифр точности. 12-я итерация алгоритма CORDIC имеет предполагаемую угловую точность$$ 2^{-11} $$:

max_cordicAlgErr_bits = log2(max_cordicAlgErr);
fprintf('Iterations: %2d, CORDIC algorithmic error in bits: %g\n',...
    [[8,12]; max_cordicAlgErr_bits(:)']);

Iterations:  8, CORDIC algorithmic error in bits: -7.01691
Iterations: 12, CORDIC algorithmic error in bits: -11.0149

Следующий код показывает величину максимальной алгоритмической ошибки полиномиального приближения для порядков 3, 5, и 7:

poly_algErr    = [zfltRef;zfltRef;zfltRef] - [zfltp3;zfltp5;zfltp7];
max_polyAlgErr = max(abs(poly_algErr'));
fprintf('Order: %d, Polynomial approximation algorithmic error in real-world-value: %g\n',...
    [3:2:7; max_polyAlgErr(:)']);

Order: 3, Polynomial approximation algorithmic error in real-world-value: 0.00541647
Order: 5, Polynomial approximation algorithmic error in real-world-value: 0.000679384
Order: 7, Polynomial approximation algorithmic error in real-world-value: 9.16204e-05

Логарифмическая основа 2 ошибки показывает количество двоичных цифр точности.

max_polyAlgErr_bits = log2(max_polyAlgErr);
fprintf('Order: %d, Polynomial approximation algorithmic error in bits: %g\n',...
    [3:2:7; max_polyAlgErr_bits(:)']);

Order: 3, Polynomial approximation algorithmic error in bits: -7.52843
Order: 5, Polynomial approximation algorithmic error in bits: -10.5235
Order: 7, Polynomial approximation algorithmic error in bits: -13.414

figno = figno + 1;
fidemo.fixpt_atan2_demo_plot(figno, theta, cordic_algErr, poly_algErr)

Преобразование алгоритма аппроксимации полинома Чебышева с плавающей точкой для фиксированной точки

Примите, что размеры слова ввода и вывода ограничиваются к 16 битам оборудованием, и 5-й порядок Полином Чебышева используется в приближении. Поскольку динамический диапазон входных параметров xY и y/x все в [-1, +1], можно избежать переполнения путем выбора типа входных данных фиксированной точки со знаком с размером слова 16 битов и дробной длиной 14 битов. Коэффициенты полинома чисто дробны и в (-1, +1), таким образом, мы можем выбрать их типы данных как подписанную фиксированную точку с размером слова 16 битов и дробной длиной 15 битов (лучшая точность). Алгоритм устойчив, потому что$$ (y/x)^{n} $$ в [-1, +1], и умножение коэффициентов и$$ (y/x)^{n} $$ в (-1, +1). Таким образом динамический диапазон не вырастет, и из-за предопределенных типов данных с фиксированной точкой, переполнение не ожидается.

Подобно алгоритму CORDIC, четырем квадрантным основанным на полиномиальном приближении atan2 алгоритм выходные параметры оценил углы в$$ [-\pi, \pi] $$. Поэтому мы можем выбрать выходную длину части 13 битов, чтобы избежать переполнения и предоставить динамический спектр [-4, +3.9998779296875].

Основное приближение Полинома Чебышева с плавающей точкой арктангенса на интервале [-1, +1] реализовано как chebyPoly_atan_fltpt локальная функция в poly_atan2.m файл.

   function z = chebyPoly_atan_fltpt(y,x,N,constA,Tz,RoundingMethodStr)
   tmp = y/x;
   switch N
       case 3
           z = constA(1)*tmp + constA(2)*tmp^3;
       case 5
           z = constA(1)*tmp + constA(2)*tmp^3 + constA(3)*tmp^5;
       case 7
           z = constA(1)*tmp + constA(2)*tmp^3 + constA(3)*tmp^5 + constA(4)*tmp^7;
       otherwise
           disp('Supported order of Chebyshev polynomials are 3, 5 and 7');
   end

Основное приближение Полинома Чебышева фиксированной точки арктангенса на интервале [-1, +1] реализовано как chebyPoly_atan_fixpt локальная функция в poly_atan2.m файл.

   function z = chebyPoly_atan_fixpt(y,x,N,constA,Tz,RoundingMethodStr)
   z = fi(0,'numerictype', Tz, 'RoundingMethod', RoundingMethodStr);
   Tx = numerictype(x);
   tmp = fi(0, 'numerictype',Tx, 'RoundingMethod', RoundingMethodStr);
   tmp(:) = Tx.divide(y, x); % y/x;
   tmp2 = fi(0, 'numerictype',Tx, 'RoundingMethod', RoundingMethodStr);
   tmp3 = fi(0, 'numerictype',Tx, 'RoundingMethod', RoundingMethodStr);
   tmp2(:) = tmp*tmp;  % (y/x)^2
   tmp3(:) = tmp2*tmp; % (y/x)^3
   z(:) = constA(1)*tmp + constA(2)*tmp3; % for order N = 3
   if (N == 5) || (N == 7)
       tmp5 = fi(0, 'numerictype',Tx, 'RoundingMethod', RoundingMethodStr);
       tmp5(:) = tmp3 * tmp2; % (y/x)^5
       z(:) = z + constA(3)*tmp5; % for order N = 5
       if N == 7
           tmp7 = fi(0, 'numerictype',Tx, 'RoundingMethod', RoundingMethodStr);
           tmp7(:) = tmp5 * tmp2; % (y/x)^7
           z(:) = z + constA(4)*tmp7; %for order N = 7
       end
   end

Универсальные четыре квадрантных atan2 вычисление с помощью приближения Полинома Чебышева реализовано в poly_atan2.m файл.

   function z = poly_atan2(y,x,N,constA,Tz,RoundingMethodStr)
   if nargin < 5
       % floating-point algorithm
       fhandle = @chebyPoly_atan_fltpt;
       Tz = [];
       RoundingMethodStr = [];
       z = zeros(size(y));
   else
       % fixed-point algorithm
       fhandle = @chebyPoly_atan_fixpt;
       %pre-allocate output
       z = fi(zeros(size(y)), 'numerictype', Tz, 'RoundingMethod', RoundingMethodStr);
   end
   % Apply angle correction to obtain four quadrant output
   for idx = 1:length(y)
      % first quadrant
      if abs(x(idx)) >= abs(y(idx))
          % (0, pi/4]
          z(idx) = feval(fhandle, abs(y(idx)), abs(x(idx)), N, constA, Tz, RoundingMethodStr);
      else
          % (pi/4, pi/2)
          z(idx) = pi/2 - feval(fhandle, abs(x(idx)), abs(y(idx)), N, constA, Tz, RoundingMethodStr);
      end
      if x(idx) < 0
          % second and third quadrant
          if y(idx) < 0
              z(idx) = -pi + z(idx);
          else
             z(idx) = pi - z(idx);
          end
      else % fourth quadrant
          if y(idx) < 0
              z(idx) = -z(idx);
          end
      end
   end

Выполнение полного анализа ошибок алгоритма полиномиального приближения

Подобно алгоритму CORDIC полная ошибка алгоритма полиномиального приближения состоит из двух частей - алгоритмическая ошибка и ошибка квантования. Алгоритмическая ошибка алгоритма полиномиального приближения анализировалась и по сравнению с алгоритмической ошибкой алгоритма CORDIC в предыдущем разделе.

Вычислите ошибку квантования

Вычислите ошибку квантования путем сравнения полиномиального приближения фиксированной точки с полиномиальным приближением с плавающей точкой.

Квантуйте входные параметры и коэффициенты с конвергентным округлением:

inXfix = fi(fi(inXflt,  1, 16, 14,'RoundingMethod','Convergent'),'fimath',[]);
inYfix = fi(fi(inYflt,  1, 16, 14,'RoundingMethod','Convergent'),'fimath',[]);
constAfix3 = fi(fi(constA3, 1, 16,'RoundingMethod','Convergent'),'fimath',[]);
constAfix5 = fi(fi(constA5, 1, 16,'RoundingMethod','Convergent'),'fimath',[]);
constAfix7 = fi(fi(constA7, 1, 16,'RoundingMethod','Convergent'),'fimath',[]);

Вычислите максимальную величину ошибки квантования Floor округление:

ord    = 3:2:7; % using 3rd, 5th, 7th order polynomials
Tz     = numerictype(1, 16, 13); % output data type
zfix3p = fidemo.poly_atan2(inYfix,inXfix,ord(1),constAfix3,Tz,'Floor'); % 3rd order
zfix5p = fidemo.poly_atan2(inYfix,inXfix,ord(2),constAfix5,Tz,'Floor'); % 5th order
zfix7p = fidemo.poly_atan2(inYfix,inXfix,ord(3),constAfix7,Tz,'Floor'); % 7th order
poly_quantErr = bsxfun(@minus, [zfltp3;zfltp5;zfltp7], double([zfix3p;zfix5p;zfix7p]));
max_polyQuantErr_real_world_value = max(abs(poly_quantErr'));
max_polyQuantErr_bits = log2(max_polyQuantErr_real_world_value);
fprintf('PolyOrder: %2d, Quant error in bits: %g\n',...
    [ord; max_polyQuantErr_bits]);

PolyOrder:  3, Quant error in bits: -12.7101
PolyOrder:  5, Quant error in bits: -12.325
PolyOrder:  7, Quant error in bits: -11.8416

Вычислите полную ошибку

Вычислите полную ошибку путем сравнения полиномиального приближения фиксированной точки со встроенным atan2 функция. Идеальным ссылочным выходом является zfltRef. Полная ошибка 7-го полиномиального приближения порядка во власти ошибки квантования, которая происходит из-за конечной точности входных данных, коэффициентов и округляющихся эффектов от операций вычислений с фиксированной точкой.

poly_err = bsxfun(@minus, zfltRef, double([zfix3p;zfix5p;zfix7p]));
max_polyErr_real_world_value = max(abs(poly_err'));
max_polyErr_bits = log2(max_polyErr_real_world_value);
fprintf('PolyOrder: %2d, Overall error in bits: %g\n',...
    [ord; max_polyErr_bits]);

PolyOrder:  3, Overall error in bits: -7.51907
PolyOrder:  5, Overall error in bits: -10.2497
PolyOrder:  7, Overall error in bits: -11.5883

figno = figno + 1;
fidemo.fixpt_atan2_demo_plot(figno, theta, poly_err)

Эффект округления режимов в полиномиальном приближении

По сравнению с алгоритмом CORDIC с 12 итерациями и 13-битной дробной длиной в угловом аккумуляторе, пятый порядок приближение Полинома Чебышева дает подобное распоряжение ошибки квантования. В следующем примере, Nearestвокруг и Convergent округляющиеся режимы дают меньшие ошибки квантования, чем Floor округление режима.

Максимальная величина ошибки квантования Floor округление

poly5_quantErrFloor = max(abs(poly_quantErr(2,:)));
poly5_quantErrFloor_bits = log2(poly5_quantErrFloor)

poly5_quantErrFloor_bits =

 -12.324996933210334

Для сравнения вычислите максимальную величину ошибки квантования Nearest округление:

zfixp5n = fidemo.poly_atan2(inYfix,inXfix,5,constAfix5,Tz,'Nearest');
poly5_quantErrNearest = max(abs(zfltp5 - double(zfixp5n)));
poly5_quantErrNearest_bits = log2(poly5_quantErrNearest)
set(0, 'format', origFormat); % reset MATLAB output format

poly5_quantErrNearest_bits =

 -13.175966487895451

Вычисление atan2(y,x) Используя интерполяционные таблицы

Существует много основанных на интерполяционной таблице подходов, которые могут использоваться, чтобы реализовать приближения арктангенса фиксированной точки. Следующее является недорогим подходом на основе одной интерполяционной таблицы с действительным знаком и простой линейной интерполяцией ближайшего соседа.

Один основанный на интерполяционной таблице подход

atan2 метод fi объект в Fixed-Point Designer™ аппроксимирует MATLAB® встроенный atan2 с плавающей точкой функция, с помощью одного основанного на интерполяционной таблице подхода с простой линейной интерполяцией ближайшего соседа между значениями. Этот подход допускает маленькую интерполяционную таблицу с действительным знаком и использует простую арифметику.

Используя одну интерполяционную таблицу с действительным знаком упрощает расчет индекса и полную арифметику, требуемую достигнуть очень хорошей точности результатов. Эти упрощения дают к относительно быстродействию, а также относительно низким требованиям к памяти.

Понимание основанного на интерполяционной таблице ATAN2 Реализация

Размер интерполяционной таблицы и точность

Два важных конструктивных соображения интерполяционной таблицы являются ее размером и ее точностью. Не возможно составить таблицу для каждого возможного$$ y/x $$ входного значения. Также не возможно быть совершенно точным из-за квантования значений интерполяционной таблицы.

Как компромисс, atan2 метод Fixed-Point Designer fi возразите использует 8-битную интерполяционную таблицу в качестве части ее реализации. 8-битная таблица является только 256 элементами долго, таким образом, она мала и эффективна. Восемь битов также соответствуют размеру байта или слова на многих платформах. Используемый в сочетании с линейной интерполяцией и 16-битным выходом (значение интерполяционной таблицы) точность, интерполяционная таблица с побитовой адресацией 8 обеспечивает очень хорошую точность, а также производительность.

Обзор реализации алгоритма

Чтобы лучше изучить реализацию Fixed-Point Designer, сначала рассмотрите симметрию atan2(y,x) с четырьмя квадрантами функция. Если вы всегда вычисляете арктангенс в первом октанте x-y пробела (т.е. между углами 0 и радианами пи/4), то можно выполнить коррекцию октанта на получившемся углу для любого y и x значений.

Как часть фрагмента предварительной обработки, рассматриваются знаки и относительные величины y и x, и деление выполняется. На основе знаков и величин y и x, вычисляется только одно из следующих значений: y/x, x/y,-y/x,-x/y,-y/-x,-x/-y. Результат без знака, который, как гарантируют, будет неотрицательным и чисто дробный, вычисляется, на основе априорного знания знаков и величин y и x. 16-битная дробная фиксированная точка без знака используется в этом значении.

8 старших значащих битов (MSBs) сохраненного представления беззнаковых целых чисел просто дробного результата фиксированной точки без знака затем используются, чтобы непосредственно индексировать 8-битное (длина 256) значение интерполяционной таблицы, содержащее угловые значения между 0 и радианы пи/4. Два поиска по таблице выполняется, один в вычисленном табличном местоположении индекса lutValBelow, и один в следующем местоположении индекса lutValAbove:

idxUint8MSBs = bitsliceget(idxUFIX16, 16, 9);
zeroBasedIdx = int16(idxUint8MSBs);
lutValBelow  = FI_ATAN_LUT(zeroBasedIdx + 1);
lutValAbove  = FI_ATAN_LUT(zeroBasedIdx + 2);

Остающиеся 8 младших значащих битов (LSBs) idxUFIX16 используются, чтобы интерполировать между этими двумя табличными значениями. Значения LSB обработаны как нормированный масштабный коэффициент с 8-битным дробным типом данных rFracNT:

rFracNT      = numerictype(0,8,8); % fractional remainder data type
idxFrac8LSBs = reinterpretcast(bitsliceget(idxUFIX16,8,1), rFracNT);
rFraction    = idxFrac8LSBs;

Два значения интерполяционной таблицы, с остатком (rFraction) значение, используются, чтобы выполнить простую линейную интерполяцию ближайшего соседа. Действительное умножается, используется, чтобы определить взвешенное различие между двумя точками. Это приводит к простому вычислению (эквивалентный одному продукту и двум суммам), чтобы получить интерполированный результат фиксированной точки:

temp = rFraction * (lutValAbove - lutValBelow);
rslt = lutValBelow + temp;

Наконец, на основе исходных знаков и относительных величин y и x, выходной результат формируется с помощью простой логики коррекции октанта и арифметики. Первый октант [0, пи/4] угловые результаты значения добавлены или вычтены с константами, чтобы сформировать откорректированный октантом угол выходные параметры.

Вычисление арктангенса фиксированной точки Используя ATAN2

Можно вызвать atan2 функционируйте непосредственно использующий фиксированную точку или входные параметры с плавающей точкой. Основанный на интерполяционной таблице алгоритм используется в фиксированной точке atan2 реализация:

zFxpLUT = atan2(inYfix,inXfix);

Вычислите полную ошибку

Можно вычислить полную ошибку путем сравнения фиксированной точки основанное на интерполяционной таблице приближение со встроенным atan2 функция. Идеальным ссылочным выходом является zfltRef.

lut_err = bsxfun(@minus, zfltRef, double(zFxpLUT));
max_lutErr_real_world_value = max(abs(lut_err'));
max_lutErr_bits = log2(max_lutErr_real_world_value);
fprintf('Overall error in bits: %g\n', max_lutErr_bits);

Overall error in bits: -12.6743

figno = figno + 1;
fidemo.fixpt_atan2_demo_plot(figno, theta, lut_err)

Сравнение полной ошибки между реализациями фиксированной точки

Как был сделан ранее, можно вычислить полную ошибку путем сравнения приближения (приближений) фиксированной точки со встроенным atan2 функция. Идеальным ссылочным выходом является zfltRef.

zfixCDC15      = cordicatan2(inYfix, inXfix, 15);
cordic_15I_err = bsxfun(@minus, zfltRef, double(zfixCDC15));
poly_7p_err    = bsxfun(@minus, zfltRef, double(zfix7p));
figno = figno + 1;
fidemo.fixpt_atan2_demo_plot(figno, theta, cordic_15I_err, poly_7p_err, lut_err)

Сравнение затрат на алгоритмы аппроксимации фиксированной точки

Фиксированная точка алгоритм CORDIC требует следующих операций:

  • 1 поиск по таблице на итерацию

  • 2 сдвига на итерацию

  • 3 сложения на итерацию

Алгоритм аппроксимации Полинома Чебышева фиксированной точки N-го порядка требует следующих операций:

  • 1 деление

  • (N+1) умножение

  • (N-1)/2 сложения

Упрощенный один алгоритм интерполяционной таблицы с линейной интерполяцией ближайшего соседа требует следующих операций:

  • 1 деление

  • 2 поиска по таблице

  • 1 умножение

  • 2 сложения

В приложениях реального мира, выбирая алгоритм для вычисления арктангенса фиксированной точки обычно зависит от требуемой точности, стоимости и аппаратных ограничений.

close all; % close all figure windows

Ссылки

  1. Джек Э. Волдер, Тригонометрический Вычислительный Метод CORDIC, Транзакции IRE на Электронно-вычислительных машинах, Volume EC 8, сентябрь 1959, стр 330-334.

  2. Рэй Андрэка, обзор алгоритма CORDIC для основанных на FPGA компьютеров, Продолжений 1998 шестых международных симпозиумов ACM/SIGDA по Программируемым пользователем вентильным матрицам, 22-24 февраля 1998, стр 191-200.

%#ok<*NOPTS>
Документация Fixed-Point Designer
Поддержка
Для просмотра документации необходимо авторизоваться на сайте