Forum: Digitale Signalverarbeitung / DSP / Machine Learning IIR Filter in C implementieren


von Christian (Gast)


Angehängte Dateien:

Lesenswert?

Hallo,

will einen IIR Filter programmieren und eventuell später einmal auf 
einem Controller/DSP realisieren.

Bei dem Filter soll es sich um einen Tiefpassfilter handeln, der mir 
alle Frequenzen von 0 bis ca. 50HZ "durchlässt" und alles über 50Hz 
stark dämpft.
Um einen scharfen Knick hinzubekommen dachte ich an eine Art 
Tschebyscheff Filter.

Auf der Suche nach einem geeigneten Algorithmus bin ich auf dieses tolle 
Forum gestossen. Insbesondere auf dieses Topic

Beitrag "IIR Filter mit Matlab"

Dort wird folgender Algorithmus vorgestellt:

double z1=0.0,z2=0.0,r,y;
double x[]={45., 49. ,51., 45.};
double b0=0.25,b1=0.5,b2=0.25;
double a1= -1.88859 , a2=0.894476;

int k;
z1=0;z2=0;
for(k=0;k<4;k++)
{
    r=a1*z1+a2*z2;
    y=b0*(x[k]-r)+b1*z1+b2*z2;
    z2=z1;
    z1=x[k]-r;
    printf("%f \n",y);
}

Kann aber die Formel nicht ganz nachvollziehen
y=b0*(x[k]-r)+b1*z1+b2*z2;

Laut Blockschaltbild (Anhang) müsste sie doch eigentlich anders 
aussehen??
Wurde für die Berechnung ein anderes Blockschaltbild verwendet?
Wenn ja, wann muss man welches verwenden?

MFG

von Hans W. (Firma: Wilhelm.Consulting) (hans-)


Lesenswert?

bitte nicht erschlagen bin noch nicht ganz wach :)

das block-diagramm darf IMHO auch wirklich nur als solches verstanden 
werden...

nach dem FIR filter müsste folgende folge rauskommen:

y'=x[n]*b0+x[n-1]*b1

und nach der rückkopplung:

y=y'[n]+y'[n-1]*a1

ich find das so irgendwie übersichtlicher..

der source beschreibt wenns nach mir geht ein anders verhalten...

73

von Christian (Gast)


Angehängte Dateien:

Lesenswert?

Hi Hans,

bin genau Deiner Meinung.

Ich habe als Formel:

y(n) = x(n) * b0 + x(n-1) * b1 + y(n-1) * a1

Durch geschicktes Umstellen kann man dann noch ein Verzögerungsglied 
einsparen (Anhang)

Kann es denn überhaupt verschiedene Formeln geben, die Charakteristik 
eines Filters wird doch einzig und alleine durch die Dimensionierung der 
Parameter a1 bis an und b0 bis bm beschrieben.

Ich stelle mir das so vor, dass der Algorithmus (Blockschaltbild) immer 
gleich bleiben muss und lediglich die Koeffizienten geändert werden 
müssen je nachdem welche Frequenzen man haben möchte

MFG

von Detlef _. (detlef_a)


Lesenswert?

Für eine gegebene Differenzengleichung gibt es einige 
Implementationsmöglichkeiten, zwei hatte Christian ja schon gezeigt. Der 
Algorithmus berechnet die 2. kanonische Form entsprechend der zweiten 
geposteten Filterstruktur, allerdings 2. Ordnung, dargestellt ist ein 
Filter 1. Ordnung.

Stichworte sind da 'IIR' 'Direktform' (direkte Impelementation anhand 
der Differenzengleichung) oder 'kanonische Form' (minimale Anzahl der 
Speicher).

Google oder 
http://de.wikipedia.org/wiki/Filter_mit_unendlicher_Impulsantwort oder 
Arild Lacroix: 'Digitale Filter', Oldenbourg Verlag, sagt da mehr.

Cheers
Detlef

von Christian (Gast)


Lesenswert?

Hallo Detlef_a

danke für den interessanten Link.
Ähnliches hatte ich eben im Netz gefunden und einen Kaskadierte 
IIR-Filter 2. Ordnung (SOS) in C programmiert. Das erscheint mir alles 
nicht schwer.
Mein Plan war nun mit folgendem Link:

http://www.dsptutor.freeuk.com/IIRFilterDesign/IIRFilterDesign.html

