Документация

Вычисление синуса и косинуса Используя алгоритм CORDIC

Введение

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 выполняет выполняющие шаги инициализации:

Разумный выбор начальных значений позволяет алгоритму непосредственно вычислять и синус и косинус одновременно. После итераций эти начальные значения приводят к следующим выходным параметрам как подходы:

Разделяемая фиксированная точка и код ядра 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 не будет постоянным, потому что отличался бы.

Для очень больших значений алгоритм 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 требует следующих операций:

Упрощенный один алгоритм интерполяционной таблицы с линейной интерполяцией ближайшего соседа требует следующих операций:

В приложениях реального мира, выбирая алгоритм для вычислений тригонометрической функции фиксированной точки обычно зависит от требуемой точности, стоимости и аппаратных ограничений.

close all; % close all figure windows

Ссылки