Exponenta Banner
exponenta event banner

Расчет арктангенса с фиксированной точкой

Открыть сценарий

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

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

Введение

cordicatan2 функция аппроксимирует MATLAB ®atan2 с использованием алгоритма на основе CORDIC. CORDIC - аббревиатура от COORDinate Rotation DIgital Computer. Алгоритм CORDIC на основе ротации Givens (см. [1,2]) является одним из наиболее аппаратных эффективных алгоритмов, поскольку требует только итеративных операций добавления сдвига. Алгоритм CORDIC исключает необходимость явных умножителей и подходит для вычисления множества функций, таких как синус, косинус, арксинус, арккосинус, арктангенс, векторная величина, деление, квадратный корень, гиперболические и логарифмические функции.

Режим вычисления векторизации 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 (режим векторизации) выглядит следующим образом (для случая скаляра x, y, и 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 с помощью команды fiaccel MATLAB ®. Как правило, выполнение сгенерированной функции 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 и MATLAB atan2 функция, входные аргументы состоят из обоих x и y координаты вместо отношения y/x.

Для устранения ошибки квантования в приведенном ниже сравнении используются реализации с плавающей запятой алгоритмов аппроксимации полиномов КОРДИКА и Чебышева. Алгоритмическое сравнение ошибок показывает, что увеличение числа итераций 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-го порядка. Потому что динамический диапазон входов x, y и 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, Round и 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 объект в Designer™ Fixed-Point аппроксимирует плавающую точку построения MATLAB ®atan2 функция, использующая подход на основе одной таблицы подстановки с простой линейной интерполяцией ближайшего соседа между значениями. Этот подход позволяет создать небольшую таблицу поиска с действительным значением и использует простую арифметику.

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

Общие сведения о таблице подстановки ATAN2 Внедрение

Размер и точность таблицы подстановки

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

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

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

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

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

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

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

Остальные 8 младших значащих битов (LSB) 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, pi/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. Jack E. Volder, The CORDIC Trigonometric Computing Technique, IRE Transactions on Electronic Computers, том EC-8, сентябрь 1959, стр. 330-334.

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

%#ok<*NOPTS>

Документация конструктора фиксированных точек

Поддержка