mir die Koeffizienten ausrechnen zu lassen und mein Prog damit zu 
füttern.
So wie ich das auf Wiki. verstehe, muss ich die Koeffizienten aber erst 
umrechnen damit ich sie für den kaskadierten Filter nutzen kann, stimmt 
das so?
Sorry, hatte bisher mit Regelungstechnik nicht viel am Hut.

Wie schon erwähnt habe ich einen kaskadierten IIR-Filter 2. Ordnung 
(SOS) bereits programmiert. Gibt es noch eine Möglichkeit das Prog 
unkompliziert zu testen? Ich bräuchte einfach nur ein paar Zahlenwerte 
für die Koeffizienten und für den Eingang und könnte dann meine 
Ausgangswerte vergleichen. Das wäre denke die einfachste Möglichkeit.

MFG

von Detlef _. (detlef_a)


Lesenswert?

>>So wie ich das auf Wiki. verstehe, muss ich die Koeffizienten aber erst
>>umrechnen damit ich sie für den kaskadierten Filter nutzen kann, stimmt
>>das so?

Welche Ordnung hat das Filter denn, das Du implementieren willst? Falls 
es zweiter Ordnung ist braucht man für nen SOS nichts umzurechnen.

Die Filterübertragungsfkt. sind gebrochen rationale Polynome. Nen SOS 
realisiert Ordnung 2 im Zähler/Nenner. Falls Deine Übertragungsfkt. 
Polynome  höherer Ordnung im Zähler/Nenner hat, muß Du sie in Polynome 
2. Ordnung zerlegen.

Den Test machst Du so, wie in dem von Dir zitierten thread beschrieben: 
Impulsantwort per Hand oder mit Matlab/octave ausrechnen, dann kucken ob 
C dasgleiche rauskriegt. Der Code in dem thread funktioniert für 2. 
Ordnung auch super, der is nämlich von mir ;)  .

Cheers
Detlef

von Christian (Gast)


Lesenswert?

Hi Detlef_a

Vielen Dank Du hast mir schon sehr geholfen.

Die Ordnung meines Filters steht ehrlich geagt nicht fest. Ich wollte 
eben ein wenig experimentieren und später dann anhand der Ergebnisse und 
Rechenzeit abwägen für welche Ordnung ich mich entscheide.

Dein Code funktioniert in der Tat super, haben ihn eben mal getestet. 
Die Methode mit dem z finde ich sehr elegant :)

Erst mal Zu dem Code, den kann ich ja bestimmt mit dem Kuhschwanzprinzip 
um x-Ordnungen erweitern.
Stelle mir das so vor: (4. Ordnung)

int k;
z1=0;z2=0;z3=0;z4=0;
for(k=0;k<n_Eingangswerte;k++)
{
    r=a1*z1+a2*z2+a3*z3+a4*z4;
    y=b0*(x[k]-r)+b1*z1+b2*z2+b3*z3+b4*z4;
    z4=z3;
    z3=z2;
    z2=z1;
    z1=x[k]-r;
    printf("%f \n",y);
}

