Forum: Digitale Signalverarbeitung / DSP / Machine Learning Realisierung eines IIR Butterworth Tiefpassfilters


von Steve (Gast)


Lesenswert?

Hallo,

ich möchte gerne einen IIR Butterworth Tiefpassfilter mit MATLAB 
entwerfen. Er soll 2. Ordnung und normierte Grenzfrequenz Wn = 0.1 
beispielsweise haben.
Leider kann ich nicht auf die integrierten Funktionen butter bzw filter 
zurückgreifen, da in der Firma die Toolboxen nicht installiert sind und 
ich somit nicht darauf zugreifen kann. Deshalb muss ich den Filter 
selbst programmieren, nur leider fehlt mir dazu komplett der Ansatz.

Die zu filternden Werte stehen in einem Array mit je nach Messung 
unterschiedlicher Anzahl der Messwerte.

Wie könnte man vorgehen? Gibt es Beispielscripte? Ich habe schon ne 
weile gesucht aber nix brauchbares gefunden.

Wär super wenn jemand helfen könnte!

Steve

von Timo (Gast)


Lesenswert?

Probiere Octave. Ist ein freier Matlab Clone, bei dem auch die meisten 
Toolboxen zumindest teilweise implementiert wurden.

von Detlef _. (detlef_a)


Lesenswert?

Anbei die Koeffizienten und der Algorithmus. Damit sollte sich das 
programmieren lassen.

Cheers
Detlef







> [fb,fa]=butter(2,0.1)

fb =

   0.02008336556421   0.04016673112842   0.02008336556421


fa =

   1.00000000000000  -1.56101807580072   0.64135153805756

>> help filter

 FILTER One-dimensional digital filter.
    Y = FILTER(B,A,X) filters the data in vector X with the
    filter described by vectors A and B to create the filtered
    data Y.  The filter is a "Direct Form II Transposed"
    implementation of the standard difference equation:

    a(1)*y(n) = b(1)*x(n) + b(2)*x(n-1) + ... + b(nb+1)*x(n-nb)
                          - a(2)*y(n-1) - ... - a(na+1)*y(n-na)

    If a(1) is not equal to 1, FILTER normalizes the filter
    coefficients by a(1).

    FILTER always operates along the first non-singleton dimension,
    namely dimension 1 for column vectors and non-trivial matrices,
    and dimension 2 for row vectors.

    [Y,Zf] = FILTER(B,A,X,Zi) gives access to initial and final
    conditions, Zi and Zf, of the delays.  Zi is a vector of length
    MAX(LENGTH(A),LENGTH(B))-1 or an array of such vectors, one for
    each column of X.

    FILTER(B,A,X,[],DIM) or FILTER(B,A,X,Zi,DIM) operates along the
    dimension DIM.

von Steve (Gast)


Lesenswert?

Danke für eure Anregungen. Jedoch hab ich mich oben wohl etwas 
unglücklich ausgedrückt. Die Ordnung 2 soll immer fix sein, jedoch 
sollte man über ein GUI die Grenzfrequenz verändern können, also nicht 
fest 0.1 sein, sondern variabel.
Die Funktion "filter" funktioniert, alles was ich also benötige sind die 
Filterkoeffizienten. Die müssten doch irgendwie über eine 
Differenzengleichung zu bestimmen sein oder? Nur wie lässt sich das 
realisieren?

Sorry nochmal für die unkorrekte Ausdrucksweise oben!

Steve

von nop(); (Gast)


Lesenswert?

Ja, sicher. Etwas Mathematik ist noetig um neue Koeffizienten zu 
rechnen. Nun ist die Frage ob und wie sich das auf dem Prozessor zu 
implementieren lohnt. Die Moeglichkeiten sind :
1)analytischer ansatz, berechnung der koeffizienten und laden.
2)fuer einen Bereich ein paar Werte rechnen, und in einem ROM speichern

von Detlef _. (detlef_a)


