Hallo Ich habe gerade Problem mit meinem IIR Filter, und zwar er ist nicht stabil Er hat einige Pole, die außerhalb der Einheitskreis und damit, instabil Ganz allgemeine Fragen: welche Möglichkeiten hat man, um diesen instabilen Filter zu stabilisieren ? Viele Grüße
Was meinst Du? Du kannst das instabile Filter so in einen Regelkreis einbauen, dass der insgesamt stabil wird. Ansonsten müssen die Pole schon im Einheitskreis bleiben, ansonsten wird das Filter instabil. Cheers Detlef
Detlef _a schrieb: > Was meinst Du? > Du kannst das instabile Filter so in einen Regelkreis einbauen, dass der > insgesamt stabil wird. Ansonsten müssen die Pole schon im Einheitskreis > bleiben, ansonsten wird das Filter instabil. > > Cheers > Detlef Kann man z.B die Übertragungsfunktion von diesem instabilen Filter zerlegen, in einer minimalphasigen System und ein Allpass. Dann ist der Filter stabil oder Gruß
Man kann ein nicht-Minimalphasensystem in ein Minimalphasensystem und einen Allpass zerlegen. Dabei geht es aber um Nullstellen ausserhalb des Einheitskreises, Du hast aber Polstellen ausserhalb. Cheers Detlef
Anh Le schrieb: > Ganz allgemeine Fragen: welche Möglichkeiten hat man, um diesen > instabilen Filter zu stabilisieren ? Filterkoeffizienten ändern ;-)
Detlef _a schrieb: > Man kann ein nicht-Minimalphasensystem in ein Minimalphasensystem und > einen Allpass zerlegen. Dabei geht es aber um Nullstellen ausserhalb des > Einheitskreises, Du hast aber Polstellen ausserhalb. > > Cheers > Detlef Danke für deine Antwort Konkret ist so folgendes Ich habe einen FIR Filter mit der Übertrafungsfunktion H(z) und dies soll ich invertieren als H'(z) = 1/ H(z) und ich soll jetzt H'(z) stabilisieren. H'(z) hat einige Pole, die außerhalb des Einheitskreises liegen. Kann ich meine obige Ansatz verwenden. Wäre sehr nett, wenn du kurz erklärst, warum kann ich oder warum kann ich nicht Viele Grüße
Die Nullen von H(z) werden ja zu Polen von H'(z) und zwar mit Kehrwert des Betrags der ursprünglichen Null. Vllt. trifft es ja die Dir vorgegebene spec. wenn du die Nullen von H(z), die innerhalb des Einheitskreises liegen, am Kreis spiegeltst, also Winkel beibehalten und Kehrwert des Betrags nehmen. Dann sind alle Nullen ausserhalb und fallen nach Transformation zu H'(z) in den Einheitskreis, das Ding wird also stabil. Cheers Detlef
Moin Und noch eine Frage: Kann man jedes zeitdiskrete System zerlegen oder nur stabiles nicht minimalphasiges System Danke Viele Gruesse
Hallo ich habe dasselbe Problem wie oben schon beschrieben. Ich versuche derzeit einen IIR-Filter höherer Ordnung auf FPGA-Ebene zu implementieren. Die Eckdaten des Filters wie Abtastfrequenz Grenzfrequenz, Sperr-und Durchlassdämpfung sind bekannt. Die dazugehörigen Filterkoeffizienten wurden mittels Matlab ermittelt. Die einzelnen Komponenten wie der Addierer, Multiplizierer und Verzögerer habe ich schon in VHDL entworfen. Allerdings ist es bekannt, das IIR-Filter dazu neigen instabil zu werden, sofern Rundungsfehler der Berechnungen oder der Koeffizienten vorhanden sind. Ich habe jetzt schon meine Koeffizienten auf VHDL-Eben 30 Bit groß gemacht und frage mich ob dies immer noch nicht genau genug ist. Ich versuche derzeit Lösungen zu finden wie ich meinen Filter stabilisieren kann. Würde mich auf ein paar Ideen freuen Gruß
> Ich habe einen FIR Filter mit der Übertrafungsfunktion H(z) und dies
soll ich invertieren als H'(z) = 1/ H(z) und ich soll jetzt H'(z)
stabilisieren.
H'(z) hat einige Pole, die außerhalb des Einheitskreises liegen.
Ist das nicht schon vom Ansatz her falsch?
Du machst ein Filter das bei bestimmten Frequenzen unendliche
Verstärkung hat. Wie soll das in der Anwendung jemals funktionieren ohne
zu schwingen?
>>> einen IIR-Filter höherer Ordnung auf FPGA-Ebene zu implementieren.
Höhere Ordnung ist schwierig. Du zerlegst das Filter in
hintereinandergeschaltete Teilfilter 2.Ordnung, das ist stabiler, besser
zu übersehen und unempfindlicher gegen nichtlineare Effekte, z.B.
Grenzzyklen. Zur Paarung der Pole mit den Nullstellen gibt es auch
einiges zu sagen, da muss Du mal kucken.
Cheers
Detlef
Hallo Detlef, Danke für deine Antwort. Das was du angesprochen hast habe ich berücksichtigt. Ich habe den Filter in Direktform2, SOS implementiert. Sprich in Kaskaden Form (Mehrere Sektionen 2.Ordnung). Ich weiß auch dass wenn alle Sektionen für sich stabil sind, dass das gesamte System stabil ist. Ich schaffe es aber nicht einen einfachen Bandpass 4.Ordnung mit folgenden Eckdaten zu stabilisieren (Implementierung auf FPGA-Ebene). Eckdaten: Fs=4882.8125Hz Fc1=500Hz Fc2=800Hz Entwurfs mittels FDATool Matlab mit Einstellung Elliptic,Equiripple. Zusätzlich weiß ich nicht wie sich die Auswirkung der Instabilität auf der Hardware bemerkbar macht. Ich sehe am Oszilloskop ein Ausgangssignal obwohl ich am Eingang des AD-Wandlers nichts Einspeise. Ich gehe davon aus das dies das instabile Verhalten des Filters ist. Vielleicht könntest du auch mal genauer die Sache mit der Pol-und Nullstellen Paarung genauer erläutern. Vielen Dank im Voraus Gruß
Quantisiere doch erst mal die Koeffizienten, und kontrolliere ob das Filter mit den quantisierten Koeffizienten stabil ist bzw. überhaupt noch den gewünschten Frequenzgang hat. Gerade elliptische Filter reagieren sehr empfindlich auf Quantisierung. Erst wenn das alles passt kannst du dir Gedanken um "Online"-Quantisierungseffekte machen (Quantisierungsrauschen, Grenzzyklen). Sicherlich ist ein Oszilloskop aber nicht das geeignete Debugging-Instrument für ein FPGA. Eine bitexakte Nachbildung des Filters in MATLAB oder ein VHDL-Simulator wären da sinnvoller.
Diesen Vorgang habe ich schon durchgeführt. Folgendes Habe ich gemacht: -Filterkoeffizienten per Matlab bestimmt. -Koeffizienten von Dezimal in Dual Festkommazahl (Vektor 30 Bit) umgerechnet -Dual Festkommazahl wieder in Dezimal Kommazahl umgerechnet (abgerundete Filterkoeffizienten) Die abgerundeten Filterkoeffizienten habe ich in Matlab eingegeben und geprüft ob der Filter stabil ist. Matlab liefert mir als Ergebnis das der Filter Stabil ist. Trotzdem Funktioniert das leider nicht.
>> Matlab liefert mir als Ergebnis das der Filter Stabil ist. Wie macht es das, indem es die Pole mit den gerundeten Koeffizienten berechnet? Du musst, wie schon gesagt, das Filter bitgenau simulieren. Das geht in Matlab, indem Du die integer-Rechnung des FPGAs nachvollziehst, mit allen Rundungen der Zwischenergebnisse. Ausgangssignal ohne Eingangssignal kann Grenzzyklus sein, wahrscheinlicher ist aber ein Implementierungsfehler ;-) . Mit der richtigen Paarung der Pole und Nullen kannst Du die benötigte Dynamik (Bitbreite) runterdrücken, auch gut zu simulieren. mail mal die Filter. Cheers Detlef
die abgerundeten Filterkoeffizienten berechne ich selber. Ich habe mir dazu ein C-Programm geschrieben das mir die originalen Filterkoeffizienten aus einer Datei einliest, die Berechnung der gerundeten Filterkoeffizienten durchführt und die gerundeten Filterkoeffizienten wieder in eine Datei ausliest. Somit kann ich die gerundeten Filterkoeffizienten schneller implementieren ohne das per Hand durchzuführen. Habe auch verifiziert ob die gerundeten Werte von dem Programm übereinstimmen. Um die Stabilität hier zu prüfen nutze ich wieder das FDATool. Der Unterschied ist nur das ich einen Filter importiere. Ich gebe dem FDATool die SOS und G Matrix vor mit den gerundeten Koeffizienten und lasse mir das Pol und Nullstellendiagramm anzeigen. Das mit der nach Implementierung in Matlab ist eine gute Idee. Ich werde mal versuchen den Filter in Simulink bis in Bitebene aufzubauen und zu Simulieren. So hattest du es bestimmt gemeint oder…? Was würdest du gerne von den Filtern sehen? Etwa den Entwurf von dem FDATool? Gruß
>>Was würdest du gerne von den Filtern sehen? Zwei SOS sind 2*5 Zahlen, zehn integer Werte mit ihrer Skalierung, die hast Du doch implementiert (oder kann FDATool VHDL exportieren?). SOS Repräsentierung des FDATools kenne ich nicht, G-Matrix (gain?) auch nicht, vermutlich liegt die Lösung in einer korrekten Bedienung des tools. Du hast C drauf, darin kann man integer Filter auch gut simulieren. Cheers Detlef
Hallo Detlef, ich erhalte für den oben genannten Filter von Matlab folgende Koeffizienten: G-Matrix (Gain Faktoren) G1=0,170467974908347 G2=0,170467974908347 G3=1 SOS-Matrix Filterkoeffizienten Sektion1 Zähler: Nenner: B1=1 A1=1 B2=0 A2=-1,38131865455268 B3=-1 A3=0,788859239856839 Filterkoeffizienten Sektion2 Zähler Nenner: B4=1 A4=1 B5=0 A5=-0,983403466910449 B6=-1 A6=0,734577208191707 Die Nenner Koeffizienten werden alle noch mit -1 multipliziert. Das bedeutet dass sich ihr Vorzeichen ändert. Das mit der Skalierung habe ich noch nicht ganz verstanden. Soweit ich verstanden habe darf die Gesamtverstärkung nicht größer als 1 werden. Falls ja wie kann das realisiert werden. Gruß
Angehängte Dateien:
-
Frequnzgang.jpg
25 KB -
freqdiff.jpg
23 KB -
impulsdiff.jpg
25 KB
Das Ding ist ein absolut gutmütiger Bandpass, Beträge der Pole sind
0.888 bzw. 0.857. Wenn die Beträge 0.99 wären könnte das schwierig
werden, aber die Pole sind Lichtjahre vom Einheitskreis weg.
Frequenzgang.jpg zeigt den Frequenzgang des Filters. Dann habe ich mal
grob auf 13Bit quantisiert, siehe angehängtes script, freqdiff ist die
Differenz des quantisierten zum nicht-quantisierten, -60dB Unterschied.
impulsdiff.jpg ist der Unterschied in den Impulsantworten, alles gut.
Das Ding kannst Du gut in integer rechnen, da oszilliert nix, bei Dir
steckt noch ein Fehler in der Implementierung.
>>Die Nenner Koeffizienten werden alle noch mit -1 multipliziert.
Nein, A(1) sowieso nicht und wenn die A(2) A(3) das Vorzeichen wechseln
wirds instabil.
Vllt. liegen Deine Probleme daran?
Cheers
Detlef
clear
%G-Matrix (Gain Faktoren)
G1=0.170467974908347;
G2=0.170467974908347;
G3=1;
fb1=[1 0 -1];
fa1=[1 -1.38131865455268 0.788859239856839];
fb2=[1 0 -1];
fa2=[1 -0.983403466910449 0.734577208191707];
H1=freqz(fb1,fa1);
H2=freqz(fb2,fa2);
H=G1*G2*H1.*H2;
plot(20*log10(abs(H)));
fff=8192;
fb1r=round(fb1*fff)/fff;
fa1r=round(fa1*fff)/fff;
fb2r=round(fb2*fff)/fff;
fa2r=round(fa2*fff)/fff;
H1r=freqz(fb1r,fa1r);
H2r=freqz(fb2r,fa2r);
Hr=G1*G2*H1r.*H2r;
figure(1);
plot(1:length(H),20*log10(abs(H)),'b.-',1:length(Hr),20*log10(abs(Hr)),'
r.-');
figure(2);
plot(20*log10(abs(H-Hr)))
n=200;
irep=filter(fb1,fa1,[1 zeros(1,n)]);
irep=filter(fb2,fa2,irep);
irepr=filter(fb1r,fa1r,[1 zeros(1,n)]);
irepr=filter(fb2r,fa2r,irepr);
figure(3);
plot(1:length(irep),irep,'b.-',1:length(irepr),irepr,'r.-')
plot(1:length(irep),irep-irepr)
Angehängte Dateien:
-
Filterstruktur_D2SOS.png
6,6 KB
Hallo Detlef, danke für deine Schnelle Antwort. Ich denke auch dass es ein Implementierungsfehler sein müsste. Ich habe mich falsch ausgedrückt bezüglich des Vorzeichenwechsels. A0 ist in dieser Form immer 1. Im Anhang findest du die dazugehörige Filterstruktur einer Sektion. Wie du sehen kannst, findet bei A1 und A2 noch ein Vorzeichenwechsel statt. Vielleicht liegt da ja mein Verständnisproblem. Bin gerade dabei mich im Bereich der Skalierung der Koeffizienten und Eingangsdaten einzulesen. Weiß aber noch nicht genau wie das funktionieren soll. Soweit ich jetzt weiß verhindert die Skalierung einen unerwünschten Overflow. Allerdings müssen hier gescheite Skalierungsfaktoren gewählt werden. Da ich noch ein Anfänger in diesem Bereich bin komm ich nur schleppend voran. Deshalb freue ich mich für jeden kleinen Denkansatz. Gruß
>>Vielleicht liegt da ja mein Verständnisproblem. Ja, für diese Struktur nimmst Du die Koeffs. wie sie sind, nix Vorzeichenwechsel. Ist der Klassiker für digitale Filter. >>Soweit ich jetzt weiß verhindert die Skalierung einen unerwünschten Overflow. Ja genau. Hängt davon was Dein FPGA nativ kann, vllt. 18x18Bit=36Bit? Das würde dicke reichen, das funzt schon für 13Bit Koeffizienten: Du skalierst 1=8192, multiplizierst und shiftest dann 13Bit runter. Wenn Deine Eingangsweite z.B. 16Bit breit sind benötigst Du 13+16=29 Bit Rechendynamik, passt in 36 Bit Rechnung rein (ähm, grob gerechnet, ein Bit kommt für Vorzeichen hinzu, ein Koeff. ist größer 1, nochmal ein Bit dazu, aber immer noch genug Dynamik übrig). Das Filter selber ist unkritisch, das muss gehen. Cheers Detlef
>> Hängt davon was Dein FPGA nativ kann, vllt. 18x18Bit=36Bit?
So wie ich weiß besitzt das Board zwölf 18-Bit Multiplizierer. Ich
verwende einen 12 Bit AD-und DA-Wandler. Somit ist mein
Eingangswertebereich 2047 bis -2048 (Zweier Komplement). Heißt das, dass
ich 1 Bit für das Vorzeichen, 11 Bit vor dem Komma und 13 Bit nach dem
Komma brauche (für alle Koeffizienten und Eingangsvektoren)? Könntest du
vielleicht die Quantisierung für einen Koeffizienten durchführen dann
würde ich das besser verstehen.
Gruß
>>. Heißt das, dass ich 1 Bit für das Vorzeichen, 11 Bit vor dem Komma
Ja, 12 Bit Wandler, 1 Bit Vorzeichen und 11 Bit Rest, das ist das x(n)
entsprechend Deinem Strukturbild.
Jetzt rechnest Du für ein SOS zwei Formeln entsprechend Deinem
Strukturbild:
(I) d(n)=x(n)-a1*d(n-1)-a2*d(n-2)
(II) y(n)=d(n)-d(n-2) (b0=1,b1=0,b2= -1)
Jetzt zur Dynamik. Bei Quantisierung der Koeffs auf 13 Bit (geht auch
mehr Bits, nur ein Beispiel):
A2=-1,38131865455268
A2*8192 ~ -11316=a1 (hat 15Bit)
A3=0,788859239856839
A3*8192 ~ 6462=a2 (hat 13 Bit)
G1=0.170467974908347 ~ 1396/8192
Also (in C notation)
(I) d(n)=x(n)-(11316*d(n-1)+6462*d(n-2)) >>13;
(II)y(n)=(1396*(d(n)-d(n-2)))>>13;
Das wäre mal nen erster hack. Das geht vllt. noch bißchen besser, weil
man in Formel (I) zuviele Bits für 18*18 Bit Multiplizierer wegwirft,
aber das sollte so erstmal funzen, allerdings ist es schon spät,
möglichweise ist nen bug drin.
Aber mal probieren, das geht ja in C einfach: integer Rechnung
implementieren, passband signal reinstecken, das muß rauskommen,
stopband Siganl muss unterdrückt werden.
Cheers
Detlef
Würde es aber nach dieser Berechnungsmethode nicht zu sehr starken Fehlern kommen? Oder habe ich das falsch verstanden? Ich führe mal ein einfaches Rechenbeispiel durch. Angenommen meine Eingangsweite wäre 4-Bit. Das würde heißen, dass der Wertebereich von 7 bis –8 beträgt. Ich skaliere auf 1=16 Angenommen ich habe einen Koeffizient der hat den Wert 1.75 und muss nun einen Multiplikation mit dem Eingangswert von 2 durchführen. Binärdarstellung: 1.75=1.1100 +2= 0010 Multiplikation: 0010x1.1100= 0011.1000 Rechtschifft um vier ergibt 0011 und Dezimal 3. Eigentlich ergibt der Wert 3.5. Da dieser Wert bei jeder Rückkopplung in die Berechnung eingeht, machen wir doch einen Fehler. Gruß
Yo, so ist das, die Rechnung stimmt. Du zeigst auch lediglich den Rundungsfehler: 2*1.75 ~ 3. Der Quantisierungsfehler für die Koeffizienten ist bei Deinem Beispiel nicht drin, weil 1.75*16 eine ganze Zahl ist, aber z.B. 0.788859239856839*8192 ist keine. Aber keine Sorge, das funktioniert schon, das ist keine Hexerei, die Fehler und flaws kriegst Du durch genügend Bitbreite kontrolliert. Cheers Detlef
Ich brauchte Spass und habe das erste Teilsystem G1=0.170467974908347; fb1=[1 0 -1]; fa1=[1 -1.38131865455268 0.788859239856839]; mal in integer gerechnet, C-Code angehängt. freq1 und freq2 zeigen Amplitudengang und Phasengang des Systems, mit Matlab berechnet. Die Eingangs/Ausgangssignale habe ich mit C berechnet und mit Matlab angezeigt. Bei der Resonanzfrequenz 112/512 ( 1 entspricht Nyquist Frequenz) erwarte ich +4dB ~ Faktor 1.6 und 0 Grad Phasenverschiebung, passt, siehe Bild freq3 Bei Nyquist/2 erwarte ich ca. -11dB ~ 0.28 und Phasenverschiebung von fast -pi, passt, siehe Bild freq4. Das funzt gut so, mach das so mal aufm FPGA! Cheers Detlef #include <conio.h> #include <stdint.h> #include <stdio.h> #include <math.h> /*******************************************/ int biquad(int i) /*******************************************/ { static int32_t q2=0,q1=0; int32_t r,rr; #define A1 ((int32_t)((-1.38131865455268*8192))) #define A2 ((int32_t)((0.788859239856839*8192))) #define GG ((int32_t)((0.170467974908347*8192))) r=i-((A1*q1+A2*q2)>>13); rr=q2; q2=q1; q1=r; return((GG*(q1-rr))>>13); } /****************************************************************/ int main() { int k=0,m=0,i; FILE *fo,*fod,*foc; if((foc=fopen("e:\\interim\\rein.txt","r"))==NULL){ printf("geht nicht auf \n");return(0);} if((fod=fopen("e:\\interim\\raus.txt","w"))==NULL){ printf("geht nicht auf \n");return(0);} do{ k=fscanf(foc,"%d \n",&i); if(k==EOF) break; m++; //sscanf(c,"%e\n",&t1); printf("%d %d \n",m,i); fprintf(fod,"%d\n",biquad(i)); } while(1); return(1); }
Thread beobachten
Seitenaufteilung abschaltenBitte 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.