Документация

Вычисление 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, чтобы минимизировать компонент вектора невязок. Для каждой итерации, если координата вектора невязок положительна, вращающее устройство CORDIC вращается по часовой стрелке (использующий отрицательный угол); в противном случае это вращается против часовой стрелки (использующий положительный угол). Если угловой аккумулятор инициализируется к 0, в конце итераций, накопленный угол поворота является углом исходного входного вектора.

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

угловой аккумулятор

где, если, и в противном случае;

, и общее количество итераций.

Как подходы:

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

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

Введение

Функция cordicatan2 вычисляет четыре квадрантных арктангенса элементов X и Y, где. cordicatan2 вычисляет арктангенс с помощью алгоритма режима CORDIC векторизации, согласно вышеупомянутому уравнения CORDIC.

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

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

После итераций эти начальные значения приводят к

Разделяемая фиксированная точка и код ядра 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 не будет постоянным, потому что отличался бы.

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

Пример

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

Инициализируйте входной вектор с угловыми степенями, значение = 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 подходы, усиление вращающего устройства CORDIC приближается 1.64676. В этом примере вход был на модульном круге, таким образом, начальное значение вращающего устройства равняется 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

Обратите внимание на то, что подходы 0 и подходы, потому что.

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

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

Вычислите алгоритмическую ошибку 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 двоичных цифр), таким образом, значение ошибки является меньше, чем

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) является полиномом приближения. В целом можно получить полиномы очень близко к оптимальному путем приближения заданной функции с точки зрения Полиномов Чебышева и отключения полинома в желаемой степени.

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

где

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

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

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

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

Сравнение алгоритмической ошибки 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 имеет предполагаемую угловую точность:

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 битов (лучшая точность). Алгоритм устойчив, потому что в [-1, +1], и умножение коэффициентов и в (-1, +1). Таким образом динамический диапазон не вырастет, и из-за предопределенных типов данных с фиксированной точкой, переполнение не ожидается.

Подобно алгоритму CORDIC четыре квадрантных основанных на полиномиальном приближении алгоритма atan2 выводят, оценил углы в. Поэтому мы можем выбрать выходную длину части 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 в Fixed-Point Designer™ аппроксимирует MATLAB® встроенная функция atan2 с плавающей точкой, с помощью одного основанного на интерполяционной таблице подхода с простой линейной интерполяцией ближайшего соседа между значениями. Этот подход допускает маленькую интерполяционную таблицу с действительным знаком и использует простую арифметику.

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

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

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

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

Как компромисс, метод 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 требует следующих операций:

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

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

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

close all; % close all figure windows

Ссылки

%#ok<*NOPTS>