Forum: Digitale Signalverarbeitung / DSP Biquadfilter instabil


Announcement: there is an English version of this forum on EmbDev.net. Posts you create there will be displayed on Mikrocontroller.net and EmbDev.net.
von M. H. (bambel2)


Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Guten Tag,

ich spiele gerade etwas mit der CMSIS DSP Library des Cortex-M4 herum.
Ich habe ein Biquadfilter mit folgenden Koeffizienten:

Zählerpolynom:
b0 = 0.16667
b1 = 0.11111
b2 = 0.16667

Nennerpolynom:
a0 = 1
a1 = -0.88889
a2 = 0.33333

Laut Octave sind die Nullstellen dafür:
  -0.33333 + 0.94281i
  -0.33333 - 0.94281i

und die Polstellen:
   0.44444 - 0.36851i
   0.44444 + 0.36851i

Die Polstellen liegen also schön weit im Einheitskreis der z-Ebene.
Den Frequenzgang und den Pol-Nullstellenplot habe ich angehängt.

Ich habe nun die Filterkoeffizienten, wie oben angegeben in meinen 
Mikrocontroller gehackt und das Ganze laufen lassen. NAch einigen 
Tausend Samples, läuft allerdings das Filter komplett weg und schlägt an 
die +/- inf Grenze der 32bit floating point Zahlen. Der Eingangswert ist 
eine Floating Point Zahl die aktuell relativ konstant ist. Wert: 3050 
(+- Rauschen von 10).

Das Filter selbst ist mir nicht so wichtig. Das war nur zum Testen.
Allerdings verstehe ich nicht, warum das Filter sich sofort ins Nirvana 
schwingt. Die Polstellen sollten soweit vom Einheitskreis entfernt 
liegen, dass selbst die Floating-Point Quantisierung das hier nicht 
nennenswert beeinflusst, oder?

Es wäre nett, wenn mir jemand sagen kann, warum das Ganze nicht stabil 
ist. Ich habe mit der praktischen Auslegung von Filtern keine/kaum 
praktische Erfahrungen.

Habt ihr noch praktische Tipps, wie ich beispielsweise mit Octave ein 
stabiles IIR Filter erzeugen kann? Mein Ansatz ist da:
1) Mit entsprechenden Funktionen, bspw: butter(...) ein Filter zu 
erstellen
2) mit zf2sos das Ganze in eine Kaskade aus Biquad-Filtern umbauen.

Allerdings habe ich damit bisher wenig Erfolg :)

: Bearbeitet durch User
von Helmut S. (helmuts)


Bewertung
0 lesenswert
nicht lesenswert
Das Problem liegt bestimmt nicht an den Filterkoeffizienten.
Vielleicht hast du einen Fehler bei der Umsetzung der 
Filterkoeffizienten in die Software gemacht oder du schickst ab und zu 
falsche Zahlen in das Filter.

: Bearbeitet durch User
von M. H. (bambel2)


Bewertung
0 lesenswert
nicht lesenswert
Hallo,

vielen Dank für deine Antwort.
Ich nutze den "fertigen Code" von ARM, den sie in einer Library zur 
Verfügung stellen. 
https://www.keil.com/pack/doc/CMSIS/DSP/html/group__BiquadCascadeDF1.html#ga812eb9769adde4fed6d5ff4968a011d3

Ich habe mich noch nicht daran gemacht, das selbst zu coden, da ich 
dachte, die Library wird das schon optimierter machen, als ich. Werde es 
aber morgen gleich versuchen.

Ich habe hier noch folgendes gelesen:

Dergute W. schrieb:
> Allgemein wird's unangenehmer (=muss genauer gerechnet werden), wenn
> Pole oder Nullstellen des Filters in der Naehe des Einheitskreises sind.

Mein Stand war bisher, dass die Nullstellen für die Stabilität 
unkritisch sind. Meine liegen direkt auf dem Einheitskreis. Kann das ein 
Problem sein?

