Простейший БИХ-фильтр, Не пойму как коэффиценты сделать целочисельными Опции

[condor]

Надо на 16-разрядном контроллере сделать ФНЧ.
Пытаюсь в МАТЛАБе смоделировать хотя бы фильтр второго порядка в виде канонического звена.
Коэффициенты получил с помощью fdatool.
С дробными коэффициентами все правильно работает, а как сделать их челочисельными?
Везде в литературе об этом вскользь упоминается как о само собой разумеющемся.
Если их просто умножить на 2^14, то из-за рекурсии при вычислении
w(n) = x(n) - a1*w(n-1) - a2*w(n-2)
w(n) быстро уходит в бесконечность. (На всякий случай: w - это промежуточное значение в каноническом звена; x - входной сигнал; ai - коэф.).
Подскажите, а то я сегодня туго соображаю

[condor]

Я все это делал.
Квантовал в fdatool.
Написано Stable: Yes

Quantized Numerator:
1
1.99591064453125
1
Quantized Denominator:
1
0.06439208984375
0.3143310546875
Quantized Gain:
0.3074951171875

С этими квантованными коэффициентами все работает отлично.
Но если их умножить на 2^14, то получим.
a1 = 0,0644*2^14 = 1055;
a2 = 0,3143*2^14 = 5150;
b1 = 1.99591*2^14 = 32701;
b2 = 1*2^14 = 16384;
sk = 0,3075*2^14 = 5038
Такие же коэффициенты генерятся в Сишном хедере.

w(n) = sk*x(n) - a1*w(n-1) - a2*w(n-2) ;
w(n) = 5038*x(n) - 1055*w(n-1) - 5150*w(n-2);
Начальные w я принял равными 0, а x(n) - какое-нибудь постоянное значение.
Теперь даже проделав вручную нескольно итераций видно, что w(n) нестабильно.
Думал что в формуле опечатка, но ведь с дробными кавантованными все работает, и даже АЧХ примерно совпадает.

Вот МАТЛАБовский код. Может я где-то ошибся, но кажется тут особо не где.

Nsamp = 2000;
wk = 0;
w1 = 1;
w2 = 1;
sk = 0.3075 * 2^14;
a1 = 0.0644 * 2^14;
a2 = 0.3143 * 2^14;
b0 = 1;
b1 = 1.996;
b2 = 1;

freq = 10;
t = 0:0.0005:1;
t (end) = [];

x = sin (2*pi*freq*t);

for n = 1:Nsamp
xn = x(n);
wk = sk*x(n) - a1*w1 - a2*w2;
y(n) = b0*wk + b1*w1 + b2*w2;
w2 = w1;
w1 = wk;
end

[Doka]

как-то так:

Код

Nsamp = 2000;
wk = 0;w1 = 1;w2 = 1;
sk = 0.3075 * 2^14;
a1 = 0.0644 * 2^14;
a2 = 0.3143 * 2^14;
b0 = 1;
b1 = 1.996;
b2 = 1;

freq = 10;
t = 0:0.0005:1;
t (end) = [];
x = sin (2*pi*freq*t);

NN=1:Nsamp;
for n = NN
    wk = sk*x(n) - a1*w1 - a2*w2;
    wk = wk *2^-14;      %    <----- добавил только это (остальное - косметика)
    y(n) = b0*wk + b1*w1 + b2*w2;
    w2 = w1;      w1 = wk;
end

plot (NN,x,NN,y);

upd: т.е. в базисе программирования МК: взятие старших разрядов от рекурсивной части (там где вы применили целочисленные коэфф.) - уж коли вы с целочисленной арифметикой работаете .
Си-компилятор не знает - целые/дробные у вас числа или какие другие
аппаратный умножитель - тоже: как трактовать результат мат.операций - на совести программиста..

Сообщение отредактировал Doka - Jul 26 2006, 16:09

[zhorro]

Общая формула фильтрации

y(k) = a0*x(k)+a1*x(k-1)+...+b1*y(k-1)+b2*y(k-2)+...

так вот, просто так взять да и умножить a0..aN на какое-то число все испортится, а вот если помножить и правую и левую часть на одно и тоже большое число (ну там на те же 2^14) то вроде как все останется не изменным.
Но! Тогда и левая часть помножится на 2^14 !

2^14*y(k) = 2^14*a0*x(k)+2^14*a1*x(k-1)+...+2^14*b1*y(k-1)+2^14*b2*y(k-2)+...

Если обозвать 2^14*aK = aK', 2^14*bK = bK' и взять целую часть то для вычисления выхода фильтра нужно проделать следующую процедуру

y(k) = (a0'*x(k)+a1'*x(k-1)+...+b1'*y(k-1)+b2'*y(k-2)+...) / 2^14

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