Разработайте алгоритмы фиксированной точки

В этом примере показано, как разработать и проверить простой алгоритм фиксированной точки.

Простой пример разработки алгоритмов

Этот пример показывает разработку и верификацию простого алгоритма фильтра фиксированной точки. Мы выполним следующие шаги:

1) Реализуйте алгоритм фильтра второго порядка и симулируйте в с двойной точностью, с плавающей точкой.

2) Оснастите код, чтобы визуализировать динамический диапазон выхода и состояния.

3) Преобразуйте алгоритм в фиксированную точку путем изменения типа данных переменных - сам алгоритм не изменяется.

4) Сравните и постройте фиксированную точку и результаты с плавающей точкой.

Определения переменной с плавающей точкой

Мы разрабатываем наш алгоритм в с двойной точностью, с плавающей точкой. Мы будем использовать фильтр lowpass второго порядка, чтобы удалить высокие частоты во входном сигнале.

b = [ 0.25 0.5      0.25    ]; % Numerator coefficients
a = [ 1    0.09375  0.28125 ]; % Denominator coefficients
% Random input that has both high and low frequencies.
s = rng; rng(0,'v5uniform');
x = randn(1000,1);
rng(s); % restore RNG state
% Pre-allocate the output and state for speed.
y = zeros(size(x));
z = [0;0];

Независимый от типа данных алгоритм

Это - фильтр второго порядка, который реализует стандартное разностное уравнение:

y(n) = b(1)*x(n) + b(2)*x(n-1) + b(3)*x(n-2) - a(2)*y(n-1) - a(3)*y(n-2)

for k=1:length(x)
    y(k) =  b(1)*x(k) + z(1);
    z(1) = (b(2)*x(k) + z(2)) - a(2)*y(k);
    z(2) =  b(3)*x(k)         - a(3)*y(k);
end

% Save the Floating-Point Result
ydouble = y;

Визуализируйте динамический диапазон

Для того, чтобы преобразовать в фиксированную точку, мы должны знать область значений переменных. В зависимости от сложности алгоритма эта задача может быть простой или довольно сложной. В этом примере известна область значений входного значения, так выбор соответствующего типа данных с фиксированной точкой прост. Мы сконцентрируемся на выходе (y) и состояния (z), поскольку их область значений неизвестна. Чтобы просмотреть динамический диапазон выхода и состояний, мы изменим код немного, чтобы оснастить его. Мы создадим два объекта NumericTypeScope и просмотрим динамический диапазон выхода (y) и состояния (z) одновременно.

Инструмент код с плавающей точкой

% Reset states
z = [0;0];

hscope1 = NumericTypeScope;
hscope2 = NumericTypeScope;
for k=1:length(x)
    y(k) =  b(1)*x(k) + z(1);
    z(1) = (b(2)*x(k) + z(2)) - a(2)*y(k);
    z(2) =  b(3)*x(k)         - a(3)*y(k);
    % process the data and update the visual.
    step(hscope1,z);
end
step(hscope2,y);

Анализируйте информацию в осциллографе

Давайте сначала анализировать информацию, отображенную для переменной z (состояние). От гистограммы мы видим, что динамический диапазон находится между (].

По умолчанию осциллограф использует размер слова 16 битов с нулевым терпимым переполнением. Это приводит к типу данных numerictype (верный, 16, 14), поскольку нам нужны по крайней мере 2 целочисленных бита, чтобы избежать переполнения. Можно получить больше информации о статистических данных от Входных данных и Получившихся панелей Типа. От панели Входных данных мы видим, что данные имеют и положительные и отрицательные величины и следовательно количество со знаком, которое отражается в предложенном numerictype. Кроме того, максимальное значение данных 1.51, который может быть представлен предложенным типом.

Затем давайте посмотрим на переменную y (выход). Из графика гистограммы мы видим, что динамический диапазон находится между (].

Определения переменной фиксированной точки

Мы преобразуем переменные в фиксированную точку и запускаем алгоритм снова. Мы включим логгирование, чтобы видеть переполнение и потери значимости, введенные выбранными типами данных.

% Turn on logging to see overflows/underflows.
FIPREF_STATE = get(fipref);
reset(fipref)
fp = fipref;
default_loggingmode = fp.LoggingMode;
fp.LoggingMode = 'On';
% Capture the present state of and reset the global fimath to the factory
% settings.
globalFimathAtStart = fimath;
resetglobalfimath;
% Define the fixed-point types for the variables in the below format:
%   fi(Data, Signed, WordLength, FractionLength)
b = fi(b, 1, 8, 6);
a = fi(a, 1, 8, 6);