Vielen Dank

Edit:

ich habe zum Test das Filter durch ein Moving average filter mit
b0 = 0.1
b1 = b2 = 0

a0 = 1
a1 = 0.9
a2 = 0

ausgetauscht. Das läuft wie erwartet und mittelt munter den 
Eingangswert. Die Funktion scheint also zu funktionieren. Werde es 
trotzdem mal händisch nachcoden.

: Bearbeitet durch User
von Helmut S. (helmuts)


Bewertung
0 lesenswert
nicht lesenswert
In dem Link steht

Some design tools use the difference equation
    y[n] = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] - a1 * y[n-1] - a2 * 
y[n-2]

Probier es mal mit umgekehrtem Vorzeichen für a1 und a2.

von Dergute W. (derguteweka)


Bewertung
0 lesenswert
nicht lesenswert
Moin,

M. H. schrieb:
> Mein Stand war bisher, dass die Nullstellen für die Stabilität
> unkritisch sind. Meine liegen direkt auf dem Einheitskreis.

Wuerde mich auch wundern, wenn die Nullstellen fuer die Stabilitaet ein 
Problem waeren. Was problematisch werden koennte, waere die 
Sperrdaempfung in der Umgebung der Nullstellen, wenn die eben durch 
irgendwelche Rundungs/Rechenfehler nicht genau auf'm Einheitskreis 
liegen, sondern etwas daneben. Dann ist da halt die Daempfung nur noch 
"hoch", nicht mehr unendlich.

Gruss
WK

von Detlef _. (detlef_a)


Bewertung
0 lesenswert
nicht lesenswert
Helmut S. schrieb:
> In dem Link steht
>
> Some design tools use the difference equation
>     y[n] = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] - a1 * y[n-1] - a2 *
> y[n-2]
>
> Probier es mal mit umgekehrtem Vorzeichen für a1 und a2.

Glaube ich auch. Die lib wird von umgedrehten Vorzeichen der a's 
ausgehen. Komisch ist allerdings, dass das Ding nach einigen Tausend 
samples explodiert. Das sollte eigentlich schneller gehen.

Cheers
Detlef

von Helmut S. (helmuts)


Bewertung
0 lesenswert
nicht lesenswert
Detlef _. schrieb:
> Helmut S. schrieb:
>> In dem Link steht
>>
>> Some design tools use the difference equation
>>     y[n] = b0 * x[n] + b1 * x[n-1] + b2 * x[n-2] - a1 * y[n-1] - a2 *
>> y[n-2]
>>
>> Probier es mal mit umgekehrtem Vorzeichen für a1 und a2.
>
> Glaube ich auch. Die lib wird von umgedrehten Vorzeichen der a's
> ausgehen. Komisch ist allerdings, dass das Ding nach einigen Tausend
> samples explodiert. Das sollte eigentlich schneller gehen.
>
> Cheers
> Detlef

Mit fällt da noch etwas anderes ein.
Falls du C++ verwendest: Legst du eventuell das Biquad-Objekt vor jedem 
Aufruf neu an und dir geht am Ende der Speicher aus?

von M. H. (bambel2)


Bewertung
0 lesenswert
nicht lesenswert
Detlef _. schrieb:
> Glaube ich auch. Die lib wird von umgedrehten Vorzeichen der a's
> ausgehen.

Tut sie nicht. Der Hinweis gilt dafür, dass manche Design-Tools das 
invertiert ausgeben.

Mein Moving average Filter
M. H. schrieb:
> ich habe zum Test das Filter durch ein Moving average filter mit
> b0 = 0.1
> b1 = b2 = 0
>
> a0 = 1
> a1 = 0.9
> a2 = 0
>
> ausgetauscht.

hätte sonst auch nicht getan.

Helmut S. schrieb:
> Falls du C++ verwendest
Reines C. Mit aktivierter FPU

