Этот пример показывает, как использовать алгоритм 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, чтобы минимизировать компонент вектора невязок. Для каждой итерации, если координата вектора невязок положительная, вращатель CORDIC вращается по часовой стрелке (с использованием отрицательного угла); в противном случае он вращается против часовой стрелки (с помощью положительного угла). Если аккумулятор угла инициализирован равным 0, в конце итераций накопленный угол поворота является углом исходного входного вектора.
В режиме векторизации уравнения CORDIC:
- аккумулятор угла
где если и в противном случае;
, и общее количество итераций.
По мере приближения:
Как объяснено выше, арктангенс может быть непосредственно вычислен с помощью вращателя CORDIC в режиме векторизации с угловым аккумулятором, инициализированным в нуль, т.е. и.
CORDICATAN2
КодВведение
The cordicatan2
функция вычисляет четыре квадрантных арктангенса элементов x и y, где. cordicatan2
вычисляет арктангенс с помощью алгоритма CORDIC в режиме векторизации, согласно вышеприведенным уравнениям CORDIC.
Инициализация
The 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 с вычислениями 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 двоичных цифр), поэтому величина ошибки меньше, чем
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)
- аппроксимирующий полином. В целом можно получить полиномы очень близкую к оптимальной путем аппроксимации заданной функции с точки зрения Полиномов Чебышева и отсечения полинома в нужной степени.
Приближение арктангенса через интервал [-1, + 1] с помощью многочлена Чебышёва первого рода суммируется в следующей формуле:
где
Поэтому аппроксимация Полинома Чебышева 3-го порядка
Аппроксимация Полинома Чебышева 5-го порядка
Полином Чебышева 7-го порядка
Можно получить выход четырех квадрантов через коррекцию угла на основе свойств функции arctangent.
В целом более высокие степени полиномиального приближения дают более точные конечные результаты. Однако более высокие степени полиномиального приближения также увеличивают сложность алгоритма и требуют больших операций 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 имеет оцененную точность угла:
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 битов (лучшая точность). Алгоритм является устойчивым, потому что находится в пределах [-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)
Использование интерполяционных таблицСуществует много подходов, основанных на интерполяционной таблице, которые могут использоваться для реализации арктических приближений с фиксированной точкой. Ниже приведен дешевый подход, основанный на одной реальной интерполяционной таблице и простой линейной интерполяции по ближайшему соседу.
The atan2
метод fi
объект в Fixed-Point Designer™ аппроксимирует atan2
с плавающей точкой MATLAB ® builtin функция, с использованием подхода, основанного на одной интерполяционной таблице, с простой линейной интерполяцией между значениями по ближайшему соседу. Этот подход позволяет создать небольшую реальную интерполяционную таблицу и использует простую арифметику.
Использование одной реальной интерполяционной таблицы упрощает расчет индекса и общую арифметику, необходимую для достижения очень хорошей точности результатов. Эти упрощения обеспечивают относительно высокую скорость эффективности, а также относительно низкие требования к памяти.
ATAN2
РеализацияРазмер и точность интерполяционной таблицы
Двумя важными факторами проекта интерполяционной таблицы являются ее размер и ее точность. Невозможно создать таблицу для всех возможных входных значений. Также невозможно быть совершенно точным из-за квантования значений интерполяционной таблицы.
Как компромисс, 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] складываются или вычитаются с помощью констант, чтобы сформировать выходы угла с поправкой на октант.
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
Jack E. Volder, The CORDIC Trigonometric Computing Technique, IRE Transactions on Electronic Computers, Volume EC-8, September 1959, pp. 330-334.
Рэй Андрака, Обзор алгоритма CORDIC для компьютеров на основе FPGA, Материалы шестого международного симпозиума ACM/SIGDA 1998 года по программируемым массивам, 22-24 февраля 1998 года, стр. 191-200.
%#ok<*NOPTS>