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
Probiere Octave. Ist ein freier Matlab Clone, bei dem auch die meisten Toolboxen zumindest teilweise implementiert wurden.
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.
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
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
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
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
Ähm, shit, ich hab nen fetten Fehler gemacht. War wohl doch schon zu spät, gestern. Da arbeite ich nochmal nach. Cheers Detlef
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
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?
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
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
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
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
>>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
Mit Google-Account einloggen
Noch kein Account? Hier anmelden.