Вычисление синуса и косинуса с фиксированной точкой

Этот пример показывает, как использовать как основанные на CORDIC, так и основанные на интерполяционной таблице алгоритмы, предоставленные Fixed-Point Designer™, чтобы аппроксимировать синус MATLAB ® (SIN) и косинус (COS) функции. Эффективные алгоритмы синуса и косинуса с фиксированной точкой имеют решающее значение для многих встраиваемых приложений, включая управление двигателем, навигацию, обработку сигналов и беспроводную связь.

Вычисление синуса и косинуса с помощью алгоритма CORDIC

Введение

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:

$$ z_{i+1} = z_{i} - d_{i}*atan(2^{-i}) $$

$$ x_{i+1} = x_{i} - y_{i}*d_{i}*2^{-i} $$

$$ y_{i+1} = y_{i} + x_{i}*d_{i}*2^{-i} $$

где если$$ d_{i} = -1 $$ и $$ z_{i} < 0 $$$$ +1 $$в противном случае;

$$ i = 0, 1, ..., N-1 $$, и$$ N $$ общее количество итераций.

Это обеспечивает следующий результат при$$ N $$ подходах:$$ +\infty $$

$$ z_{N} = 0 $$

$$ x_{N} = A_{N}(x_{0}\cos{z_{0}} - y_{0}\sin{z_{0}}) $$

$$ y_{N} = A_{N}(y_{0}\cos{z_{0}} + x_{0}\sin{z_{0}}) $$

Где:

$$ A_{N} = \prod_{i=0}^{N-1}{\sqrt{1+2^{-2i}}} $$.

В режиме вращения алгоритм CORDIC ограничивается углами поворота между$$ -\pi/2 $$ и. $$ \pi/2 $$Чтобы поддерживать углы за пределами этой области значений, cordiccexp, cordicsincos, cordicsin, и cordiccos функции используют квадрантную коррекцию (включая возможное дополнительное отрицание) после завершения итераций CORDIC.

Понимание CORDICSINCOS Код синуса и косинуса

Введение

The cordicsincos функция вычисляет синус и косинус входных углов в области значений [-2 * pi, 2 * pi) с помощью алгоритма CORDIC. Эта функция принимает$$ \theta $$ угол (радианы) и количество итераций в качестве входных параметров. Функция возвращает приближения синуса и косинуса.

Выходы расчета CORDIC масштабируются усилением вращателя. Это усиление учитывается путем предварительного масштабирования начального$$ 1 / A_{N} $$ постоянного значения.

Инициализация

The cordicsincos функция выполняет следующие шаги инициализации:

  • Интерполяционная таблица входного угла inpLUT установлено в atan(2 .^ -(0:N-1)).

  • $$ z_{0} $$ задается значение$$ \theta $$ входного параметра.

  • $$ x_{0} $$ установлено в.$$ 1 / A_{N} $$

  • $$ y_{0} $$ установлено в нуль.

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

$$ x_{N} \approx cos(\theta) $$

$$ y_{N} \approx sin(\theta) $$

Общий код ядра 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$$ A_{n} $$ не будет постоянным, потому что$$ n $$ он будет варьироваться.

Для очень больших значений, $$ n $$алгоритм 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 Реализация

Размер и точность интерполяционной таблицы

Двумя важными факторами проекта интерполяционной таблицы являются ее размер и ее точность. Невозможно создать таблицу для всех возможных входных значений. $$ u $$Также невозможно быть совершенно точным из-за квантования$$ sin(u) $$ или $$ cos(u) $$значений интерполяционной таблицы.

В качестве компромисса, 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 объекты включают первое приведение входных входов угла с фиксированной точкой$$ u $$ (в радианах) к предварительно заданному типу данных в области значений [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

Ссылки

  1. Jack E. Volder, The CORDIC Trigonometric Computing Technique, IRE Transactions on Electronic Computers, Volume EC-8, September 1959, pp330-334.

  2. Рэй Андрака, Обзор алгоритма CORDIC для компьютеров на основе FPGA, Труды 1998 ACM/SIGDA шестой международный симпозиум по программируемым массивам ворот для поля, 22-24, 1998 февраля, pp191-200