PT1 Filter mit sehr kleine Koeffizienten und Rundung

Die Standardlösung ist, mit höherer Wortbreite (Genauigkeit) zu rechnen, 
das hilft Dir aber nicht weiter. Wie wäre es, wenn Du Dein PT1 mit viel 
niedriger Samplerate laufen lässt, dann werden die Koeffizienten ja 
wieder größer.
z.B. Faktor 1024 langsamer: Immer 1024 Eingangswerte aufaddieren und das 
Ergebnis durch 1024 Teilen. Dann das Ergebnis ins langsame PT1 Filter 
einspeisen.

Du kannst das Filter auch in Integer rechnen. Bei einem 32-Bit-Integer
mit Vorzeichen ist die Auflösung immerhin um den Faktor 128 höher, was
etwa 2 zusätzlichen Dezimalstellen entspricht.

Ob 32 Bit reichen, hängt davon ab, welche Auflösung die Ausgangswerte
haben müssen und wie fein k eingestellt werden können muss.

asd schrieb:
> Gibt es eine Standardlösung für diese problem, in der Art dass man auf
> eine Literaturstelle verweisen kann als "state of the art"?
Rechne den Regler in long. Dann hast du immerhin 9 signifikante 
Stellen...

EDIT: Fenster zu lange offen --> Zweiter...  ;-)

PTx schrieb:
> Deshalb nehm ich dann ein PT1 mit 1h Zeitkonstante, das mir alle 1sec
> einen Eingangswert nimmt, und schalte dem ein PT1 mit 1s Zeitkonstante
> und 1ms Abtastzeit davor.
>
> Das schnelle PT1 Filter wirkt hier also eher wie ein Antialiasing
> Filter.

Das Antialiasing ist bei dieser Dimensionierung leider unvollständig.
Wenn das Eingangssignal bspw. eine Frequenzkomponente von 1Hz enthält,
wird diese durch das Filter mit τ=1s nur um 3db gedämpft. Durch die
anschließende Abtastung mit 1Hz kann deswegen ein unerwünschter
Gleichanteil entstehen, der durch das nachgeschaltete Filter trotz
seiner großen Zeitkonstante von 1h nicht weggefiltert wird.

Walter L. schrieb:
> warum nicht mal mit der exakten Lösung des PT1-Gliedes arbeiten?
>
> for(;;)
> {
> y = exp(-t/TAU)* y0 + (1-exp(-t/TAU)) * u;
> y0=y;
> }

Ersetze exp(-t/TAU) durch 1-k, y durch out, y0 durch alt, u durch in und
fasse auf der rechten Seite der Gleichung alle Summanden, die den Faktor
k enthalten zusammen, dann erhältst du

  out = alt + k*(in-alt)

also genau die Gleichung, die der TE verwendet hat und deren Berechnung
mit 32-Bit-Floats zu große Rundungsfehler aufweist.


Auch wenn ich persönlich für solche Anwendungen die Integer-Arithmetik
bevorzuge (konstante Auflösung über den gesamten Wertebereich, Fehler
besser abschätzbar, auf CPUs ohne FPU oft effizienter usw.), habe ich
mich mal in einem Berechnungsverfahren für 32-Bit-Float- Arithmetik
versucht, das denselben Filtertyp wie der des TE, aber mit geringeren
Abweichungen vom Ideal realisiert.

Das Filter des TE (zeitdiskreter Tiefpass 1. Ordnung) wird durch
folgende Formel beschrieben, die in jedem Zyklus neu berechnet wird:

  y' = y + k · (x - y)


Das Problem:

Macht das Eingangssignal bspw. einen Sprung und bleibt dann konstant auf
einem non 0 verschiedenen Wert x0 stehen, nähert sich das Ausgangssignal
asymptotisch diesem Wert. In jedem Zyklus wird der Ausgangswert y um
Δy=k·(x-y) erhöht. Je mehr sich y x0 annähert, desto kleiner wird Δy.
Irgendwann ist Δy relativ zu y so klein, dass die Addition auf Grund der
begrenzten Auflösung der Zahlendarstellung zu keiner Änderung von y mehr
führt. Zu diesem Zeitpunkt ist Δy zwar winzig, trotzdem hat y zum
Zielwert x0 noch einen Abstand von x0-y=Δy/k, was bei sehr kleinem k
(beim TE ist k=1e-6) ein nicht zu vernachlässigender Fehler sein kann.