Aber wie kann ich mir die Koeffizienten der Ordnungen >2 für diese 
Methode (Blockschaltbild) ausrechnen lassen?
Matlab hätte ich zur Verfügung, kann auch bisschen programmieren aber 
mit den ganzen Filterbefehlen blicke ich nicht durch :(

MFG

von Christian (Gast)


Lesenswert?

Achso, noch ein kleiner Zusatz.

Dein Code beschreibt ja die Direct-Form 2 (DF2).
Diesen kann ich ja relativ einfach erweitern und meine Koeffizienten im 
Applet ausrechnen lassen und es sollte funktionieren.

Bei höheren Ordnungen soll man ja aufgrund der Stabilität auf den 
Kaskadierte IIR-Filter zurückgreifen, der quasi immer Filter 2. Ordnung 
hintereinander schaltet. Dort muss ich die Koeffizienten dann umrechnen. 
Hast Du da vielleicht ein Matlab Beispiel, das wäre hilfreich.

Ab welcher Ordnung ist es denn besser den Kaskadierten Filter zu 
benutzen?
Eigentlich kann ich den doch immer nehmen?

MFG

von Detlef _. (detlef_a)


Lesenswert?

Christian wrote:

> Ab welcher Ordnung ist es denn besser den Kaskadierten Filter zu
> benutzen?
> Eigentlich kann ich den doch immer nehmen?
>
Ja, kaskadierte Filter 2.Ordnung geht immer. Ist aber eher relevant für 
fixed point Berechnung, in float kannst Du auch höhere Ordnung direkt 
rechnen.

> Dort muss ich die Koeffizienten dann umrechnen.
> Hast Du da vielleicht ein Matlab Beispiel, das wäre hilfreich.
>

Nullstellenbestimmung eines Zähler/Nennerpolynoms, dann zwei 
Nullstellen/Polstellen zusammenfassen zum Polynom 2. Ordnung. Wenn 
konjugiert komplexe Nullstellen/Pole rauskommen, immer die kombinieren, 
sonst gibts keine reellen Koeffizienten.

Zaehlerpolynom= [1 2 3 4 5];
Nennerpolynom=[6 7 8 9 10];

rz=roots(Zaehlerpolynom);
rn=roots(Nennerpolynom);

zaehler_kaskade1=poly([rz(1) rz(2)])
nenner_kaskade1 =poly([rn(1) rn(2)])
zaehler_kaskade2=poly([rz(3) rz(4)])
nenner_kaskade2 =poly([rn(3) rn(4)])

Cheers
Detlef

Ach so, ersten post nicht gesehen:
>Die Methode mit dem z finde ich sehr elegant
Naja, das ist eher anschaulich. Im wahren Signalverarbeiterleben kopiert 
man die Inhalte der z nicht um, sondern läßt nen wrappenden pointer über 
nen array laufen. In der Regel sind die b's ja auch symmetrisch, das 
nutzt der Code nicht aus.

von Christian W. (christian_83)


Lesenswert?

Hi,

mein Prog. ist soweit fertig. Alles funktioniert mit theoretischen 
Teststreams.
Sobald ich Zeit habe werde ich mal einen Testaufbau machen und reale 
Werte mit nem Oszi aufnehmen und in meinem Prog. verarbeiten.

Irgendwann würde ich dann noch gerne einen passenden Controller oder DSP 
(abhängig von den ersten praktischen Tests) auswählen wollen der dann 
alles vollständig übernehmen soll.
Kennt einer ne gute Seite wo man ne Art Vergleichstabelle von 
Controllern/DSPs hat. Hatte mal ne Seite gefunden wo man sich die 
Parameter des persönlichen Traumcontrollers zusammenstellen konnte und 
der spuckte dann den passenden Typ heraus. Leider habe ich den Link 
nicht mehr :(

Eine Frage noch. Wie bereits erwähnt habe ich mir meine 
Filterkoeffizienten per Software ausspucken lassen. Die befinden sich in 
utopischen Größenordnungen wie 4.03335632290*10^-9. Wie wird sowas 
verarbeitet. Stelle mir die Multiplikation solcher Werte eher langsam 
vor. Bedarf es für diese Anwendung vielleicht sogar auf jeden Fall einen 
DSP der dafür ausgelegt ist??

MFG

von Mike C. (Firma: Arbeitslos) (mikecontroller)


Lesenswert?

Bei VisualDSP++ gibt es eine Funktion die IIR Filter für mehrere 
Ordnungen realiesiert. In Sektionen zu je 2.Ordnung:
Beitrag "IIR Filter mit Matlab"

Ich habe folgendes im Internet gefunden (S.4): 
http://users.etech.haw-hamburg.de/users/ITLabor/DV_Lab_R885/Laboraufgaben_Labtasks/DVP_I6/LAB_IIR_de.pdf

Matlab rechnet z.B. in Gleitkommazahlen. Für die Koeffizienten sollte 
man bei der Genauigkeit im Breich von float ode besser double bleiben. 
Also runden. Würdest du die Koeffizienten nicht genau angeben, könnte es 
sein das du damit ausserhalb des Einheitskreises kommst. Dann ist der 
Filter instabil. Der DSP kann ja nur diskret Werte im Einheitskreis 
ablegen.

Geeignet sind ADSP-BF5xx Reihe so wie ich das rausgelesen habe.

PS: Ja, ich weiß der threat ist uralt.

von Mike C. (Firma: Arbeitslos) (mikecontroller)


Lesenswert?

Kann man die koeffizienten auch von einem DSP berechnen lassen. Also ein 
variabler Tiefpassfilter ist in meinem Sinn. Geht das, oder zu 
aufwendig?

Beitrag #5095930 wurde vom Autor gelöscht.
Bitte melde dich an um einen Beitrag zu schreiben. Anmeldung ist kostenlos und dauert nur eine Minute.
Bestehender Account
Schon ein Account bei Google/GoogleMail? Keine Anmeldung erforderlich!
Mit Google-Account einloggen
Noch kein Account? Hier anmelden.