В этом примере показано, как использовать алгоритм CORDIC, полиномиальное приближение, и интерполяционная таблица приближается, чтобы вычислить фиксированную точку, четыре квадрантных обратных касательные. Эти реализации являются приближениями к встроенной функции MATLAB® atan2
. Эффективный алгоритм арктангенса фиксированной точки, чтобы оценить угол очень важен для многих приложений, включая управление робототехники, отслеживание частоты в радиосвязях, и многое другое.
atan2(y,x)
Используя алгоритм CORDICВведение
cordicatan2
функция аппроксимирует MATLAB® atan2
функция, с помощью основанного на CORDIC алгоритма. CORDIC является акронимом для Координатного Компьютера Вращения. Основанный на вращении алгоритм CORDIC Givens (см. [1,2]) является одним из большей части оборудования эффективные алгоритмы, потому что это только требует итеративных операций shift-add. Алгоритм CORDIC избавляет от необходимости явные множители и подходит для вычисления множества функций, таков как синус, косинус, arcsine, arccosine, арктангенс, векторная величина, разделитесь, квадратный корень, гиперболические и логарифмические функции.
CORDIC векторизация режима расчета
CORDIC векторизация уравнений режима широко используются, чтобы вычислить atan(y/x)
. В векторизации режима вращающее устройство CORDIC вращает входной вектор к положительной Оси X, чтобы минимизировать компонент вектора невязок. Для каждой итерации, если координата вектора невязок положительна, вращающее устройство CORDIC вращается по часовой стрелке (использующий отрицательный угол); в противном случае это вращается против часовой стрелки (использующий положительный угол). Если угловой аккумулятор инициализируется к 0, в конце итераций, накопленный угол поворота является углом исходного входного вектора.
В векторизации режима уравнения CORDIC:
угловой аккумулятор
где, если, и в противном случае;
, и общее количество итераций.
Как подходы:
Как объяснено выше, арктангенс может быть непосредственно вычислен с помощью вращающего устройства режима CORDIC векторизации с угловым аккумулятором, инициализированным, чтобы обнулить, i.e., и.
CORDICATAN2
КодВведение
cordicatan2
функция вычисляет четыре квадрантных арктангенса элементов X и Y, где. cordicatan2
вычисляет арктангенс с помощью алгоритма режима CORDIC векторизации, согласно вышеупомянутым уравнениям CORDIC.
Инициализация
cordicatan2
функция выполняет выполняющие шаги инициализации:
установлен в начальную букву X входных значений.
установлен в начальную букву Y входное значение.
обнуляется.
После итераций эти начальные значения приводят к
Разделяемая фиксированная точка и код ядра 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, как гарантируют, будет сходиться, но не всегда монотонно. Как будет показан в следующем примере, промежуточные итерации иногда вращают вектор ближе к положительной Оси 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 итераций (i.e., с точностью до 11 двоичных цифр), таким образом, величина ошибки меньше
cordic_algErr_bits = log2(cordic_algErr_real_world_value)
cordic_algErr_bits = -11.038839889583048
Отношение между количеством итераций и точностью
Если ошибка квантования доминирует над полной ошибкой, i.e., ошибка квантования больше алгоритмической ошибки, увеличивание общего числа итераций не значительно уменьшит полную ошибку фиксированной точки алгоритм 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-й порядок приближение Полинома Чебышева
Можно получить четыре квадрантных выхода посредством угловой коррекции на основе свойств функции арктангенса.
В общем случае более высокие степени полиномиального приближения приводят к более точным конечным результатам. Однако более высокие степени полиномиального приближения также увеличивают сложность алгоритма и требуют большего количества операций 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
вокруг
и 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 пробела (i.e., между углами 0 и радианами пи/4), затем можно выполнить коррекцию октанта на получившемся углу для любого y и x значений.
Как часть фрагмента предварительной обработки, рассматриваются знаки и относительные величины y и x, и деление выполняется. На основе знаков и величин y и x, вычисляется только одно из следующих значений: y/x, x/y,-y/x,-x/y,-y/-x,-x/-y. Результат без знака, который, как гарантируют, будет неотрицательным и чисто дробный, вычисляется, на основе априорного знания знаков и величин y и x. 16-битная дробная фиксированная точка без знака используется для этого значения.
8 старших значащих битов (MSBs) сохраненного представления беззнаковых целых чисел просто дробного результата фиксированной точки без знака затем используются, чтобы непосредственно индексировать 8-битное (длина 256) значение интерполяционной таблицы, содержащее угловые значения между 0 и радианы пи/4. Два поиска по таблице выполняется, один в вычисленном табличном местоположении индекса lutValBelow
, и один в следующем местоположении индекса lutValAbove
:
idxUint8MSBs = bitsliceget(idxUFIX16, 16, 9); zeroBasedIdx = int16(idxUint8MSBs); lutValBelow = FI_ATAN_LUT(zeroBasedIdx + 1); lutValAbove = FI_ATAN_LUT(zeroBasedIdx + 2);
Остающиеся 8 младших значащих битов (LSBs) idxUFIX16 используются, чтобы интерполировать между этими двумя табличными значениями. Значения LSB обработаны как нормированный масштабный коэффициент с 8-битным дробным типом данных rFracNT
:
rFracNT = numerictype(0,8,8); % fractional remainder data type
idxFrac8LSBs = reinterpretcast(bitsliceget(idxUFIX16,8,1), rFracNT);
rFraction = idxFrac8LSBs;
Два значения интерполяционной таблицы, с остатком (rFraction) значение, используются, чтобы выполнить простую линейную интерполяцию ближайшего соседа. Действительное умножается, используется, чтобы определить взвешенное различие между двумя точками. Это приводит к простому вычислению (эквивалентный одному продукту и двум суммам), чтобы получить интерполированный результат фиксированной точки:
temp = rFraction * (lutValAbove - lutValBelow); rslt = lutValBelow + temp;
Наконец, на основе исходных знаков и относительных величин y и x, выходной результат формируется с помощью простой логики коррекции октанта и арифметики. Первый октант [0, пи/4] угловые результаты значения добавлены или вычтены с константами, чтобы сформировать откорректированный октантом угол выходные параметры.
ATAN2
Можно вызвать atan2
функционируйте непосредственно использующий фиксированную точку или входные параметры с плавающей точкой. Основанный на интерполяционной таблице алгоритм используется для фиксированной точки atan2
реализация:
zFxpLUT = atan2(inYfix,inXfix);
Вычислите полную ошибку
Можно вычислить полную ошибку путем сравнения фиксированной точки основанное на интерполяционной таблице приближение со встроенным atan2
функция. Идеальным ссылочным выходом является zfltRef
.
lut_err = bsxfun(@minus, zfltRef, double(zFxpLUT));
max_lutErr_real_world_value = max(abs(lut_err'));
max_lutErr_bits = log2(max_lutErr_real_world_value);
fprintf('Overall error in bits: %g\n', max_lutErr_bits);
Overall error in bits: -12.6743
figno = figno + 1; fidemo.fixpt_atan2_demo_plot(figno, theta, lut_err)
Как был сделан ранее, можно вычислить полную ошибку путем сравнения приближения (приближений) фиксированной точки со встроенным atan2
функция. Идеальным ссылочным выходом является zfltRef
.
zfixCDC15 = cordicatan2(inYfix, inXfix, 15); cordic_15I_err = bsxfun(@minus, zfltRef, double(zfixCDC15)); poly_7p_err = bsxfun(@minus, zfltRef, double(zfix7p)); figno = figno + 1; fidemo.fixpt_atan2_demo_plot(figno, theta, cordic_15I_err, poly_7p_err, lut_err)
Фиксированная точка алгоритм CORDIC требует следующих операций:
1 поиск по таблице на итерацию
2 сдвига на итерацию
3 сложения на итерацию
Алгоритм аппроксимации Полинома Чебышева фиксированной точки N-го порядка требует следующих операций:
1 деление
(N+1) умножение
(N-1)/2 сложения
Упрощенный один алгоритм интерполяционной таблицы с линейной интерполяцией ближайшего соседа требует следующих операций:
1 деление
2 поиска по таблице
1 умножение
2 сложения
В приложениях реального мира, выбирая алгоритм для вычисления арктангенса фиксированной точки обычно зависит от требуемой точности, стоимости и аппаратных ограничений.
close all; % close all figure windows
Джек Э. Волдер, Тригонометрический Вычислительный Метод CORDIC, Транзакции IRE на Электронно-вычислительных машинах, Volume EC 8, сентябрь 1959, стр 330-334.
Рэй Андрэка, обзор алгоритма CORDIC для основанных на FPGA компьютеров, Продолжений 1998 шестых международных симпозиумов ACM/SIGDA по Программируемым пользователем вентильным матрицам, 22-24 февраля 1998, стр 191-200.
%#ok<*NOPTS>