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
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
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
