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

Вычисление синуса и косинуса Используя алгоритм 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 (режим вращения) фрагмент ядра можно следующим образом (для случая скалярного xY, и 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

Ссылки