IIR-Filter mit verschwindendem Fehler

[Walter T. (nicolas)]

Hallo zusammen,

ich will ein IIR-Filter erster Ordnung implementieren. Randbedingung 
ist, dass die Implementierung in Integer erfolgt, und dass das Ziel ist, 
dass sich der Rundungsfehler nicht akkumuliert, sondern sowohl der 
Fehler bei der Impulsantwort als auch bei der Sprungantwort dauerhaft 
verschwindet.

Oder anders gesagt: Die Differenz zwischen der Summe der Eingangswerte 
und Ausgangswerte wächst nicht an.

Jetzt habe ich allerdings gerade einen Knoten im Kopf, was die 
Formulierung des Fehlerterms angeht (siehe Anhang). Kann mir jemand 
einen Tipp geben?

[Analogmann (Gast)]

Was soll der Fehlerterm bringen? Willst du das Ergebnis am Ende 
irgendwann einmal korrigieren? Der IIR soll doch sicher ein permanentes 
Ergebnis ausgeben.

[Wolfgang (Gast)]

Analogmann schrieb:
> Was soll der Fehlerterm bringen? Willst du das Ergebnis am Ende
> irgendwann einmal korrigieren?

Solange der Erwartungswert für den Fehlertherm Null ist, braucht man 
nichts zu korrigieren und auch nicht bis zu irgendeinem Ende warten.

[Walter T. (nicolas)]

Wolfgang schrieb:
> Solange der Erwartungswert für den Fehlertherm Null ist, braucht man
> nichts zu korrigieren

Mal sehen: Ohne Fehlerterm ist bei einem Einheitssprung als Eingang 
u(t<0) = 0, u(t>=0) = 1  und alpha = 0.9 ist die Summe über alle u(k) 
unendlich. Die Summe über alle y(k) ist Null. Also wächst der Fehler y-u 
ohne Fehlerterm gegen minus unendlich. Um einen Fehlerterm komme ich 
also nicht herum. Ich weiß nur gerade nicht, wie ich das richtig 
herleite.

[c-heater (Gast)]

int filtIn = 0;
long filtLong = 0;
int  filtOut = 0;


//Filter erster Ordnung:
filtLong += (filtIn-filtOut);
filtOut = filtLong / 128; // Für Filterzeit: tau= 128 * Taskzeit

[c-heater (Gast)]

Nachtrag: filtLong ist deshalb "long", weil die Werte um die Faktor 
Filterzeit, hier also 128 größer sind als die Ein-und Ausgangswerte

[K. S. (the_yrr)]

Walter T. schrieb:
> Mal sehen: Ohne Fehlerterm ist bei einem Einheitssprung als Eingang
> u(t<0) = 0, u(t>=0) = 1  und alpha = 0.9 ist die Summe über alle u(k)
> unendlich. Die Summe über alle y(k) ist Null.

dann musst du (würde ich so sehen) die u(k) auch skalieren. alpha 
skalierst du ja mit 2^16. bei Int ist der kleinste Schritt halt 1, 
probier das ganze mal mit 2^8 oder so als Sprung. Wenn du bei Float den 
Sprung klein genug machst verschwindet der auch wieder durch 
Rundungsfehler. niemand würde erwarten dass ein Filter mit Float einen 
Sprung von 0.0000...x verarbeiten kann.

selbst wenn deine Fehlerkorrektur funktioniert, ist u(k) = y(k + zeit 
bis fehlerkorrektur greift). damit gewinnst du nichts, du müsstest u 
(und damit y) anders quantisieren bzw. auch hoch skalieren.

[Walter T. (nicolas)]

c-heater schrieb:
> [...]

Ja, so sah mein zweiter Ansatz auch aus. Der Nachteil an dieser 
Formulierung besteht darin, dass der Fehler für Summe(y(t)) immer noch 
maximal 128-1 (bei mir 20-1) ist. Irgendwie habe ich naiv erwartet, dass 
der Fehler auf Plusminus 1 zu bekommen ist.

P.S.: Hoppla....bei den "exponnentiellen" Filtern (IIR 1. Ordnung) sind 
mir die "n" ein wenig verruscht.

[c-heater (Gast)]

Hallo,
nein, da entsteht kein Fehler. Das Filter läuft nach einer gewissen Zeit 
auf genau den Eingangswert.
Hier als Bild und Tabelle das Beispiel mit Eingangswerten 0,100,50,71,19 
und Filterzeit 20. Wobei in C sich Filterzeiten 16,32,64 anbieten, weil 
keine Division erforderlich ist.