Angehängte Dateien:

Lesenswert?

Hallo Steve,

die Pole eines Butterworth der Ordnung zwei liegen ziemlich gut auf ner 
Parabel 4.Ordnung.

Deine Koeffizienten b des Filters lauten ja immer gleich, nämlich [1 2 
1].
Die Rückkoppelkoeffizienten kannst Du hinreichend genau folgendermaßen 
berechnen:

r sei die Grenzfrequenz, also z.B. 0.1

(1) y= -0.06680920400991*r^4-0.34511218278966*r^2+0.41379906532123
(2) a= (z-r+i*y)*(z-r-i*y), i=sqrt(-1)

a ist nen reelles Polynom in z mit Deinen gesuchten Koeffizienten als 
Faktoren von z^1 und z^0.

Das sollte für ne Wald- und Wiesenaufgabe eigentlich so hinhauen. Anbei 
noch das Matlab script, mit dem ich die Herleitung von (1) gebastelt 
habe und nen Bild des Abstandes der wahren und der angenäherten 
pole-locations.

Cheers
Detlef

clear
w=linspace(0,1,1000);
w=w(2:end-1);
rr=[];
for(k=1:length(w))
    [fb,fa]=butter(2,w(k));
    r=roots(fa);
    rr=[rr r(1)];
end;
x=real(rr);
y=imag(rr);
x=x(:);y=y(:);
M=[x.*x.*x.*x  x.*x.*x x.*x x ones(length(x),1) ];
coff=inv(M'*M)*M'*y;
%plot(,y,'r.',x,M*coff,'b.')
plot(x,y-M*coff,'r.')
return

von Steve (Gast)


Lesenswert?

Hallo Detlev,

besten Dank für deinen Post oben! Jetzt kommen wir der Sache doch schon 
gewaltig näher. Super!

Ich hätte jedoch noch eine Frage:

Bei der Berechnung der Koeffizienten unter (2) steig ich nicht ganz 
durch. Wie kann man das in Matlab implementieren, dass man die 
Koeffizienten bekommt? Multipliziere ich aus, ergibt sich für mich

a = z^2 - 2r*z + (r^2 + y^2)

Angenommen die Grenzfrequenz r sei 0.1 ergibt sich für y unter (1):
y = 0.4103

--> Einsetzten in a:
a = z^2 - 0.2*z + 0.1784

Nach dieser Berechnung lauten die Koeffizienten für A also A= [1 -0.2 
0.1784] oder?

Berechnet man die mit MATLAB und der Funktion Butter ergibt sich jedoch 
für A=[1 -1.5610 0.6414].

Es sind also deutliche Abweichungen. Dementsprechend funktioniert die 
Filterung auch nicht so wie das sein sollte.

Steh ich grad voll auf dem Schlauch, hab ich mich verrechnet oder was 
könnte da los sein?

Vielen Dank nochmals für dein super script!

Steve


von Detlef _. (detlef_a)


Lesenswert?

Ähm, shit, ich hab nen fetten Fehler gemacht. War wohl doch schon zu 
spät, gestern. Da arbeite ich nochmal nach.

Cheers
Detlef

von Detlef _. (detlef_a)


Lesenswert?

Das muß so:

r sei (1-2*Grenzfrequenz).

r=(1-2*0.1)=0.8
y=-0.06680920400991*0.8^4-0.34511218278966*0.8^2+0.41379906532123=
0.16556221837339

a=z^2-1.6*z+0.8^2+y^2
=a-1.6z+0.66741084815272

Matlab sagt -1.56 und 0.6414

Ist ein wenig schlechter approximiert, aber vielleicht reichts ja von 
der Genauigkeit.

Cheers
Detlef

von Steve (Gast)


Lesenswert?

Dankeschön!

Aber wie kommst du denn auf die Koeffizienten -1.56 und 0.6414? Die 
roots von a=z^2-1.6*z+0.8^2+0.16556^2 liegen bei mir laut MATLAB im 
konjugiert komplexen Bereich von

0.8000 + 0.1655i und
0.8000 - 0.1655i.

Sind somit nicht reell wie eigentlich vermutet. An was kann das noch 
liegen?

von Detlef _. (detlef_a)


Lesenswert?

Matlab sagt für butter(2,0.1) als Rückkoppelkoeffizienten des Filters 
-1.56 und 0.6414. Die Wurzeln des Polynoms (die Filterpole) liegen bei 
0.7805 +/- 0.1793 i. Die Wurzel treten immer konjugiert komplex auf, 
sonst werden die Koeffizienten ja nicht reel.

Meine vorgeschlagene Näherung für die Koeffizienten lautet -1.6 und 
0.667. Die Pole des Filters liegen, wie du richtig schreibst, also bei 
0.8 +/- 0.1655 i .

Habe ich Deine Frage damit richtig verstanden?

Cheers
Detlef

von Steve (Gast)


Lesenswert?

Hallo Detlev,

alles klar. Vielen vielen Dank für Deine Hilfe! Das ganze klappt so 
jetzt eigentlich ganz gut! Würdest du eine relativ einfache Möglichkeit 
sehen die approximierten Koeffizienten noch etwas genauer an die 
wirklichen ranzubringen? Ist dann auch meine letzte Frage und dann hast 
du Ruhe von mir;-) Hast mir eh schon so viel weitergeholfen!

