Цифровой фильтр в ассемблере, оптимизация данных Опции

[IGK]

Мне не удалось достичь необходимого динамического диапазона фильтра для следующих условий:
ФВЧ2, Fc=2Гц, Fsmp=606 Гц, Uin=1.14, аккумулятор 32 бит, ассемблер. Для ц все в порядке...

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

1.Для входных значений >1000 все в порядке, если меньше - "уезжает" вверх частота среза,
и для значений Uin<50 выход фильтра =0 (для контрольной частоты 3Гц д.быть 0.914 от 1.0).

2.Требуется повышение точности вычислений. 16 бит данных и 32 бита аккумулятора недостаточно. Повышение точности коэф.ничего не дает.

3.Возможно, поможет переход на direct-form II или отказ от стандартной формы вычислений (параллельная)?

Я решил перейти на арифметику 24данные*16коэф=40аккумулятор.
Как это сделать правильно? Я просто масштабировал на 2^8 входные данные - получил 24 бита. В таком виде и храню промежуточные значения.
Результат удовлетворительный. Для компенсации увеличения размера программы можно вычисления делать в каскадной форме, не обращая внимания на нормирующий коэффициент (только для этого фильтра).

Если не трудно - посмотрите свежим взглядом на модель, не ошибся ли где. А то, может, гоняю третий день ошибку и думаю, что все в порядке...

Пояснения к модели:
1.Вид подогнан к трансляции на ассемблер.
2.Для анализа использовал Фильтр Солюшн - там удобнее, имхо, чем в Матлабе, проверять коэффициенты и получать графики.
3.Если d=2^0 - 16*16=32, ecли d=2^8 - 24*16=40.
4.График 1: зеленым входной сигнал ампл.Asig, синим выходной.
5.График 2 показывает использование аккумулятора во время вычислений и нормирован к 2^39 (со знаком будет 2^40). Если больше 1.0 -переполнение.
6.График 3 показывает поток приведенных данных и нормирован к 2^23 (со знаком будет 2^24).

Вот м-файл.

-------------------------------------------------------------------------------------
clear
clc
clf
%Параметры сигнала ===================================================
Nsmp=3000; %Кол-во выборок
Fsmp=606; %Частота дискретизации Гц
Asig=16383; %Амплитуда сигнала
Fsig=3; %Частота сигнала
m=2^14; %Масштабирование коэффициентов
d=2^8; %Масштабирование данных 16->24 bit

%Входной сигнал 16bit ------------------------------------------------
k=1:Nsmp;
x(k) = fix(Asig*sin((k*Fsig*2*pi)/(Fsmp)));

%Коэффициенты ========================================================

%-- HPF2 2Hz stand.---------------------------------------------------
a0=fix(0.9854440*m); a1=-fix(1.970888*m); a2=fix(0.9854440*m);
b0=1*m; b1=-fix(1.970676*m); b2=fix(0.9710999*m);
%-- HPF2 2Hz cascade -------------------------------------------------
%a0=1*m; a1=-2*m; a2=1*m;
%b0=1*m; b1=-fix(1.970676*m); b2=fix(0.9710999*m);

%=== Узлы ============================================================
y0=0; y1=0; y2=0;
x0=0; x1=0; x2=0;

%=== Фильтр DF I =====================================================
for k = 1:1:Nsmp;
x2=x1; acc=x2*a2; ta1(k)=acc/(2^31*d); %r
x1=x0; acc=acc+x1*a1; ta2(k)=acc/(2^31*d); %g
x0=x(k)*d; acc=acc+x0*a0; ta3(k)=acc/(2^31*d); %b
y2=y1; acc=acc-y2*b2; ta4(k)=acc/(2^31*d); %m
y1=y0; acc=acc-y1*b1;
y0=fix(acc/m); ty1(k)=y0/(2^15*d);
y(k)=fix(y0/d); %Выход 16 бит
end;

%Графики =============================================================
k = 500:1:Nsmp; %Диапазон вывода
subplot(3,1,1);
hold on; grid on;
plot(k,x(k),'g'); %Входной сигнал
plot(k,y(k),'b');

%Тестовые сигналы ----------------------------------------------------
subplot(3,1,2);
hold on; grid on; %Контроль аккумулятора
plot(k,ta1(k),'r');
plot(k,ta2(k),'g');
plot(k,ta3(k),'b');
plot(k,ta4(k),'m');
ylabel('ACC / 2^3^1 / 2^3^9');

subplot(3,1,3);
hold on; grid on;
plot(k,ty1(k),'b'); %Промежуточные значения Y
ylabel('Y0 / 2^1^5 / 2^2^3');