[c-heater (Gast)]

Bei

filtOut += (filtIn - filtOut) / 128;

hättest Du allerdings das von Dir beschriebene Problem. Nicht aber, wenn 
Du die 2-zeilige Version von oben nimmst. Der Preis dafür ist allerdings 
eine zusätzliche Variable, die möglicherweise einen breiteren Datentyp 
erfordert.

[Walter T. (nicolas)]

c-heater schrieb:
> Das Filter läuft nach einer gewissen Zeit
> auf genau den Eingangswert.

Das stimmt. Allerdings stimmen die Zwischenwerte nicht. Bei einem 
linearen Filter müsste bei den folgenden Varianten das gleiche Ergebnis 
herauskommen:

deltaU = input();
u += deltaU;
y0 += filterExp(deltaU, &Status0);
y1  = filterExp(u, &Status1);

Bei den in double implementierten Filtern ist das auch so. Aber beim 
Integer-Filter ist noch der Wurm drin.

[K. S. (the_yrr)]

Walter T. schrieb:
> Bei den in double implementierten Filtern ist das auch so. Aber beim
> Integer-Filter ist noch der Wurm drin.

wenn alpha >0.5 ist kommt der Integer Filter immer nur auf +-1 an den 
Sollwert, sowohl von oben als auch von unten.

es ist halt nicht einfach und/oder gut bei einem Filter, der mit 
Bruchteilen des Eingangswertes arbeitet, die kleinste Quantisierung (bei 
Integer ist das halt 1) als signifikante Stelle zu benutzen.

in float folgt auf die 1 die 1.0000001192092896 als nächste mögliche 
Zahl, berechne mal den Sprung in einem Filter nur in Float, da kommt 
auch nur Müll raus (da du hier die Zwischenberechnung mit 2^16 
hochskalierst, müsste man dann fairerweise die Berechnung in Double 
machen und dann nur u(k) und y(k) wieder nach float konvertieren, ändert 
aber nichts daran, dass der Filter bei alpha>0.5 nicht springen wird => 
Fehler wird unendlich)

nimm bei deinem ersten Ansatz statt u(k) mal u(k)*2^8, speicher y(k) so 
wie es aus der Formel kommt für den nächsten Schritt (für y(k-1)), und 
nimm dann nur für die Ausgabe y(k)/2^8. die 2^8 kannst du relativ 
beliebig variieren. dann funktioniert der Filter problemlos ohne 
jeglichen Fehler.

[J. S. (engineer)]

Die einfachste Möglichkeit, das Problem zu umgehen, ist nach wie vor, 
den Eingangswert des Filters zu verrauschen. Der maximale Abstand zum 
Real-Wert ist dann maximal temporär 1 digit. Mache ich seit X Jahren in 
all meinen Synthies und vielen Integer-basierten Rechnungen im FPGA so.

[Walter T. (nicolas)]

Jürgen S. schrieb:
> ist nach wie vor, den Eingangswert des Filters zu verrauschen.

Gehört habe ich davon schon, aber mich noch nie damit beschäftigt. 
Kennst Du dazu empfehlenswerte Literatur?

[c-heater (Gast)]

>Das stimmt. Allerdings stimmen die Zwischenwerte nicht.
Aber so: Jetzt wird auch der Rest berücksichtigt, der bei der 
Integer-Division ansonsten verloren geht. Ist allerdings die Frage, ob 
es nicht einfacher ist, einfach alle Eingangswerte um TAU größer zu 
skalieren.

#define TAU 128
int filtIn = 0;
long filtLong = 0;
int  filtOut = 0;
inr filtRest = 0;


//Filter erster Ordnung:
filtLong += (filtIn-filtOut);
filtOut = (filtLong+filtRest) / TAU ;
filtRest = (filtLong+filtRest ) - filtOut * TAU ;

[Walter T. (nicolas)]

c-heater schrieb:
> Aber so

Erst einmal: Danke. Das funktioniert tatsächlich so. Egal, ob das Filter 
inkrementell oder auf den Wert angewendet wird, stimmt jetzt das 
Ergebnis.

Jetzt muß ich mir noch einmal in Ruhe anschauen, worin genau der 
Unteschied zu meinem Ansatz am Montag besteht. Da habe ich die Summe aus 
filtRest und filtLong gespeichert.

Fun fact: Das Filter ist zwar ein Tiefpass, aber nur im Rahmen der 
Integer-Auflösung.