Festkomma IIR Filter

[Felix G]

Hallo Leute,

ich schreibe mir gerade ein IIR Filter, und das funktioniert auch soweit ganz gut, allerdings zur Zeit nur mit lahmen floats.


float filter(uint8_t order, float* coef_a, float* coef_b, float* buf, float x)
{
uint8_t n;

float y = coef_b[0] * x + buf[0];

for(n = 0; n < (order-1); n++)
buf[n] = buf[n+1] + coef_b[n+1] * x - coef_a[n] * y;

buf[order-1] = coef_b[order] * x - coef_a[order-1] * y;

return y;
}
Wie man sieht filtert diese Funktion die Werte quasi "on-the-fly", möchte man ein komplettes Array filtern muss sie eben in einer Schleife aufgerufen werden.


Im nächsten Schritt möchte ich das Ding jetzt auf 16-Bit Festkomma umbauen, es soll also im Endeffekt jedes "float" durch ein "int16_t" ersetzt werden.
Natürlich ist es nicht ganz so einfach, ein paar Steine versperren noch den Weg.


Mein Problem ist folgendes:
ich möchte mit dem int16_t einen Wertebereich von -1 bis 1 abdecken, es gibt (bis auf das Vorzeichen) also nur Nachkomma-Bits.
Das bedeutet natürlich auch, daß die Koeffizienten meines Filters in diesem Bereich liegen müssen... tun sie aber nicht.

Erst dachte ich es würde ausreichen einfach sämtliche Koeffizienten mit dem gleichen Faktor zu skalieren, aber Matlab hat mich da leider eines besseren belehrt.
Mit den skalierten Koeffizienten hatte mein Filter eine völlig andere Impulsantwort.


Also was tun?
wie kriege ich die Koeffizienten in einen Bereich zwischen -1 und 1 gequetscht?

[Felix G]

Also inzwischen habe ich den Fehler selbst gefunden...

prinzipiell kann ich sämtliche Koeffizienten durch einen Faktor f teilen, aber entsprechend muss ich auch y mit f multiplizieren, damit das Ergebnis stimmt.

die Zeile

y = coef_b[0] * x + buf[0];

muss also so aussehen

y = f * (coef_b[0] * x + buf[0]);

Bei meiner Festkomma-Implementierung werde ich für die Skalierung einen Faktor von 2 fest vorgeben. Das sollte für die meisten Filter ausreichen, und bringt noch etwas mehr Performance als eine Implementierung mit einstellbarem Faktor.

[Felix G]

Was ist los Leute, hat hier etwa noch Niemand so ein Filter implementiert!?

Naja, also falls mal Jemand sowas brauchen sollte, habe ich hier ein IIR Filter 2. Ordnung:
(wer nicht ganze Arrays, sondern einzelne Werte filtern möchte entfernt einfach die for-Schleife)

void filter_fr16_biquad_vect(int16_t* coef_a, int16_t* coef_b, int16_t* buf, int16_t* x, int16_t* y, uint16_t len)
{
uint16_t n;

for(n=0; n<len; n++)
{
// y(n) = x(n) * b0 + u0
y[n] = ((int32_t)x[n] * (int32_t)coef_b[0]) >> 15; // 1.15 * 1.15 = 2.30; 2.30 >> 15 = 1.15
y[n] += buf[0];
y[n] <<= 1;

// u0 = x(n) * b1 - y(n) * a1 + u1
buf[0] = (((int32_t)x[n] * (int32_t)coef_b[1]) - ((int32_t)y[n] * (int32_t)coef_a[0])) >> 15;
buf[0] += buf[1];

// u1 = x(n) * b2 - y(n) * a2
buf[1] = (((int32_t)x[n] * (int32_t)coef_b[2]) - ((int32_t)y[n] * (int32_t)coef_a[1])) >> 15;
}
}
Dabei müssen alle Werte zwischen -1 und 1 liegen, d.h. sowohl Eingangssignal als auch die Koeffizienten müssen entsprechend skaliert werden, wobei der Faktor für die Koeffizienten bei dieser Implementierung fest mit 0,5 vorgegeben ist.


So weit so gut, aber ein Problem habe ich doch noch...
betrachten wir mal ein in Matlab erzeugtes butterworth-Filter 2. Ordnung, dann fällt auf, daß die b-Koeffizienten bei niedrigen Grenzfrequenzen sehr klein werden können. Das ist natürlich schlecht, denn es verringert die Präzision erheblich, da mir massenhaft Nachkommastellen verloren gehen.