von Jan K. (jan_k)


Bewertung
0 lesenswert
nicht lesenswert
Habs gerade mal mit single float Koeffizienten in ML durchgerechnet, 
RMSE < 2 Promille bei nem Eingangssignal von 3050+randn(10000,1). Also 
wie erwartet liegt es nicht an der Quantisierung.

Zeig mal bitte den Code, bisher hat die CMSIS Implementierung keine 
Probleme gemacht. Welche Version ist das? Da es der M4 ist, gibt's da 
noch irgendwelche Probleme mit der FPU oder so? Ggf. mal hardfloat 
ausschalten.

von M. H. (bambel2)


Bewertung
0 lesenswert
nicht lesenswert
habe den Code gestern rausgemüllt. Müsste es nochmal nachcoden.
Auf jeden Fall gut, dass du es mal nachprobiert hast. Ich schreibe 
gerade ien kleines C-Programm auf dem PC, was das auch mal testen soll.

Mit der CMSIS Lib habe ich auch das Problem gehabt, dass ich eigentlich 
ein Festkomma Filter wollte, die Funktionen aber immer gehardfaultet 
sind. Kann schon sein, dass mein Kompilat da sehr verbuggt ist. habe die 
.a-Datei aus einem relativ alten projekt und schleife sie seitdem herum.

Ich mache mich daran, dass Problem minimal auf dem STM nachzustellen.

von M. H. (bambel2)


Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Also: Ich habe das Ganze nun auf dem PC implementiert.

Ergebnis: Es läuft, genau wie auf dem uC, nach einigen 100 Takten weg. 
Habe das Programm angehängt.

Wäre nett, wenn mir jemand sagen kann, wo ich den Fehler eingebaut habe.

Die ersten paar Ausgangswerte:
0,508.252289
1,395.463257
2,1173.332886
3,444.518494
4,1351.419434
5,302.570648
6,1536.957031
7,90.332680
8,1787.458496
9,-203.080811
10,2131.770508
11,-606.939880
12,2605.526367
13,-1162.675049
14,3257.430664
Nach 522 Samples oszilliert der Ausgang nur noch zwischen -inf und +inf.

Edit:
Habe das Ganze nochmal mit double als Zahlenformat implementiert. Da 
läuft das Filter auch weg.

: Bearbeitet durch User
von M. H. (bambel2)


Bewertung
1 lesenswert
nicht lesenswert
Also. Habe das Ganze nun gelöst.

Das Problem ist meine Dummheit. Ich habe das Ganze jetzt nochmals auf 
Papier ausgekaspert. Natürlich muss bei einer Funktion

die Differenzengleichung

lauten. Also war es doch das mit dem Vorzeichen der Filterkoeffizienten. 
Habe mich von der beschreibung in der CMSIS Library komplett verwirren 
lassen.
Damit ist die Sache dann wohl gelöst.

Vielen Dank an euch.

: Bearbeitet durch User
von Jan K. (jan_k)


Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Hab deinen Code in ML getestet. Deine Vorzeichen sind falsch, siehe 
https://en.wikipedia.org/wiki/Digital_biquad_filter.

Anbei ein Plot, oben single vs Referenz, unten double. Geht beides :)

Code (sorry, quick + dirty)
clear mex, mex CFilter.c -R2018a, 
b = [.16667 .11111 .16667]; a = [1 -.88889 .33333]; 
% f = dfilt.df1(b,a) % für fvtool, braucht signal processing tb
ff = filter(b,a,3050.453*ones(1000,1));
[cf1, cf2] = CFilter; 

figure, subplot(211), hold all, 
plot(cf1), plot(ff,'--'), 
ylim([2000 4000]),
legend('single')

subplot(212),hold all, 
plot(cf2), 
plot(ff,'--'),
legend('double'),
ylim([2000 4000])

edit: Kacke, zu lahm :D

: Bearbeitet durch User
von Detlef _. (detlef_a)


