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

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

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 бит с нулевым допустимым переполнением. Это приводит к типу данных числового типа (true, 16, 14), поскольку нам нужно по крайней мере 2 целочисленных бита, чтобы избежать переполнения. Дополнительную информацию о статистических данных можно получить на панелях «Входные данные» и «Полученный тип». На панели Входных данных мы видим, что данные имеют как положительные, так и отрицательные значения и, следовательно, количество со знаком, которое отражено в предлагаемом типе числа. Кроме того, максимальное значение данных составляет 1,51, которое может быть представлено предлагаемым типом.

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

По умолчанию в возможностях используется размер слова 16 бит с нулевым допустимым переполнением. Это приводит к типу данных числового типа (true, 16, 14), так как нам нужно по крайней мере 2 целочисленных бита, чтобы избежать переполнения. С помощью этого предлагаемого типа вы не видите переполнений или подтекстов.

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

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

% 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 для примера записи и использования алгоритма типа данных-agnostic.

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

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

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

Если у вас есть Designer™ Simulink ® и Fixed-Point, можно запустить эту модель, которая является эквивалентной алгоритму выше. Выход, 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 и числовому типу для получения информации о том, как задать другие параметры, такие как режим округления и режим переполнения.

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