mikrocontroller.net

Forum: Digitale Signalverarbeitung / DSP IIR Filter - Polynom


Autor: Tim B. (phantomias)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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.

Autor: Christoph Kessler (db1uq) (christoph_kessler)
Datum:

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

http://www.mikrocontroller.net/articles/Digitalfil...

Autor: Tim B. (phantomias)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
Hat schon geholfen, vielen Dank für deine Hilfe!

Autor: Michael Schlittenbauer (de1rush)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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

Autor: Frank Meier (aktenasche)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
inwiefern ist das ergebnis anders?

Autor: Michael Schlittenbauer (de1rush)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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

Autor: Detlef _a (detlef_a)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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

Autor: Michael Schlittenbauer (de1rush)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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?

Autor: Michael Schlittenbauer (de1rush)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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

Autor: Detlef _a (detlef_a)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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

Antwort schreiben

Die Angabe einer E-Mail-Adresse ist freiwillig. Wenn Sie automatisch per E-Mail über Antworten auf Ihren Beitrag informiert werden möchten, melden Sie sich bitte an.

Wichtige Regeln - erst lesen, dann posten!

  • Groß- und Kleinschreibung verwenden
  • Längeren Sourcecode nicht im Text einfügen, sondern als Dateianhang

Formatierung (mehr Informationen...)

  • [c]C-Code[/c]
  • [avrasm]AVR-Assembler-Code[/avrasm]
  • [code]Code in anderen Sprachen, ASCII-Zeichnungen[/code]
  • [math]Formel in LaTeX-Syntax[/math]
  • [[Titel]] - Link zu Artikel
  • Verweis auf anderen Beitrag einfügen: Rechtsklick auf Beitragstitel,
    "Adresse kopieren", und in den Text einfügen




Bild automatisch verkleinern, falls nötig
Bitte das JPG-Format nur für Fotos und Scans verwenden!
Zeichnungen und Screenshots im PNG- oder
GIF-Format hochladen. Siehe Bildformate.
Hinweis: der ursprüngliche Beitrag ist mehr als 6 Monate alt.
Bitte hier nur auf die ursprüngliche Frage antworten,
für neue Fragen einen neuen Beitrag erstellen.

Mit dem Abschicken bestätigst du, die Nutzungsbedingungen anzuerkennen.