Этот пример показывает, как использовать как основанные на 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