IIR Filter - Polynom

[Tim B. (phantomias)]

Hallo.

Ich habe ein Digitalfilter (IIR), welcher eine Kaskade aus 2 Blöcken je 
zweiten Grads ist (konj. kompl. NS & Pole) + ein Block ersten Grads 
(reeller Pol, NS bei -1); Multiplikationen mit 1/-1 werden nicht 
gerechnet (NS auf Einheitskreis -> im Zählerpolynom zwei Koeffizienten = 
1).

Hier mal die Koeffizienten für die Cascades und den Block ersten Grads 
in gekürzter Form:
http://ud05_188.ud05.udmedia.de/spotlight/iir.png

Aber wie baue ich nun dieses Polynom zweiten Grads auf mittels der 
Koeffizienten?

Und wie kriege ich den Block des ersten Grads? Diese 3 Teile müssen dann 
ja miteinander multipliziert werden, so dass ein Polynom 5.ten Grads 
entsteht?

Gruss Tim.

[Christoph db1uq K. (christoph_kessler)]

schreib das doch mal so hin, das man es auch lesen kann in der üblichen 
Schreibweise:

$$H(z)=\frac{b_{2}*z^{-2}+b_{1}*z^{-1}+b_{0}}{a_{2}*z^{-2}+a_{1}*z^{-1}+1} = {\frac{b_{2}+b_{1}*z+b_{0}*z^{2}}{a_{2}+a_{1}*z+z^{2}}}$$

http://www.mikrocontroller.net/articles/Digitalfilter_mit_ATmega#Von_der_.C3.9Cbertragungsfunktion_zum_Programm

[Tim B. (phantomias)]

Hat schon geholfen, vielen Dank für deine Hilfe!

[Michael S. (de1rush)]

Hallo,

würde hierzu auch gerne eine Frage stellen:

ich hab mir IIR Filterkoeffizienten für einen Bandpass [300Hz; 3400Hz] 
mittels MATLAB berechnet. Der Filter funktioniert auch mit dem Befehl
[b,a] = butter(n,Wn, 'bandpass')
filter(b, a, y)
sehr gut. y ist der Audiodateneingang. Die FFT zeigt das gewünschte 
Ergebnis.

nun will ich das ganze in C Code portieren und wollte mir den Filter 
handprogrammiert verifizieren:


xc = y; (Koeffizienteninit X = Eingang)
x0 = 0;
x1 = 0;
x2 = 0;
x3 = 0;
x4 = 0;

yc = y; (init nur um gleiche Datenbreite wie y zu erhalten)
y1 = 0;
y2 = 0;
y3 = 0;
y4 = 0;

for i= 1:N

    x4 = x3;
    x3 = x2;
    x2 = x1;
    x1 = x0;
    x0 = xc(i);

    y4 = y3;
    y3 = y2;
    y2 = y1;
    y1 = yc(i);

    yc(i) = b(1) * x0 + b(2)  x1  b(3) * x2 + b(4) * x3 + b(5) * x4...
            - a(2) * y1 - a(3) * y2 - a(4) * y3 - a(5) * y4;

end

leider kommt hier nicht das gleiche Ergebnis wie beim filter() Befehl 
raus. Ich befürchte ich mach hier nen ganz trivialen Fehler :((

Könnt ihr bitte kurz drüberschaun?

Danke

[Frank M. (aktenasche)]

inwiefern ist das ergebnis anders?

[Michael S. (de1rush)]

also ich mach ne FFT über das Ausgangssignal und eine FFT über das 
Signal welches nach der Filter() Anweisung entsteht. Alle Frequenzen 
ausserhalb [300;3400]Hz wurden ausreichend gedämpft.
Wenn ich aber das Ausgangssignal per Code im Thread filtere was 
eigentlich das gleich Ergebnis erzeugen sollte, so werden Frequenzen 
überhalb 3400Hz klar verstärkt und die geforderten Frequenzen nicht 
gedämpft.
Ich hoffe das stammt nicht von Ungenauigkeiten?!?
Ich verwende ja die gleichen Filterkoeffizienten für filter() sowie für 
die Handrechnung.

grüße

[Detlef _. (detlef_a)]

Dises Zeile ist buggy
y1 = yc(i);

y1 ist Dein altes y(n-1), das ist in Deinem Code nicht richtig 
implementiert.  Debugging macht man nicht mit ner FFT, sondern man kuckt 
mit nem kurzen, einfachen Eingangsignal, zB Dirac, ob der eigene Code 
dasselbe liefert wie 'filter'.

Cheers
Detlef

[Michael S. (de1rush)]

Hey,

jo super danke für die Hilfe, habs ausgebessert und provisorisch durch
...
    y4 = y3;
    y3 = y2;
    y2 = y1;
    y1 = temp;

    yc(i) = b(1) * x0 + b(2)  x1  b(3) * x2 + b(4) * x3 + b(5) * x4...
            - a(2) * y1 - a(3) * y2 - a(4) * y3 - a(5) * y4;

    temp = yc(i);
...


ersetzt.

nun sollte in temp immer der y(n-1)te Wert stehen. Leider ist mein 
filter nun instabil :(
mit der Testfunktion setze ich mich noch auseinander! Danke für den Tip.

Seht ihr so aus der Hüfte noch nen Fehler?

[Michael S. (de1rush)]

Hallo,

ok das Problem hat sich gelöst, hab den Filter mit ner Sprungfunktion 
geprüft und nachgerechnet... leider hat sich bei x1 ein * statt ein - 
eingeschlichen :(

Blöder Fehler,

Aber Danke für die Hilfe!!

grüße

[Detlef _. (detlef_a)]

Geht doch, Glückwunsch. Zwei Dinge: Ich würde das Filter in einer 
cascadierten kanonischen Direktform rechnen: zweimal zweiter Ordnung 
hintereinander, dann auch nur zwei Zwischenspeicher pro Teilsystem, das 
ist Standard und wenn Du das auf nem DSP rechnen willst hat der seine 
Hardware vllt. so vorgegeben. Diese Umspeicherei y4=y3; etc. ist 
extremst uncool. Array und Pointer oder Indizes verwenden.

Cheers
Detlef