В этом примере показано, как разработать и проверить простой алгоритм фиксированной точки.
Этот пример показывает разработку и верификацию простого алгоритма фильтра фиксированной точки. Мы выполним следующие шаги:
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 (выход). Из графика гистограммы мы видим, что динамический диапазон находится между (
].
По умолчанию осциллограф использует размер слова 16 битов с нулевым терпимым переполнением. Это приводит к типу данных numerictype (верный, 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 для примера записи и использования агностического типом данных алгоритма.
Мы теперь построим ответ величины результатов и фиксированной точки с плавающей точкой и ответ фильтра, чтобы видеть, ведет ли фильтр себя как ожидалось, когда это преобразовано в фиксированную точку.
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® и 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);