Erzeuge ich beispielsweise ein Filter mit dem Befehl butter(2, 0.008), dann sind die b-Koeffizienten (umgerechnt in das 16-Bit Festkomma-Format) gerademal 3, 5 und 3. Rechnet man diese wieder in floats um liegt der Fehler für b0 und b2 bei satten 15%, entsprechend beschissen sieht das gefilterte Signal dann auch aus.


Jemand eine Idee, was ich da tun kann?

[cmock]

du mußt ja nicht für alle zahlen das selbe festkommaformat verwenden, wenn du deine b-koeffizienten zb 8 bit weiter nach links schiebst, mußt du halt nach der multiplikation das ergebnis um 8 bit weiter nach rechts schieben...

das jetzt nur einmal als nicht-zuendegedachter ansatz.

cm.

[Felix G]

Ich fürchte das hilft nicht...


in der ersten Zeile berechne ich x * b0 und shifte das Ergebnis um 15 nach rechts, da ja sowohl x als auch b0 1.15 Festkommazahlen sind. Wäre b mit 2^8 hochskaliert müsste ich, um nach der Berechnung das richtige Ergebnis zu bekommen, nochmal um 8 nach rechts shiften. Ich könnte wetten, daß das Endergebnis in beiden Fällen völlig identisch ist.

[cmock]

zum glück hab ich das mit dem nicht-zuendegedacht dazugeschrieben, denn es ist in der tat sinnlos :-)

aber es will mir auch nix besseres einfallen; wenn, jetzt mal in realen zahlen gerechnet, die größenordnung deiner input- und output-signale um das 20000fache größer ist als die b-koeffizienten, dann stößt du an die grenzen des mit 16 bit darstellbaren...

cm.

[Felix G]

Noch schlimmer ist eigentlich, daß sich die a und b Koeffizienten so stark unterscheiden, sonst könnte ich ja einfach alle hochskalieren...

b0 = 0.00015514842347
b1 = 0.00031029684695
b2 = 0.00015514842347

a1 = -1.96446058020523
a2 = 0.96508117389914

Die a Koeffizienten sind also gut 6000x größer als die b Koeffizienten, aber so wie ich das sehe kann ich die Koeffizienten nicht unabhängig voneinander skalieren, ohne das Ergebnis zu beeinflussen.

Oder gibt es da vielleicht doch eine Lösung ?
Also die Gleichung für mein Filter sähe ja in diesem Fall etwa so aus:

y(n) = i * (b0*x(n)+b1*x(n-1)+b2*x(n-2)) - j * (a1*y(n-1) + a2*y(n-2))

Wenn i=j ist habe ich keine Probleme, dann muss ich den ganzen Kram nur durch i teilen (genau das mache ich ja schon, ich teile alle Koeffizienten durch 2, weil a1 sich sonst nicht darstellen lässt). Aber was soll ich tun wenn i und j nicht gleich sind?

[hacker]

Hast du dir dazu schon die App Note von Atmel angeschaut? Diese Probleme werden darin auch angesprochen.

http://www.atmel.com/dyn/resources/prod_documents/doc2527.pdf

[Felix G]

Danke für den Tipp,

leider steht da auch nicht viel mehr als ich eh schon wusste, was dafür spricht, daß mit 16 Bit tatsächlich nicht mehr möglich ist (bei den beiden Beispiel-Filtern sind die Koeffizienten ja recht unkritisch). Das bedeutet ich müsste mehr Bits für meine Koeffizienten einplanen...

also mal am Beispiel von b0 (0,000077574211735)

b0 1.15 = 0,000091552734375 (18% Fehler)
b0 1.19 = 0,000078201293945 (0,8% Fehler)
b0 1.23 = 0,000077605247498 (0,04% Fehler)

Naja, wenigstens muss ich nicht bis 32 Bit gehen sondern kann mich auf 24 beschränken (oder 20, aber so oder so brauch ich ein Byte mehr). Wenn ich für das Eingangssignal bei 16 Bit bleibe wären meine Zwischenergebnisse 40Bit breit, wobei ich die unteren 24 Bit wegwerfen kann.


Ok, so werd ichs wohl implementieren...
leider macht das mein Filter langsamer, da ich mehr Multiplikationen und Additionen brauche und ich vermutlich zumindest die Koeffizienten im SRAM halten muss (bei 16 Bit hätte ich alles in Register quetschen können, das wär dann richtig schön schnell gewesen).