В этом примере показано, как использовать и основанные на CORDIC и основанные на интерполяционной таблице алгоритмы, предоставленные Fixed-Point Designer™, чтобы аппроксимировать синус MATLAB® (SIN
) и косинус (COS
Функции. Эффективный синус фиксированной точки и алгоритмы косинуса очень важны для многих встраиваемых приложений, включая блоки управления приводом, навигацию, обработку сигналов и радиосвязи.
Введение
cordiccexp
, cordicsincos
, cordicsin
, и cordiccos
функции аппроксимируют sin
MATLAB и
cos
функции с помощью основанного на CORDIC алгоритма. CORDIC является акронимом для Координатного Компьютера Вращения. Основанный на вращении алгоритм CORDIC Givens (см. [1,2]) является одним из большей части оборудования эффективные алгоритмы, потому что это только требует итеративных операций shift-add. Алгоритм CORDIC избавляет от необходимости явные множители и подходит для вычисления множества функций, таков как синус, косинус, arcsine, arccosine, арктангенс, векторная величина, разделитесь, квадратный корень, гиперболические и логарифмические функции.
Можно использовать вращение CORDIC вычислительный режим, чтобы вычислить синус и косинус, и также полярные-к-декартову операции преобразования. В этом режиме векторная величина и угол вращения известны и координата (X-Y), компоненты вычисляются после вращения.
Режим расчета вращения CORDIC
Алгоритм режима вращения CORDIC начинается путем инициализации углового аккумулятора желаемым углом поворота. Затем решение вращения в каждой итерации CORDIC сделано способом, который уменьшает величину остаточного углового аккумулятора. Решение вращения основано на знаке остаточного угла в угловом аккумуляторе после каждой итерации.
В режиме вращения уравнения CORDIC:
где, если, и в противном случае;
, и общее количество итераций.
Это обеспечивает следующий результат как подходы:
Где:
.
В режиме вращения алгоритм CORDIC ограничивается углами поворота между и. Поддерживать углы за пределами той области значений, cordiccexp
, cordicsincos
, cordicsin
, и cordiccos
функции используют квадрантную коррекцию (включая возможное дополнительное отрицание) после того, как итерации CORDIC будут завершены.
CORDICSINCOS
Синус и код косинусаВведение
cordicsincos
функция вычисляет синус и косинус входных углов в области значений [-2*pi, 2*pi), использование алгоритма CORDIC. Эта функция берет угол (радианы) и количество итераций как входные параметры. Функция возвращает приближения синуса и косинуса.
Расчет CORDIC выходные параметры масштабируется усилением вращающего устройства. Это усиление составляется путем предварительного масштабирования начального постоянного значения.
Инициализация
cordicsincos
функция выполняет выполняющие шаги инициализации:
Угол ввел интерполяционную таблицу inpLUT
установлен в atan(2 .^ -(0:N-1))
.
установлен в значение входного параметра.
установлен в.
обнуляется.
Разумный выбор начальных значений позволяет алгоритму непосредственно вычислять и синус и косинус одновременно. После итераций эти начальные значения приводят к следующим выходным параметрам как подходы:
Разделяемая фиксированная точка и код ядра CORDIC с плавающей точкой
Код MATLAB для алгоритма CORDIC (режим вращения) фрагмент ядра можно следующим образом (для случая скалярного x
Y
, и z
). Этот тот же код используется и в фиксированной точке и в операциях с плавающей точкой:
function [x, y, z] = cordic_rotation_kernel(x, y, z, inpLUT, n) % Perform CORDIC rotation kernel algorithm for N kernel iterations. xtmp = x; ytmp = y; for idx = 1:n if z < 0 z(:) = z + inpLUT(idx); x(:) = x + ytmp; y(:) = y - xtmp; else z(:) = z - inpLUT(idx); x(:) = x - ytmp; y(:) = y + xtmp; 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, как гарантируют, будет сходиться, но не всегда монотонно. Как будет показан в следующем примере, промежуточные итерации иногда приводят к более точным результатам, чем более поздние итерации. Можно обычно достигать большей точности путем увеличения общего числа итераций.
Пример
В следующем примере итерация 5 обеспечивает лучшую оценку результата, чем итерация 6, и алгоритм CORDIC сходится в более поздних итерациях.
theta = pi/5; % input angle in radians niters = 10; % number of iterations sinTh = sin(theta); % reference result cosTh = cos(theta); % reference result y_sin = zeros(niters, 1); sin_err = zeros(niters, 1); x_cos = zeros(niters, 1); cos_err = zeros(niters, 1); fprintf('\n\nNITERS \tERROR\n'); fprintf('------\t----------\n'); for n = 1:niters [y_sin(n), x_cos(n)] = cordicsincos(theta, n); sin_err(n) = abs(y_sin(n) - sinTh); cos_err(n) = abs(x_cos(n) - cosTh); if n < 10 fprintf(' %d \t %1.8f\n', n, cos_err(n)); else fprintf(' %d \t %1.8f\n', n, cos_err(n)); end end fprintf('\n');
NITERS ERROR ------ ---------- 1 0.10191021 2 0.13966630 3 0.03464449 4 0.03846157 5 0.00020393 6 0.01776952 7 0.00888037 8 0.00436052 9 0.00208192 10 0.00093798
Постройте ошибку приближения CORDIC на столбчатом графике
figure(1); clf; bar(1:niters, cos_err(1:niters)); xlabel('Number of iterations','fontsize',12,'fontweight','b'); ylabel('Error','fontsize',12,'fontweight','b'); title('CORDIC approximation error for cos(pi/5) computation',... 'fontsize',12,'fontweight','b'); axis([0 niters 0 0.14]);
Постройте результаты X-Y для 5 итераций
Niter2Draw = 5; figure(2), clf, hold on plot(cos(0:0.1:pi/2), sin(0:0.1:pi/2), 'b--'); % semi-circle for i=1:Niter2Draw plot([0 x_cos(i)],[0 y_sin(i)], 'LineWidth', 2); % CORDIC iteration result text(x_cos(i),y_sin(i),int2str(i),'fontsize',12,'fontweight','b'); end plot(cos(theta), sin(theta), 'r*', 'MarkerSize', 20); % IDEAL result xlabel('X (COS)','fontsize',12,'fontweight','b') ylabel('Y (SIN)','fontsize',12,'fontweight','b') title('CORDIC iterations for cos(pi/5) computation',... 'fontsize',12,'fontweight','b') axis equal; axis square;
cordicsin
Создайте 1 024 точки между [-2*pi, 2*pi),
stepSize = pi/256; thRadDbl = (-2*pi):stepSize:(2*pi - stepSize); thRadFxp = sfi(thRadDbl, 12); % signed, 12-bit fixed-point values sinThRef = sin(double(thRadFxp)); % reference results
Сравните фиксированную точку CORDIC по сравнению с аккуратными функциональными результатами с двойной точностью
Используйте 12-битные квантованные входные параметры и варьируйтесь количество итераций от 4 до 10.
for niters = 4:3:10 cdcSinTh = cordicsin(thRadFxp, niters); errCdcRef = sinThRef - double(cdcSinTh); figure; hold on; axis([-2*pi 2*pi -1.25 1.25]); plot(thRadFxp, sinThRef, 'b'); plot(thRadFxp, cdcSinTh, 'g'); plot(thRadFxp, errCdcRef, 'r'); ylabel('sin(\Theta)','fontsize',12,'fontweight','b'); set(gca,'XTick',-2*pi:pi/2:2*pi); set(gca,'XTickLabel',... {'-2*pi', '-3*pi/2', '-pi', '-pi/2', ... '0', 'pi/2', 'pi', '3*pi/2','2*pi'}); set(gca,'YTick',-1:0.5:1); set(gca,'YTickLabel',{'-1.0','-0.5','0','0.5','1.0'}); ref_str = 'Reference: sin(double(\Theta))'; cdc_str = sprintf('12-bit CORDICSIN; N = %d', niters); err_str = sprintf('Error (max = %f)', max(abs(errCdcRef))); legend(ref_str, cdc_str, err_str); title(cdc_str,'fontsize',12,'fontweight','b'); end
Вычислите ошибку LSB для N = 10
figure; fracLen = cdcSinTh.FractionLength; plot(thRadFxp, abs(errCdcRef) * pow2(fracLen)); set(gca,'XTick',-2*pi:pi/2:2*pi); set(gca,'XTickLabel',... {'-2*pi', '-3*pi/2', '-pi', '-pi/2', ... '0', 'pi/2', 'pi', '3*pi/2','2*pi'}); ylabel(sprintf('LSB Error: 1 LSB = 2^{-%d}',fracLen),'fontsize',12,'fontweight','b'); title('LSB Error: 12-bit CORDICSIN; N=10','fontsize',12,'fontweight','b'); axis([-2*pi 2*pi 0 6]);
Вычислите уровень шума
fft_mag = abs(fft(double(cdcSinTh))); max_mag = max(fft_mag); mag_db = 20*log10(fft_mag/max_mag); figure; hold on; plot(0:1023, mag_db); plot(0:1023, zeros(1,1024),'r--'); % Normalized peak (0 dB) plot(0:1023, -62.*ones(1,1024),'r--'); % Noise floor level ylabel('dB Magnitude','fontsize',12,'fontweight','b'); title('62 dB Noise Floor: 12-bit CORDICSIN; N=10',... 'fontsize',12,'fontweight','b'); % axis([0 1023 -120 0]); full FFT axis([0 round(1024*(pi/8)) -100 10]); % zoom in set(gca,'XTick',[0 round(1024*pi/16) round(1024*pi/8)]); set(gca,'XTickLabel',{'0','pi/16','pi/8'});
CORDICSINCOS
Функция с FIACCEL
Можно сгенерировать MEX-функцию из кода MATLAB с помощью MATLAB® fiaccel функция. Как правило, выполнение сгенерированной MEX-функции может улучшить скорость симуляции, несмотря на то, что фактическое улучшение скорости зависит от используемой платформы симуляции. Следующий пример показывает, как ускорить фиксированную точку cordicsincos
функция с помощью fiaccel
.
fiaccel
функционируйте компилирует код MATLAB в MEX-функцию. Этот шаг требует создания временной директории и полномочий записи в этой директории.
tempdirObj = fidemo.fiTempdir('fi_sin_cos_demo');
Когда вы объявляете, что количество итераций константа (например, 10
) использование coder.newtype('constant',10)
, скомпилированная угловая интерполяционная таблица также будет постоянной, и таким образом не будет вычислена в каждой итерации. Кроме того, когда вы вызываете cordicsincos_mex
, вы не должны будете давать ему входной параметр для количества итераций. Если вы передадите в количестве итераций, MEX-функция будет ошибка.
Тип данных входных параметров определяет ли cordicsincos
функция выполняет фиксированную точку или вычисления с плавающей точкой. Когда MATLAB генерирует код для этого файла, код только сгенерирован для определенного типа данных. Например, если входной параметр THETA является фиксированной точкой, то только фиксированная точка сгенерирована.
inp = {thRadFxp, coder.newtype('constant',10)}; % example inputs for the function fiaccel('cordicsincos', '-o', 'cordicsincos_mex', '-args', inp)
Во-первых, вычислите синус и косинус путем вызова cordicsincos
.
tstart = tic; cordicsincos(thRadFxp,10); telapsed_Mcordicsincos = toc(tstart);
Затем вычислите синус и косинус путем вызова MEX-функции cordicsincos_mex
.
cordicsincos_mex(thRadFxp); % load the MEX file
tstart = tic;
cordicsincos_mex(thRadFxp);
telapsed_MEXcordicsincos = toc(tstart);
Теперь сравните скорость. Введите следующее в командной строке MATLAB, чтобы видеть улучшение скорости на вашей платформе:
fiaccel_speedup = telapsed_Mcordicsincos/telapsed_MEXcordicsincos;
Чтобы очистить временную директорию, запустите следующие команды:
clear cordicsincos_mex;
status = tempdirObj.cleanUp;
SIN
и COS
Используя интерполяционные таблицыСуществует много основанных на интерполяционной таблице подходов, которые могут использоваться, чтобы реализовать синус фиксированной точки и приближения косинуса. Следующее является недорогим подходом на основе одной интерполяционной таблицы с действительным знаком и простой линейной интерполяцией ближайшего соседа.
sin
и cos
методы fi
объект в Fixed-Point Designer аппроксимирует MATLAB® встроенный sin
с плавающей точкой и
cos
функции, с помощью основанного на интерполяционной таблице подхода с простой линейной интерполяцией ближайшего соседа между значениями. Этот подход допускает маленькую интерполяционную таблицу с действительным знаком и использует простую арифметику.
Используя одну интерполяционную таблицу с действительным знаком упрощает расчет индекса и полную арифметику, требуемую достигнуть очень хорошей точности результатов. Эти упрощения дают к относительно быстродействию и также относительно низким требованиям к памяти.
SIN
и COS
РеализацияРазмер интерполяционной таблицы и точность
Два важных конструктивных соображения интерполяционной таблицы являются ее размером и ее точностью. Не возможно составить таблицу для каждого возможного входного значения. Также не возможно быть совершенно точным из-за квантования или значений интерполяционной таблицы.
Как компромисс, Fixed-Point Designer SIN
и COS
методы FI
используйте 8-битную интерполяционную таблицу в качестве части их реализации. 8-битная таблица является только 256 элементами долго, таким образом, она мала и эффективна. Восемь битов также соответствуют размеру байта или слова на многих платформах. Используемый в сочетании с линейной интерполяцией и 16-битным выходом (значение интерполяционной таблицы) точность, интерполяционная таблица с побитовой адресацией 8 обеспечивает и очень хорошую точность и производительность.
Инициализация постоянного SIN
Значения интерполяционной таблицы
Для простоты реализации, однородности табличного значения и скорости, используется полная sinewave таблица. Во-первых, SIN
волны четвертью функция производится в 64 универсальных интервалах в области значений [0, пи/2), радианы. Выбор 16-битного дробного типа данных с фиксированной точкой со знаком для табличных значений, т.е.
tblValsNT = numerictype(1,16,15)
, приводит к лучшим результатам точности в SIN
выведите область значений [-1.0, 1.0). Значения предварительно квантуются, прежде чем они будут установлены, чтобы избежать предупреждений переполнения.
tblValsNT = numerictype(1,16,15); quarterSinDblFltPtVals = (sin(2*pi*((0:63) ./ 256)))'; endpointQuantized_Plus1 = 1.0 - double(eps(fi(0,tblValsNT))); halfSinWaveDblFltPtVals = ... [quarterSinDblFltPtVals; ... endpointQuantized_Plus1; ... flipud(quarterSinDblFltPtVals(2:end))]; fullSinWaveDblFltPtVals = ... [halfSinWaveDblFltPtVals; -halfSinWaveDblFltPtVals]; FI_SIN_LUT = fi(fullSinWaveDblFltPtVals, tblValsNT);
Обзор реализации алгоритма
Реализация Fixed-Point Designer sin
и cos
методы fi
объекты включают сначала кастинг угловых входных параметров фиксированной точки (в радианах) к предопределенному типу данных в области значений [0, 2pi]. С этой целью операция по-модулю-2pi выполняется, чтобы получить входное значение фиксированной точки inpValInRange
в области значений [0, 2pi] и бросок к лучшей двоичной точке точности масштабировал 16-битную фиксированную точку без знака numerictype(0,16,13)
:
% Best UNSIGNED type for real-world value range [0, 2*pi], % which maps to fixed-point stored integer vals [0, 51472]. inpInRangeNT = numerictype(0,16,13);
Затем мы получаем 16-битное сохраненное значение беззнаковых целых чисел от этой фиксированной точки в области значений угловое значение FI:
idxUFIX16 = fi(storedInteger(inpValInRange), numerictype(0,16,0));
Мы умножаем сохраненное целочисленное значение на постоянную нормализацию, 65536/51472. Получившееся целочисленное значение будет в полномасштабной области значений индекса uint16:
normConst_NT = numerictype(0,32,31); normConstant = fi(65536/51472, normConst_NT); fullScaleIdx = normConstant * idxUFIX16; idxUFIX16(:) = fullScaleIdx;
Лучшие 8 старших значащих битов (MSBs) этого полномасштабного 16-битного индекса без знака idxUFIX16
используются, чтобы непосредственно индексировать в 8-битную интерполяционную таблицу синуса. Два поиска по таблице выполняется, один в вычисленном табличном местоположении индекса lutValBelow
, и один в следующем местоположении индекса lutValAbove
:
idxUint8MSBs = storedInteger(bitsliceget(idxUFIX16, 16, 9)); zeroBasedIdx = int16(idxUint8MSBs); lutValBelow = FI_SIN_LUT(zeroBasedIdx + 1); lutValAbove = FI_SIN_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;
Действительное умножается, используется, чтобы определить взвешенное различие между двумя точками. Это приводит к простому вычислению (эквивалентный одному продукту и двум суммам), чтобы получить интерполированный результат синуса фиксированной точки:
temp = rFraction * (lutValAbove - lutValBelow); rslt = lutValBelow + temp;
Пример
Используя вышеупомянутый алгоритм, вот пример интерполяционной таблицы, и процесс линейной интерполяции использовался для расчета значения SIN
для входа inpValInRange = 0.425
фиксированной точки радианы:
% Use an arbitrary input value (e.g., 0.425 radians) inpInRangeNT = numerictype(0,16,13); % best precision, [0, 2*pi] radians inpValInRange = fi(0.425, inpInRangeNT); % arbitrary fixed-point input angle % Normalize its stored integer to get full-scale unsigned 16-bit integer index idxUFIX16 = fi(storedInteger(inpValInRange), numerictype(0,16,0)); normConst_NT = numerictype(0,32,31); normConstant = fi(65536/51472, normConst_NT); fullScaleIdx = normConstant * idxUFIX16; idxUFIX16(:) = fullScaleIdx; % Do two table lookups using unsigned 8-bit integer index (i.e., 8 MSBs) idxUint8MSBs = storedInteger(bitsliceget(idxUFIX16, 16, 9)); zeroBasedIdx = int16(idxUint8MSBs); % zero-based table index value lutValBelow = FI_SIN_LUT(zeroBasedIdx + 1); % 1st table lookup value lutValAbove = FI_SIN_LUT(zeroBasedIdx + 2); % 2nd table lookup value % Do nearest-neighbor interpolation using 8 LSBs (treat as fractional remainder) rFracNT = numerictype(0,8,8); % fractional remainder data type idxFrac8LSBs = reinterpretcast(bitsliceget(idxUFIX16,8,1), rFracNT); rFraction = idxFrac8LSBs; % fractional value for linear interpolation temp = rFraction * (lutValAbove - lutValBelow); rslt = lutValBelow + temp;
Вот график результатов алгоритма:
x_vals = 0:(pi/128):(pi/4); xIdxLo = zeroBasedIdx - 1; xIdxHi = zeroBasedIdx + 4; figure; hold on; axis([x_vals(xIdxLo) x_vals(xIdxHi) 0.25 0.65]); plot(x_vals(xIdxLo:xIdxHi), double(FI_SIN_LUT(xIdxLo:xIdxHi)), 'b^--'); plot([x_vals(zeroBasedIdx+1) x_vals(zeroBasedIdx+2)], ... [lutValBelow lutValAbove], 'k.'); % Closest values plot(0.425, double(rslt), 'r*'); % Interpolated fixed-point result plot(0.425, sin(0.425), 'gs'); % Double precision reference result xlabel('X'); ylabel('SIN(X)'); lut_val_str = 'Fixed-point lookup table values'; near_str = 'Two closest fixed-point LUT values'; interp_str = 'Interpolated fixed-point result'; ref_str = 'Double precision reference value'; legend(lut_val_str, near_str, interp_str, ref_str); title('Fixed-Point Designer Lookup Table Based SIN with Linear Interpolation', ... 'fontsize',12,'fontweight','b');
SIN
Создайте 1 024 точки между [-2*pi, 2*pi),
stepSize = pi/256; thRadDbl = (-2*pi):stepSize:(2*pi - stepSize); % double precision floating-point thRadFxp = sfi(thRadDbl, 12); % signed, 12-bit fixed-point inputs
Сравните SIN фиксированной точки по сравнению с результатами SIN с двойной точностью
fxpSinTh = sin(thRadFxp); % fixed-point results sinThRef = sin(double(thRadFxp)); % reference results errSinRef = sinThRef - double(fxpSinTh); figure; hold on; axis([-2*pi 2*pi -1.25 1.25]); plot(thRadFxp, sinThRef, 'b'); plot(thRadFxp, fxpSinTh, 'g'); plot(thRadFxp, errSinRef, 'r'); ylabel('sin(\Theta)','fontsize',12,'fontweight','b'); set(gca,'XTick',-2*pi:pi/2:2*pi); set(gca,'XTickLabel',... {'-2*pi', '-3*pi/2', '-pi', '-pi/2', ... '0', 'pi/2', 'pi', '3*pi/2','2*pi'}); set(gca,'YTick',-1:0.5:1); set(gca,'YTickLabel',{'-1.0','-0.5','0','0.5','1.0'}); ref_str = 'Reference: sin(double(\Theta))'; fxp_str = sprintf('16-bit Fixed-Point SIN with 12-bit Inputs'); err_str = sprintf('Error (max = %f)', max(abs(errSinRef))); legend(ref_str, fxp_str, err_str); title(fxp_str,'fontsize',12,'fontweight','b');
Вычислите ошибку LSB
figure; fracLen = fxpSinTh.FractionLength; plot(thRadFxp, abs(errSinRef) * pow2(fracLen)); set(gca,'XTick',-2*pi:pi/2:2*pi); set(gca,'XTickLabel',... {'-2*pi', '-3*pi/2', '-pi', '-pi/2', ... '0', 'pi/2', 'pi', '3*pi/2','2*pi'}); ylabel(sprintf('LSB Error: 1 LSB = 2^{-%d}',fracLen),'fontsize',12,'fontweight','b'); title('LSB Error: 16-bit Fixed-Point SIN with 12-bit Inputs','fontsize',12,'fontweight','b'); axis([-2*pi 2*pi 0 8]);
Вычислите уровень шума
fft_mag = abs(fft(double(fxpSinTh))); max_mag = max(fft_mag); mag_db = 20*log10(fft_mag/max_mag); figure; hold on; plot(0:1023, mag_db); plot(0:1023, zeros(1,1024),'r--'); % Normalized peak (0 dB) plot(0:1023, -64.*ones(1,1024),'r--'); % Noise floor level (dB) ylabel('dB Magnitude','fontsize',12,'fontweight','b'); title('64 dB Noise Floor: 16-bit Fixed-Point SIN with 12-bit Inputs',... 'fontsize',12,'fontweight','b'); % axis([0 1023 -120 0]); full FFT axis([0 round(1024*(pi/8)) -100 10]); % zoom in set(gca,'XTick',[0 round(1024*pi/16) round(1024*pi/8)]); set(gca,'XTickLabel',{'0','pi/16','pi/8'});
Фиксированная точка алгоритм CORDIC требует следующих операций:
1 поиск по таблице на итерацию
2 сдвига на итерацию
3 сложения на итерацию
Упрощенный один алгоритм интерполяционной таблицы с линейной интерполяцией ближайшего соседа требует следующих операций:
2 поиска по таблице
1 умножение
2 сложения
В приложениях реального мира, выбирая алгоритм для вычислений тригонометрической функции фиксированной точки обычно зависит от требуемой точности, стоимости и аппаратных ограничений.
close all; % close all figure windows
Джек Э. Волдер, Тригонометрический Вычислительный Метод CORDIC, Транзакции IRE на Электронно-вычислительных машинах, Volume EC 8, сентябрь 1959, pp330-334.
Рэй Андрэка, обзор алгоритма CORDIC для основанных на FPGA компьютеров, Продолжений 1998 шестых международных симпозиумов ACM/SIGDA по Программируемым пользователем вентильным матрицам, 22-24 февраля 1998, pp191-200