Steve

von Detlef _. (detlef_a)


Lesenswert?

Ja, da kann man noch besser fitten. Die Imaginärteile Deines konjugiert 
komplexen Polpaars beschreiben in Abhängigkeit von Deiner Grenzfrequenz, 
wie gesagt, ziemlich genau diese Parabel 4.Ordnung. Die Realteile des 
Polpaares errechne ich zu (1-2*Grenzfrequenz), also als Gerade. Das ist 
sehr grob, die Kurve hat noch son S drin. Da kann man auch noch was 
fitten, dann wird das besser.

Cheers
Detlef

von Steve (Gast)


Lesenswert?

Hi Detlev,

vielen Dank nochmals! Das klingt logisch, ich werd mal versuchen die 
Kurve besser zu approximieren. Aber dann hätte ich doch noch eine 
Frage;-) Wie war da dein Ansatz für die Parabel 4.Ordnung der 
Imaginärteile des Polpaares?
Die B-Koeffizienten haben wie schon bekannt die Abhängigkeiten 
B(1)=B(3)=B(2)/2. Sind also nicht immer konstant, sondern je nach 
Grenzfrequenz verschieden.
Für die Grenzfrequenz 0.1 z.B. B(1)=B(3)= 0.0201; B(2)=0.0402.
Wie berechnen sich diese Koeffizienten abhängig von der Grenzfrequenz?

Vielen vielen Dank

Steve

von Detlef _. (detlef_a)


Lesenswert?

>>Sind also nicht immer konstant, sondern je nach Grenzfrequenz verschieden.
Ja, aber nur bis auf eine multiplikative Konstante. Es gilt immer 
B(1)=B(3)=B(2)/2. Behältst Du diese Eigenschaft bei, kannst Du die B's 
mit ner beliebigen Konstante multiplizieren, ohne die 'Form' des 
Frequenzgangs zu ändern. Du schiebst den Frequenzgang nur nach 'oben' 
oder 'unten', änderst also die Verstärkung des Filters. Matlab wird die 
vermutlich so wählen, daß das Filter Verstärkung 1 bei DC hat.

>>Wie war da dein Ansatz für die Parabel 4.Ordnung der Imaginärteile des 
Polpaares?

Ich habe Matlab die Imaginärteile des Polpaares berechnen lassen. Das 
sah ungefähr aus wie eine Parabel. Dann habe ich diese Parabel least 
square mit nem Polynom 4.Ordnung gefittet, siehe das Matlab script.

Cheers
Detlef

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.