x = fi(x, 1, 16, 13);
y = fi(zeros(size(x)), 1, 16, 13);
z = fi([0;0], 1, 16, 14);

Алгоритм совпадающего типа данных

for k=1:length(x)
    y(k) =  b(1)*x(k) + z(1);
    z(1) = (b(2)*x(k) + z(2)) - a(2)*y(k);
    z(2) =  b(3)*x(k)         - a(3)*y(k);
end
% Reset the logging mode.
fp.LoggingMode = default_loggingmode;

В этом примере мы переопределили переменные фиксированной точки с теми же именами как с плавающей точкой так, чтобы мы могли встроить код алгоритма для ясности. Однако это - лучшая практика, чтобы заключить код алгоритма в функцию файла MATLAB®, которая могла быть вызвана или с переменными или с фиксированной точки с плавающей точкой. Смотрите filimitcycledemo.m для примера записи и использования агностического типом данных алгоритма.

Сравните и постройте результаты и фиксированной точки с плавающей точкой

Мы теперь построим ответ величины результатов и фиксированной точки с плавающей точкой и ответ фильтра, чтобы видеть, ведет ли фильтр себя как ожидалось, когда это преобразовано в фиксированную точку.

n = length(x);
f = linspace(0,0.5,n/2);
x_response = 20*log10(abs(fft(double(x))));
ydouble_response = 20*log10(abs(fft(ydouble)));
y_response = 20*log10(abs(fft(double(y))));
plot(f,x_response(1:n/2),'c-',...
    f,ydouble_response(1:n/2),'bo-',...
    f,y_response(1:n/2),'gs-');
ylabel('Magnitude in dB');
xlabel('Normalized Frequency');
legend('Input','Floating point output','Fixed point output','Location','Best');
title('Magnitude response of Floating-point and Fixed-point results');

h = fft(double(b),n)./fft(double(a),n);
h = h(1:end/2);
clf
hax = axes;
plot(hax,f,20*log10(abs(h)));
set(hax,'YLim',[-40 0]);
title('Magnitude response of the filter');
ylabel('Magnitude in dB')
xlabel('Frequency');

Заметьте, что высокие частоты во входном сигнале ослабляются фильтром lowpass, который является ожидаемым поведением.

Постройте ошибку

clf
n = (0:length(y)-1)';
e = double(lsb(y));
plot(n,double(y)-ydouble,'.-r', ...
     [n(1) n(end)],[e/2 e/2],'c', ...
     [n(1) n(end)],[-e/2 -e/2],'c')
text(n(end),e/2,'+1/2 LSB','HorizontalAlignment','right','VerticalAlignment','bottom')
text(n(end),-e/2,'-1/2 LSB','HorizontalAlignment','right','VerticalAlignment','top')
xlabel('n (samples)'); ylabel('error')

Simulink®

Если у вас есть Simulink® и Fixed-Point Designer™, можно запустить эту модель, которая является эквивалентом алгоритма выше. Выход, y_sim является переменной фиксированной точки, равной переменной y вычисленный выше в коде MATLAB.

Как в коде MATLAB, параметры фиксированной точки в блоках могут быть изменены, чтобы совпадать с фактической системой; они собирались совпадать с кодом MATLAB в примере выше. Дважды кликните на блоках, чтобы видеть настройки.

if fidemo.hasSimulinkLicense

    % Set up the From Workspace variable
    x_sim.time = n;
    x_sim.signals.values = x;
    x_sim.signals.dimensions = 1;

    % Run the simulation
    out_sim = sim('fitdf2filter_demo', 'SaveOutput', 'on', ...
        'SrcWorkspace', 'current');

    % Open the model
    fitdf2filter_demo

    % Verify that the Simulink results are the same as the MATLAB file
    isequal(y, out_sim.get('y_sim'))

end

ans =

  logical

   1

Предположения, Сделанные для этого Примера

Для того, чтобы упростить пример, мы взяли математические параметры по умолчанию: раунд-к-самому-близкому, насыщайте на переполнении, продуктах полной точности и суммах. Мы можем изменить все эти параметры, чтобы совпадать с фактической системой.

Настройки были выбраны в качестве начальной точки в разработке алгоритмов. Сохраните копию этого файла MATLAB, начните проигрывать с параметрами и смотрите, какие влияния они оказывают на выход. Как алгоритм ведет себя с различным входом? Смотрите справку для fi, fimath, и numerictype для получения информации о том, как установить другие параметры, такие как округление режима и режима переполнения.

close all force;
bdclose all;
% Reset the global fimath
globalfimath(globalFimathAtStart);
fipref(FIPREF_STATE);

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