Этот пример показывает, как использовать как основанные на CORDIC, так и основанные на интерполяционной таблице алгоритмы, предоставленные Fixed-Point Designer™, чтобы аппроксимировать синус MATLAB ® (SIN
) и косинус (COS
) функции. Эффективные алгоритмы синуса и косинуса с фиксированной точкой имеют решающее значение для многих встраиваемых приложений, включая управление двигателем, навигацию, обработку сигналов и беспроводную связь.
Введение
The cordiccexp
, cordicsincos
, cordicsin
, и cordiccos
функции аппроксимируют MATLAB sin
и cos
функционирует с использованием основанного на CORDIC алгоритма. CORDIC - это аббревиатура для COordinate Rotation DIgital Computer. Алгоритм CORDIC на основе вращения Givens (см. [1,2]) является одним из наиболее аппаратно эффективных алгоритмов, потому что он требует только итерационных операций shift-add. Алгоритм CORDIC устраняет необходимость в явных умножителях и подходит для вычисления множества функций, таких как синус, косинус, арксин, арксозин, арктангенс, векторные величины, деление, квадратный корень, гиперболические и логарифмические функции.
Можно использовать режим вычисления вращения CORDIC для вычисления синуса и косинуса, а также полярно-декартовых операций преобразования. В этом режиме известны векторы величины и угол поворота, и компоненты координат (X-Y) вычисляются после поворота.
Расчет вращения CORDIC
Алгоритм CORDIC режима вращения начинается с инициализации углового аккумулятора с желаемым углом поворота. Затем решение о повороте в каждой итерации CORDIC выполняется способом, который уменьшает величину аккумулятора остаточного угла. Решение о повороте основывается на знаке остаточного угла в аккумуляторе угла после каждой итерации.
В режиме вращения уравнения CORDIC:
где если и в противном случае;
, и общее количество итераций.
Это обеспечивает следующий результат при подходах:
Где:
.
В режиме вращения алгоритм CORDIC ограничивается углами поворота между и. Чтобы поддерживать углы за пределами этой области значений, cordiccexp
, cordicsincos
, cordicsin
, и cordiccos
функции используют квадрантную коррекцию (включая возможное дополнительное отрицание) после завершения итераций CORDIC.
CORDICSINCOS
Код синуса и косинусаВведение
The cordicsincos
функция вычисляет синус и косинус входных углов в области значений [-2 * pi, 2 * pi) с помощью алгоритма CORDIC. Эта функция принимает угол (радианы) и количество итераций в качестве входных параметров. Функция возвращает приближения синуса и косинуса.
Выходы расчета CORDIC масштабируются усилением вращателя. Это усиление учитывается путем предварительного масштабирования начального постоянного значения.
Инициализация
The 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
Создайте 1024 точки между [-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
.
The 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
Использование интерполяционных таблицСуществует много подходов, основанных на интерполяционной таблице, которые могут использоваться для реализации синусоидальных и косинусоидных приближений с фиксированной точкой. Ниже приведен дешевый подход, основанный на одной реальной интерполяционной таблице и простой линейной интерполяции по ближайшему соседу.
The sin
и cos
методы fi
объект в Fixed-Point Designer аппроксимирует sin
MATLAB ® builtin с плавающей точкой и
cos
функций, используя подход, основанный на интерполяции интерполяционной таблицы с простой линейной интерполяцией между значениями по ближайшему соседу. Этот подход позволяет создать небольшую реальную интерполяционную таблицу и использует простую арифметику.
Использование одной реальной интерполяционной таблицы упрощает расчет индекса и общую арифметику, необходимую для достижения очень хорошей точности результатов. Эти упрощения обеспечивают относительно высокую скорость эффективности, а также относительно низкие требования к памяти.
SIN
и COS
РеализацияРазмер и точность интерполяционной таблицы
Двумя важными факторами проекта интерполяционной таблицы являются ее размер и ее точность. Невозможно создать таблицу для всех возможных входных значений. Также невозможно быть совершенно точным из-за квантования или значений интерполяционной таблицы.
В качестве компромисса, Fixed-Point Designer SIN
и COS
методы FI
использовать 8-битную интерполяционную таблицу как часть их реализации. 8-битная таблица имеет длину всего 256 элементов, поэтому она является маленькой и эффективной. Восемь биты также соответствуют размеру байта или слова на многих платформах. Используемая в сочетании с линейной интерполяцией и 16-битным выходом (значение интерполяционной таблицы) точность, интерполяционная таблица 8 бита обеспечивает как очень хорошую точность, так и эффективность.
Инициализация постоянной SIN
Значения интерполяционной таблицы
Для простоты реализации, однородности значений таблицы и скорости используется полный таблицу синевав. Во-первых, четверть волны SIN
функция дискретизируется с 64 равномерными интервалами в области значений [0, pi/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);
Обзор реализации алгоритма
Реализация sin
Fixed-Point Designer и
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 самых значимых битов (MSB) этого полномасштабного беззнакового 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 наименее значимых битов (LSB) 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
Создайте 1024 точки между [-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
Jack E. Volder, The CORDIC Trigonometric Computing Technique, IRE Transactions on Electronic Computers, Volume EC-8, September 1959, pp330-334.
Рэй Андрака, Обзор алгоритма CORDIC для компьютеров на основе FPGA, Труды 1998 ACM/SIGDA шестой международный симпозиум по программируемым массивам ворот для поля, 22-24, 1998 февраля, pp191-200