Bewertung
0 lesenswert
nicht lesenswert
>>>>>>>>>>>>
> Probier es mal mit umgekehrtem Vorzeichen für a1 und a2.

Glaube ich auch. Die lib wird von umgedrehten Vorzeichen der a's
ausgehen.
<<<<<<<<<<<<<<

Yo. :)))))))))

Cheers
Detlef

von Joggel E. (jetztnicht)


Bewertung
0 lesenswert
nicht lesenswert
Und was soll Floatingpoint dabei ? Welches Verhaeltnis von der groessten 
und kleinsten Zahl benoetigst du ?
Welches ist der dynamische Bereich ?

Der Witz daran : Float32 hat 6 signifikante dezimalstellen bedeutet :
20'000'000 - 1 = 20'000'000
waehrend  Int32 9 signifikannt stellen hat :
2'000'000'000 - 1 = 1'999'999'999

Einfach float zu nehmen ist unguenstig.

von Detlef _. (detlef_a)


Bewertung
0 lesenswert
nicht lesenswert
>>>> Einfach float zu nehmen ist unguenstig.

Kann man so nicht sagen. Das hängt wie immer davon ab. Wenn ich mich um 
den Wertebereich in dem mein Ergebnis liegt nicht kümmern will/kann 
nehme ich float. Es wird schon irgendwie in -/+ 10^38 liegen. Wenn es 
rauschmässig auf die Rundungsfehler ankommt muss ich integer nehmen und 
sicher sein, dass kein Overflow auftritt. Die Genauigkeit ist dann 
besser.

Früher war Signalverarbeitung mit floats nogo, inzwischen kann der STM32 
das schneller als integer.

Cheers
Detlef

: Bearbeitet durch User
von Thomas U. (thomasu)


Bewertung
0 lesenswert
nicht lesenswert
Detlef _. schrieb:
> Früher war Signalverarbeitung mit floats nogo, inzwischen kann der STM32
> das schneller als integer.

Warum sollte das so gewesen sein?
Wenn man die Randbedingungen zu Rundung und Auflösung beachtet (was man 
in beiden Fällen muss!) ist das nie ein Problem.

Man muss nur binomische Formeln kennen um die Rechnungen aufzu spalten 
und das machen LIBs.

von Christoph M. (mchris)


Bewertung
0 lesenswert
nicht lesenswert
>Detlef _. schrieb:
>> Früher war Signalverarbeitung mit floats nogo, inzwischen kann der STM32
>> das schneller als integer.

Das glaube ich nicht. Wie viele Taktzyklen braucht eine Integer-Additon 
im Vergleich zur Float-Addition inclusive laden der Register auf dem 
STM?

Antwort schreiben

Die Angabe einer E-Mail-Adresse ist freiwillig. Wenn Sie automatisch per E-Mail über Antworten auf Ihren Beitrag informiert werden möchten, melden Sie sich bitte an.

Wichtige Regeln - erst lesen, dann posten!

  • Groß- und Kleinschreibung verwenden
  • Längeren Sourcecode nicht im Text einfügen, sondern als Dateianhang

Formatierung (mehr Informationen...)

  • [c]C-Code[/c]
  • [avrasm]AVR-Assembler-Code[/avrasm]
  • [code]Code in anderen Sprachen, ASCII-Zeichnungen[/code]
  • [math]Formel in LaTeX-Syntax[/math]
  • [[Titel]] - Link zu Artikel
  • Verweis auf anderen Beitrag einfügen: Rechtsklick auf Beitragstitel,
    "Adresse kopieren", und in den Text einfügen




Bild automatisch verkleinern, falls nötig
Bitte das JPG-Format nur für Fotos und Scans verwenden!
Zeichnungen und Screenshots im PNG- oder
GIF-Format hochladen. Siehe Bildformate.

Mit dem Abschicken bestätigst du, die Nutzungsbedingungen anzuerkennen.