Lösungsansatz mit Einschränkungen:

Das Problem stellt sich nicht, wenn x0=0 ist, sich y also asymptotisch
der 0 annähert. In diesem Fall reduziert sich obige Formel zu

  Δy = -k · y

Ist k größer ist als die relative Auflösung der Zahlendarstellung (bei
32-Bit-Floats 2**-23≈1.2e-7) führt dieses Δy bei der Addition solange
zu einer Veränderung von y, bis x=0 wird.

Bei einem bekannten Sprung des Eingangssignals könnte man sich diese
Eigenschaft zunutze machen, indem man statt des Ausgangssignals y dessen
Abstand u zu x0, also u=y-x0 berechnet. Substituiert man in der Formel
ganz oben y durch x0+u, erhält man

      x0 + u' = x0 + u + k · (x - (x0 + u))
  =>       u' =      u + k · ((x - x0) - u)

Die letzte Gleichung entspricht der ursprünglichen, nur dass y durch
u=y-x0 und x durch x-x0 ersetzt worden ist. Sowohl Eingangs- als auch
Ausgangssignal werden jetzt also nicht mehr absolut, sondern auf x0
bezogen betrachtet.

Diese Vorgehensweise hat eine gewisse Ähnlichkeit mit dem Kleinsignal-
verfahren in der Analogelektronik. x-x0 und u=y-x0 sind sozusagen die
Kleinsignale des Eingangs- und Ausgangssignals und x0 ist der Arbeits-
punkt, auf den sich beide Signale beziehen.

Aus u kann man mit y=x0+u jederzeit den Ausgabewert y berechnen, der
durch den kleinen Umweg bei der Berechnung jetzt tatsächlich gegen x0
konvergiert und nicht schon vorher auf einem anderen Wert stehenbleibt.

Einen wirklichen Nutzen bringt dieser Trick leider nur bei Eingangs-
signalen, die nach jeder Änderung über einen längeren Zeitraum konstant
sind, so dass man diesen konstanten Wert in der Berechnung jeweils als
Arbeitspunkt x0 verwenden kann. Zappelt das Eingangsignal x hingegen wie
wild hin und her, kann das Ausgangssignal y beliebig stark von x
abweichen, so dass u nicht klein genug ist, um die Auslöschungsfehler
wirksam zu bekämpfen.


Die eigentliche Lösung:

Was spricht aber dagegen, den Arbeitspunkt aus dem – auf Grund der
Tiefpassfilterung – deutlich weniger zappeligen Ausgangsignal zu
gewinnen?

Nichts.

Statt auf x0 aus dem Eingangssignal bezieht man die Signal also auf ein
y0, das man zu geeigneten Zeitpunkten aus dem Ausgangsignal pickt.
Anstelle von u=y-x0 nennen wir das Ausgangskleinsignal jetzt v=y-y0, und
die Berechnungsformel lautet jetzt:

  v' = v + k · ((x - y0) - u)

Da sich y nur langsam ändert, kann der Arbeitspunkt y0 über einen
längeren Zeitraum beibehalten werden. Erst wenn v einen gewissen Wert
übersteigt, so dass bei der Berechnung Auslöschungsfehler drohen, wird
y0 auf einen neuen Wert, nämlich das aktuelle y gesetzt. Damit wird v
wieder zu 0, und das Spiel beginnt von neuem.

Auf diese Weise wird der Fehler der berechneten Ausgangswerte
deutlich verringert, und das auch bei zappeligem Eingangsqsignal.

Bei einem über längere Zeit konstanten Eingangssignal kann das Ausgangs-
signal dennoch schon vor dem Erreichen der Asymptote stehenbleiben. Das
passiert immer dann, wenn v auf Grund von Rundungsfehlern gegen einen
von der eigentlichen Asymptote verschiedenen Wert konvergiert, und v
unterhalb der Grenze bleibt, an der y0 aktualisiert wird. Der dadurch
entstehende Fehler ist zwar deutlich kleiner als im ursprünglichen
Filteralgorithmus, trotzdem möchte man ihn natürlich gerne vermeiden. Am
einfachsten geschieht dies dadurch, dass y0 unabhängig vom Verlauf von v
immer nach Ablauf einer vorgegebenen Zeitspanne (bspw. 1τ) aktualisiert
wird.


