Exponenta Banner

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

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

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

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

Введение

The cordicatan2 функция аппроксимирует MATLAB ® atan2 function, с использованием основанного на CORDIC алгоритма. CORDIC - это аббревиатура для COordinate Rotation DIgital Computer. Алгоритм CORDIC на основе вращения Givens (см. [1,2]) является одним из наиболее аппаратно эффективных алгоритмов, потому что он требует только итерационных операций shift-add. Алгоритм 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 Код

Введение

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

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

The 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 с вычислениями builtin 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.

The 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} $$

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

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

В целом более высокие степени полиномиального приближения дают более точные конечные результаты. Однако более высокие степени полиномиального приближения также увеличивают сложность алгоритма и требуют больших операций MAC и большей памяти. Чтобы соответствовать алгоритму CORDIC и MATLAB atan2 function, входные параметры состоят из обоих 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

Ошибка log base 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

Ошибка log base 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) Использование интерполяционных таблиц

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

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

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

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

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

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

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

Как компромисс, atan2 метод Fixed-Point Designer 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-битное (length-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] складываются или вычитаются с помощью констант, чтобы сформировать выходы угла с поправкой на октант.

Вычисление Arctangent с фиксированной точкой с использованием ATAN2

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

zFxpLUT = atan2(inYfix,inXfix);

Вычислим общую ошибку

Можно вычислить общую ошибку, сравнив базовое приближение интерполяционной таблицы с фиксированной точкой с builtin 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)

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

Как было сделано ранее, можно вычислить общую ошибку, сравнив приближения (приближения ) (и) с фиксированной точкой с builtin 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, Volume EC-8, September 1959, pp. 330-334.

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

%#ok<*NOPTS>

Документация Fixed-Point Designer

Поддержка