Was bringt der ganze Aufwand nun an Genauigkeit?

Ich habe die obigen Überlegungen ebenso wie das Originalfilter des TE
und das zweistufige (kaskadierte) Filter von "PTx" in C-Funktionen
umgesetzt. Im C-Code wird v in ydelta und y0 in ybase umbenannt (viele
Programmierer bestehen ja auf etwas längere Variablennamen ;-)).

Die arithmetischen Ausdrücke in den drei Filterroutinen sind mit vielen
float-Casts gespickt, um die Rundungsfehler auf der Plattform des TE,
auf der es keine 64-Bit-Floats gibt, nachzubilden.

Ein weiteres Filter entspricht dem des TE, ist aber mit 64-Bit-Float-
Arithmetik implementiert, so dass die Fehler im Vergleich zu den anderen
Verfahren vernachlässigbar sind. Das Ausgangssignal dieses Filters wird
als Referenz bei der Bewertung der anderen Verfahren verwendet.

Als Eingangssignale für den Test wird ein Sprung-, ein Rechteck- und ein
Rauschsignal generiert (s. signalformen.png im Anhang). Die Berechnung
der Ausgangssignale der einzelnen Filter läuft über 20e6 Abtastzyklen.
Die jeweiligen Differenzen zum Referenzsignal sind in abweichungen.png
dargestellt.

Das Originalfilter weicht beim Sprungsignal um bis zu 0,03 von der
Referenz ab, was der Grund dafür war, dass der TE diesem Thread
gestartet hat. Auch beim Rechtecksignal ist der Fehler noch sehr groß.
Nur beim Rauschsignal werden zufriedenstellende Ergebnisse geliefert.

Das kaskadierten Filter liefert beim Sprung- und Rauschsignal nach dem
Einschwingen gute Ergebnisse, nur direkt nach dem Start bzw. nach dem
Sprung sind die Fehler größer. Beim Rechtecksignal zeigt das Filter aber
das oben beschrieben Abtastproblem, was in einem so großen Fehler
resultiert, dass die Kurve nicht zusammen mit den beiden anderen in
einem Digramm dargestellt werden konnte.

Bei dem Filter mit der hier beschriebenen dynamischen Verschiebung des
Arbeitspunkts sind die Fehler bei den Tests so klein, man sogar die
Granularität der 32-Bit-Float-Darstellung erkennen kann. Das bedeutet,
dass man (zumindest bzgl. der getesteten Signalformen) nicht mehr viel
optimieren kann.

Ich kann aber nicht garantieren, dass die Fehler dieses Verfahrens bei
beliebigen anderen Signalformen ebenso klein sind. Genauso wie das fiese
Rechtecksignal die ansonsten guten Ergebnisse des kaskadierten Filters
zunichte macht, gibt es vielleicht auch für mein Verfahren so ein
Killersignal.

Es war für mich zwar spannend, so etwas auszudenken, zu testen und zu
vergleichen, im realen Leben fühle ich mich aber mit einem einfacheren
Verfahren, das die gewünschte Genauigkeit mit geeigneteren Datentypen
(32- oder 64-Bit-Integer oder 64-Bit-Float) erzielt, deutlich wohler.
Wer aber – aus welchen Gründen auch immer – auf 32-Bit-Float-Arithmetik
besteht, für den kann dieses Verfahren eine gute Alternative sein.

asd schrieb:
> Jetzt komme ich in Bereiche in denen der Rundungsfehler dazu führt dass
> ich einen gewissen Totbereich habe

Mache die Rechnung in int, dann hast Du keinen Totbereich:
out = sum / k
sum = sum - out + in

Bei einer konstanten Folge "in" wird nach endlicher Zeit immer "out" == 
"in".
Man nimmt für k auch gerne 2-er Potenzen, dann geht schieben statt 
dividieren.