Gleitender Mittelwert -> Algorithmus

[Maxe (maxemaxe)]

Hallo zusammen, ich brauch für die SPS einen 
Moving-Average-Filter/Gleitenden Mittelwert. Das ist ja eigentlich auch 
leicht und effizient zu programmieren, indem man dem Durchschnitt 
jeweils den neu hinzukommenden Wert anteilig addiert und den 
herausfallenden Wert anteilig subtrahiert. Jetzt frag ich mich, ob man 
dabei einen ungewollten Drift reinbekommen kann, wenn man mit 
Gleitkommazahlen arbeitet, weil es beim addieren und subtrahieren des 
gleichen Wertes evtl. zu Rundungen bzw. ungleichen Rundungsfehlern 
kommt.

Kann jemand was dazu sagen? Ist das ein reales Problem oder gleicht sich 
das statistisch sowieso aus? Oder muss man die Zahlen auf bestimmte 
Weise wählen oder Vorrunden um nicht in den Bereich der letzten 
signifikanten Bits zu kommen?

Danke!

  Average:= Average + (aInput - ValueTable[CurrentIndex])/length(ValueTable);
  ValueTable[CurrentIndex]:= aInput;

  CurrentIndex:= CurrentIndex + 1;
  if CurrentIndex > high(ValueTable) then CurrentIndex:= 0;

[N. M. (mani)]

Maxe schrieb:
> Jetzt frag ich mich, ob man dabei einen ungewollten Drift reinbekommen
> kann

Klar, kann immer passieren. Gerade bei Integration.
Aber wenn du beim reinstecken in den Filter den gleichen Fehler machst 
wie beim rausholen...

[Maxe (maxemaxe)]

N. M. schrieb:
> Aber wenn du beim reinstecken in den Filter den gleichen Fehler machst
> wie beim rausholen...

Das ist ja gerade die Frage, ob der Fehler gleich bleibt.

Ein Extrembeispiel wäre, wenn der Durchschnittswert beim "Reinstecken" 
des Wertes so groß ist, dass die kleine zugefügte Zahl komplett 
verschwindet. Also dass das kleinste Bit der Mantisse wertmäßig noch 
größer ist als der ganze zugefügte Wert. Zum Zeitpunkt, wo der Wert 
wieder aus dem Mittelwert fällt, könnte der Durchschnitt klein genug 
sein, dass sich der Wert doch wieder abbildet. Dann hätte man zwischen 
Reinstecken und Rausholen einen negativen Drift.

[Olly T. (twinpeaks)]

Das habe ich mich auch schon gefragt kann Dir aber dazu auch keine 
definitive Antwort geben. Ich vermute aber auch, dass das unter 
ungünstigen Umständen wegdriften könnte.

Aber falls die Eingabedaten eigentlich Ganzzahlen sind, z.B. wenn sie 
von einen ADC kommen, oder sich darauf abbilden lassen, mache ich das 
immer so, dass ich im Ringpuffer die Integerwerte speichere und damit 
nicht den Mittelwert sondern die Summe über die Fenstergröße berechne, 
und erst hinterher durch die Länge teile. Da dann nur Integerwerte 
addiert und subtrahiert werden, gibt es keine Rundungsfehler. Der 
Wertebereich von "Sum" muss natürlich groß genug sein.

Analog zu Deinem Beispiel:

  Sum := Sum + aInput - ValueTable[CurrentIndex];  // hier alles Integer
  ValueTable[CurrentIndex]:= aInput;
  Average := Sum / length(ValueTable);             // erst hier Float
  ...

[Marci W. (marci_w)]

Hallo Maxe,

also so ein Filter habe ich auch schon in einer SPS (Siemens S7/300) 
realisiert. Für meine Zwecke habe ich dabei einfach jeweils von n 
vorausgegangenen Messwerten den Mittelwert berechnet. Mit jedem neu 
hinzugekommenen Messwert Bin ich dann einen Messwert weiter gerückt und 
habe wieder von den letzten n Werten den Mittelwert berechnet. Problem 
dabei ist lediglich, dass man erst ab mindestens n verfügbaren 
Messwerten einen "echten" Wert bekommt.

Hat zumindest für meine Anwendung (Messkurven glätten) hervorragend 
funktioniert.

Ob das jetzt genau der Definition von "gleitendem Mittelwertfilter" 
entspricht, weiß ich nicht. Wie gesagt, es hat einwandfrei funktioniert, 
und irgendwelche Probleme mit Rechenungenauigkeiten, Auslöschungen etc. 
gab es dabei nicht.

ciao

Marci

[Marci W. (marci_w)]

Olly T. schrieb:
> dass ich im Ringpuffer die Integerwerte speichere und damit
> nicht den Mittelwert sondern die Summe über die Fenstergröße berechne,

Ja klar, wie denn sonst? Ich wäre gar nicht auf eine andere Idee, wie 
z.B. die des TO, gekommen.

ciao

Marci

[Maxe (maxemaxe)]

Olly T. schrieb:
> Aber falls die Eingabedaten eigentlich Ganzzahlen sind, z.B. wenn sie
> von einen ADC kommen, oder sich darauf abbilden lassen, mache ich das
> immer so, dass ich im Ringpuffer die Integerwerte speichere und damit
> nicht den Mittelwert sondern die Summe über die Fenstergröße berechne,
> und erst hinterher durch die Länge teile. Da dann nur Integerwerte
> addiert und subtrahiert werden, gibt es keine Rundungsfehler. Der
> Wertebereich von "Sum" muss natürlich groß genug sein.

Das wäre natürlich eine Möglichkeit. Ist auch gut zu wissen, dass nicht 
nur ich die Befürchtung habe und es die Problematik vielleicht 
tatsächlich gibt.
In meinem Fall gehen berechnete Werte in den Filter, da würde ich gerne 
bei Gleitkommazahlen (als Eingang) bleiben. Meine letzte Überlegung war, 
dem Filter einen einstellbaren Eingangsbereich (Definitionsbereich) zu 
verpassen und dann diesen Bereich geteilt durch die Anzahl Stützpunkte 
auf eine 32Bit-Ganzzahl zu projizieren, also skalieren. Die Summe dann 
als 64Bit führen und beim Auslesen auf die Gleitkommazahl 
zurückskalieren. Insgesamt halt recht aufwendig.

Marci W. schrieb:
> Für meine Zwecke habe ich dabei einfach jeweils von n
> vorausgegangenen Messwerten den Mittelwert berechnet. Mit jedem neu
> hinzugekommenen Messwert Bin ich dann einen Messwert weiter gerückt und
> habe wieder von den letzten n Werten den Mittelwert berechnet.
> [...]
> Ob das jetzt genau der Definition von "gleitendem Mittelwertfilter"
> entspricht, weiß ich nicht.

Doch schon, genau so ist der gleitende Mittelwert ja definiert. Der 
Algorithmus mit dem Hinzufügen und Herausnehmen ist eher eine 
Optimierung in Bezug auf die Rechenzeit. Wenn ich einen gleitenden 
Mittelwert mit 4000 Stützpunkten habe, müsste ich die ja alle in einem 
SPS-Zyklus aufaddieren. Mit der relativen Berechnung sind es dagegen nur 
3 Rechenoperationen. Mag sein, dass auch die Tausende Rechnoperationen 
kein Problem für die SPS sind, aber ich will mir da auch keinen 
künstlichen Flaschenhals schaffen.

> Problem
> dabei ist lediglich, dass man erst ab mindestens n verfügbaren
> Messwerten einen "echten" Wert bekommt.
Das ist ja prinzipbedingt so. Man geht da einfach von einem Anfangswert 
aus, initialisiert also das ganze Array mit einem sinnvollen 
Anfangswert.

[J. S. (engineer)]

Das so anzusetzen, ist gelinde gesagt "Murks"!

Man schleppt nämlich zusätzlich zum eigentlichen Wert noch das Integral 
des Rundungsfehlers mit und das ist nicht automatisch statistisch Null, 
sondern es ergibt sich ein beliebig tief gehendes Spektrum.

Richtig ist, vollständig zu addieren und jeweils neu zu dividieren. Dann 
gibt es nur einen temporären Rauschfehler, der erheblich kleiner ist und 
spektral auf den Bereich der betrachteten Punkte begrenzt ist.

[Marci W. (marci_w)]

Maxe schrieb:
> Mittelwert mit 4000 Stützpunkten

Ups, das ist natürlich ne Menge. Es würde mich interessieren, wo du 
derart viele Werte berücksichtigen musst.

Beim Vorrücken des Index den "rausfallenden" Wert am Anfang des Fensters 
von der Summe subtrahieren und den neu hinzugekommenen Wert addieren 
(ähnlich wie Deine Idee). Dann wieder teilen. Dann wäre die Rechnerei 
effizienter. Allerdings nur bei Ganzzahlen. Sonst hast Du wieder 
Rechenfehler, die sich aufsummieren können (da bin ich jetzt aber zu 
faul zu überlegen, ob das auch zutrifft).

Aber wie gesagt, mir kommt eine derart große Anzahl von Werten komisch 
vor.

SPS-mäßig könntest Du die Funktion in einen eigenen Baustein packen und 
diesem dann, wenn neue Messwerte ankommen, ein Startsignal beim Aufruf 
mitgeben. Im Baustein dann in jedem Zyklus nur eine kleinere Wertemenge 
addieren.
Kommt halt alles sehr auf Deine Randbedingungen an (wie schnell kommen 
die Messwerte rein, wie schnell brauchst du den gleitenden Mittelwert, 
wie schnell ist Deine SPS, was hat Sie sonst noch zu tun, max. zulässige 
Zykluszeit etc. pp.), das übliche halt.

ciao

Marci

[Maxe (maxemaxe)]

J. S. schrieb:
> Man schleppt nämlich zusätzlich zum eigentlichen Wert noch das Integral
> des Rundungsfehlers mit und das ist nicht automatisch statistisch Null,
> sondern es ergibt sich ein beliebig tief gehendes Spektrum.
Ok, dann sind wir schon zu dritt, die ein Problem in den 
Gleitkommazahlen sehen. Arbeitet man mit Ganzzahlen integrieren sich die 
Rundungsfehler ja nicht auf, sondern "fallen" wieder raus, wenn der Wert 
abgezogen wird. Man hat dann als Rundungsfehler nur noch die Summe der 
Rundungen im Ringspeicher. Ich würde vermuten, dass die sehr klein sind, 
man täuscht sich aber leicht.

Also, bei 4000 Stützstellen kann im schlechtesten Fall jeweils das 
letzte Bit falsch sein, also summiert ein Wert von 4000 oder 12 Bit. Von 
32 Bit bleiben dann noch 20 Bit verteilt auf den Eingangsbereich. Dürfte 
gut genug sein.

[Maxe (maxemaxe)]

Marci W. schrieb:
> Es würde mich interessieren, wo du
> derart viele Werte berücksichtigen musst.

Eine Regelung in der Verfahrenstechnik, das Fluid ist 3 Minuten in der 
Leitung unterwegs, also 180 Sek. Bei 10 Zyklen pro Sekunde hat man 1800 
Werte. Die Zyklenzeit könnte man wohl verlängern, und es gäbe wohl auch 
kaskadische Algorithmen für den gleitenden Mittelwert. Muss man aber 
auch alles wieder durchdenken.

[Michael D. (nospam2000)]

Wenn du das mit Integerwerten machst und die Summe der Werte speicherst 
anstatt dem Mittelwert und den Mittelwert erst für die Ausgabe 
verwendest, dann solltest du keine Rundungsfehler haben. Die 
summenvariable muss natürlich groß genug sein um den Wert auch aufnehmen 
zu können.

In etwa so

int32_t sum = 0; // muss so groß sein, dass 'count * MAX(next) reinpasst'
const int32_t count = 4;

int32_t calcMediumValue(int32_t next) {
 sum = sum - (sum / count) + next;
 return sum / count;
}

Wahrscheinlich ist es dasselbe was J. S. vorgeschlagen hat.

Übersehe ich was?

  Michael

[Max H. (nilsp)]

Für das Aufsummieren von Floating Point Werten gibt es den Kahan 
Algorithmus. Der ist wesentlich genauer weil er sich neben der Summe 
auch den laufenden Rundungsfehler merkt.

 https://en.wikipedia.org/wiki/Kahan_summation_algorithm

Lässt sich 1:1 für den gleitenden Mittelwert einsetzen. Die Subtraktion 
ersetzt Du einfach mit einer Negation und Addition.

Drift wird es immer noch geben, aber wesentlich weniger.

Was Du machen könntest wäre ständig eine zweite Instanz des Filters 
aufzubauen. Sobald die 4000 Samples erreicht sind, verwifst Du den 
aktiven Filter und schaltest auf den frischen um.

So kann nie mehr als die Drift von 4000 Samples ins System kommen weil 
Du immer neu ansetzt. Verdoppelt natürlich die Rechenzeit.

In Verbindung mit Kahan wird dieses Verfahren vermutlich sogar genauer 
sein, als in jedem Zyklus ganz naiv alle Werte zu addieren.

Wenn Du es auf die Spitze treiben willst, dann trenne zusätzlich noch 
die Additionen und Subtraktionen und berechne die Differenz erst dann, 
wenn Du sie brauchst. So vermeidest Du auch das Problem des 
Genauigkeitverlustes durch Auslöschung in der Subtraktion:

  https://de.wikipedia.org/wiki/Ausl%C3%B6schung_(numerische_Mathematik)

[Max H. (nilsp)]

Ich hab das mal eben in C zusammengehackt.

Ja, man kann vieles effizienter machen (Ringbuffer hust), aber es tut 
und sollte mit minimalem Drift arbeiten.

[Klaus K. (Gast)]

> Wenn Du es auf die Spitze treiben willst, dann trenne zusätzlich noch
> die Additionen und Subtraktionen und berechne die Differenz erst dann,
> wenn Du sie brauchst. So vermeidest Du auch das Problem des
> Genauigkeitverlustes durch Auslöschung in der Subtraktion:
>
>   https://de.wikipedia.org/wiki/Ausl%C3%B6schung_(numerische_Mathematik)

Aber nicht vergessen, vorher die sogenannte "Optimierung" des Compilers 
aus zu schalten, weil diese u.U. jegliche Genauigkeitsverbesserung 
zunuchte macht. Steht im Gleichen Artikel weiter unten:

https://en.wikipedia.org/wiki/Kahan_summation_algorithm#Possible_invalidation_by_compiler_optimization

"In principle, a sufficiently aggressive optimizing compiler could 
destroy the effectiveness of Kahan summation ...  if the compiler 
simplified expressions according to the associativity rules of real 
arithmetic"


Der bekannte intel c++ compiler ist wohl so ein Stück schlitzohriger 
Software, der es in diesem Fall mit der Optimierung (Einsparung von 
Rechenzeit) übertreibt.

Aber ob der Compiler "aggresiv optimiert" sollte der obligatorische 
Regressiontest (Vergleich Ergebnisse alter, neuer Softwarestand) zeigen.

[Max H. (nilsp)]

Da hast Du Recht, Klaus. Guter Einwand.

Ich benutz meist GCC, da ist dieser Spaß zum Glück per Default 
deaktiviert.

Wer doppelt sicher gehen will kann ja bei GCC mit -fno-associative-math 
das noch mal explizit deaktivieren.

[Justin S. (Gast)]

Was spricht dagegen, den hinzukommenden Wert (sagen wir 16-Bit) zum 
32-Bit-Akkumulator zu addieren und in den Buffer zu kopieren, und den 
herausfallenden Wert vom Akkumulator zu subtrahieren. Und dann jeweils 
den Akkumulator durch die Länge des Filters zu teilen (sinnvollerweise 
als Bitshift, bei den Registergrößen machbar bis Filterlänge 2^15)?

Operationen: eine Addition, eine Subtraktion, beide ohne Drift, ein 
Speichern, ein Lesen, ein Bitshift (mit vorheriger Addition der halben 
Filterlänge zur korrekten Rundung)?

[Michael D. (nospam2000)]

Justin S. schrieb:
> Was spricht dagegen, den hinzukommenden Wert (sagen wir 16-Bit) zum
> 32-Bit-Akkumulator zu addieren

Das ist genau das was ich oben gepostet habe, aber die Rundung habe ich 
tatsächlich vergessen. Mir ist noch nicht klar, ob ich die auch bei der 
Berechnung vom sum verwenden muss oder nur für das Ergebnis.

Ich habe für beide Variablen 32 bit Typen verwendet und eine Division 
anstatt Bitshift, da diese Optimierung Sache des Compilers ist wenn die 
Größe intelligent gewählt ist.

Hier mit Rundung:

int32_t sum = 0; // muss so groß sein, dass 'count * MAX(next) reinpasst'
const int32_t count = 4;
int32_t calcMediumValue(int32_t next) {
 sum = sum - ((sum + count/2) / count) + next; // is the '+count/2' for 5/10 rounding correct here?
 return (sum + count/2) / count; // '+count/2' for 5/10 rounding
}

  Michael

[Bruno V. (bruno_v)]

Maxe schrieb:
> Man hat dann als Rundungsfehler nur noch die Summe der Rundungen im
> Ringspeicher. Ich würde vermuten, dass die sehr klein sind, man täuscht
> sich aber leicht.

Integer-haben beim aufsummieren wie hier eine höhere Auflösung als 
Fließpunkt . Es ist nur mehr Aufwand, Integer richtig zu skalieren (also 
Fixpunkt ohne Überlauf)

> Also, bei 4000 Stützstellen kann im schlechtesten Fall jeweils das
> letzte Bit falsch sein, also summiert ein Wert von 4000 oder 12 Bit. Von
> 32 Bit bleiben dann noch 20 Bit verteilt auf den Eingangsbereich. Dürfte
> gut genug sein.

Bei float waren es nur 10 Bit. Bei Double natürlich mehr, aber wenn 32 
Bit int als Summe nicht reichen, ist es einfach, noch 16 oder 32 dazu zu 
nehmen. Dann ist int genauso gut. Dann nimm aber 4096 Stützstellen, um 
die Division rein durch schieben zu machen.

[Rolf (rolf22)]

Michael D. schrieb:
> Wenn du das mit Integerwerten machst und die Summe der Werte speicherst
> anstatt dem Mittelwert und den Mittelwert erst für die Ausgabe
> verwendest, dann solltest du keine Rundungsfehler haben.

Du hast zwei Integerdivisionen 'sum/count'. Nur in den seltenen Fällen, 
dass 'sum' ein ganzzahliges Vielfaches von 'count' ist, entsteht kein 
Rundungsfehler.

[Steve van de Grens (roehrmond)]

Rolf schrieb:
> Du hast zwei Integerdivisionen 'sum/count'. Nur in den seltenen Fällen,
> dass 'sum' ein ganzzahliges Vielfaches von 'count' ist, entsteht kein
> Rundungsfehler.

Ich denke der Michael meint, dass sich Rundungsfehler nicht immer weiter 
aufaddieren und dadurch zu einem Drift führen (das war ja die 
ursprüngliche Frage).

[Justin S. (Gast)]

Michael D. schrieb:
> Das ist genau das was ich oben gepostet habe

Nein, Du programmierst da etwas völlig anderes. Ich meine folgendes 
(ziemlich das, was Olly T. anfangs andeutete):

int32_t sum = 0;
const int32_t count = 256;
int32_t calcMediumValue(int32_t next) {
    sum -= ringbuffer[idx];
    ringbuffer[idx] = next;
    sum += next;
    if ( --idx < 0 )
        idx = count - 1;
    return (sum + count/2) / count; // '+count/2' for 5/10 rounding
}

initRingbuffer(...)

[Rolf (rolf22)]

Steve van de Grens schrieb:
> Rolf schrieb:
>> Du hast zwei Integerdivisionen 'sum/count'. Nur in den seltenen Fällen,
>> dass 'sum' ein ganzzahliges Vielfaches von 'count' ist, entsteht kein
>> Rundungsfehler.
>
> Ich denke der Michael meint, dass sich Rundungsfehler nicht immer weiter
> aufaddieren und dadurch zu einem Drift führen.

Wenn der Wert des Divisionsergebnisses xxx,5 ist, dann rundet der von 
Compiler generierte Code entweder immer auf xxx oder immer auf xxx+1. 
Dieser Fehler summiert sich in einer Richtung auf, und so entsteht eine 
Drift.

[Hans W. (Firma: Wilhelm.Consulting) (hans-)]

FYI, die übliche Moving Average Implemenierungen sind eigentlich 
"falsch". (Zu lesen: irgendwo vereinfacht)

Siehe 
https://en.m.wikipedia.org/wiki/Moving_average#Cumulative_moving_average

Es geht vor allem um AVG=(n*alter_avg+neuer_wert)/(n+1).

Wenn der AVG driften würde (egal in welche Richtung), dann würde er 
durch diese Integration/summation zu einem stabilen wert konvergieren 
(weil ja n/(n+1) kleiner 1 ist).

Es gibt da zig Abwandlungen, wie man diese Verhältnismäßig komplexe 
Formel annähern kann. Sind irgendwo alles moving-average-filter aber mit 
Vor-/Nachteilen.

73

[Maxe (maxemaxe)]

Max H. schrieb:
> Für das Aufsummieren von Floating Point Werten gibt es den Kahan
> Algorithmus.
> [...]
> Drift wird es immer noch geben, aber wesentlich weniger.

Danke, das werde ich mir auf jeden Fall merken. Allerdings möchte ich 
bei dem Filter wie es nur geht driftfrei sein. Weil es eben sehr schwer 
ist den tatsächlichen Drift einzuschätzen und sich letztlich der Regler 
nicht im Laufe der Zeit verändern darf. Es kommt letztlich nicht so auf 
eine hohe Genauigkleit an, sondern auf langfristige Zuverlässigkeit. 
Wenn der Filter im Laufe eines Jahres um ein paar Prozent abdriftet 
(k.A. ob realistisch), dann verändert sich auch das Reglerverhalten und 
das eher unvorhersehbar.

> Was Du machen könntest wäre ständig eine zweite Instanz des Filters
> aufzubauen. Sobald die 4000 Samples erreicht sind, verwifst Du den
> aktiven Filter und schaltest auf den frischen um.
Ist auch keine schlechte Idee, einfach regelmäßig wieder nachziehen. 
Allerdings gibt das dann u.U. ungleiche Zykluszeiten, was auch wieder 
unvorhersehbare Effekte auslösen kann.

--
Justin S. schrieb:
> Was spricht dagegen, den hinzukommenden Wert (sagen wir 16-Bit) zum
> 32-Bit-Akkumulator zu addieren und in den Buffer zu kopieren, [...]
Dagegen spricht, dass die Werte nicht als Ganzzahlen vorliegen, 
ansonsten ist das ja genau der Algorithmus, den ich verwende.

> (sinnvollerweise
> als Bitshift, bei den Registergrößen machbar bis Filterlänge 2^15)?
Es geht hier um eine SPS, die hat eine FPU und unterstützt auch 
Divisionen.
Und auch die absolute Genauigkeit ist beim Einsatz von Double bzw. 32bit 
Integern kein Problem. Es geht wirklich nur um einen Drift mit 
bleibender Abweichung.

--
Bruno V. schrieb:
> Maxe schrieb:
>> Man hat dann als Rundungsfehler nur noch die Summe der Rundungen im
>> Ringspeicher. Ich würde vermuten, dass die sehr klein sind, man täuscht
>> sich aber leicht.
>
> Integer-haben beim aufsummieren wie hier eine höhere Auflösung als
> Fließpunkt . Es ist nur mehr Aufwand, Integer richtig zu skalieren (also
> Fixpunkt ohne Überlauf)
Wenn ich den Integer-Weg gehe, soll sich das trotzdem über 
Gleitkomma-Parameter in das Programm einfügen. Einfach um die Logik des 
übergeordneten Programms nicht zu verkomplizieren.

[Maxe (maxemaxe)]

Was haltet ihr davon, die letzten paar Stellen der Gleitkommamantisse 
vor dem Reinstecken abzuschneiden? Soviel abschneiden, dass man beim 
"Reinstecken" und Rausholen keinen Rundungsfehler mehr haben kann. Bei 
4000 Stützstellen müsste man wohl 12 Bit abschneiden (wenn ich mich 
nicht täusche). Man könnte bspw. eine eingehende Gleitkommazahl zunächst 
in Single Precision umwandeln, also auf 23 Bit Mantisse und anschließend 
wieder auf Double (52 Bit Mantisse), die letzten 29 Bit wären also 0. 
Ich denke, es dürfte dann zu keinem bleibenden Drift mehr kommen.
Oder übersehe ich was?

[Michael (michaelp)]

Ich habe das schon mal so gelöst:

out = summe_mittelwert + in[0] - in[-4000]
summe_neu = summe_neu + in[0]

alle 4000 Werte:
  summe_mittelwert = summe_neu
  summe_neu = 0

mittelwert = summe_mittelwert / 4000

Damit wird alle 4000 Werte der Rundungsfehler zurückgesetzt.
Die Zykluszeit ist immer nahezu gleich.

[Bruno V. (bruno_v)]

Maxe schrieb:
> Oder übersehe ich was?

Welchen Sinn macht es, fixpoint in gleitkomma nachzubilden.

Gib einfach Deine Wertebereiche an, dann sind konkrete Aussagen möglich. 
Höchstwahrscheinlich ist es in Integer nicht aufwändiger. Und in Double 
zumindestens genauer als in int32. Es summieren sich da auch praktisch 
keine Fehler auf, wenn es in int gepasst hätte.

[Rainer W. (rawi)]

Rolf schrieb:
> Wenn der Wert des Divisionsergebnisses xxx,5 ist, dann rundet der von
> Compiler generierte Code entweder immer auf xxx oder immer auf xxx+1.

Float kennt gar keine xxx,5. Float Zahlen sind Binärzahlen, d.h. 
entweder ist ein Bit gesetzt oder nicht.

Alles andere ist nur ein Thema der Wandlung zur Anzeige als Dezimalzahl

[Rainer W. (rawi)]

Maxe schrieb:
> Jetzt frag ich mich, ob man dabei einen ungewollten Drift reinbekommen
> kann, wenn man mit Gleitkommazahlen arbeitet, weil es beim addieren
> und subtrahieren des gleichen Wertes evtl. zu Rundungen bzw.
> ungleichen Rundungsfehlern kommt.

Zum Addieren und Subtrahieren müssen deine Zahlen sowieso auf den 
gleichen Exponenten normiert werden, d.h. die unteren Bits des 
reinlaufenden und des rausgeschmissenen Wertes gehen verloren, je mehr, 
um so länger dein Array ist. Was also versprichst du dir von der 
Float-Rechnerei?

[Maxe (maxemaxe)]

Michael schrieb:
> Ich habe das schon mal so gelöst:
>
> out = summe_mittelwert + in[0] - in[-4000]
> summe_neu = summe_neu + in[0]
>
> alle 4000 Werte:
>   summe_mittelwert = summe_neu
>   summe_neu = 0
>
> mittelwert = summe_mittelwert / 4000
>
> Damit wird alle 4000 Werte der Rundungsfehler zurückgesetzt.
> Die Zykluszeit ist immer nahezu gleich.

Perfekt! (Deshalb auch das Vollzitat.)

Das ist auch ein transparenter Algorithmus, d.h. man sieht gleich ob die 
eigene Implementierung auch wirklich funktioniert.

--

Bruno V. schrieb:
> Welchen Sinn macht es, fixpoint in gleitkomma nachzubilden.
Naja, im Großen hätte man immernoch den gesamten Gleitkommabereich. Und 
der Drift wäre weg. Aber Michaels Vorgehensweise ist definitiv besser.

> Gib einfach Deine Wertebereiche an, dann sind konkrete Aussagen möglich.
Wertebereich -500 Watt bis 2000 Watt. Nützt doch nichts, soll ja 
universell einsetzbar sein.

--

Rainer W. schrieb:
> Zum Addieren und Subtrahieren müssen deine Zahlen sowieso auf den
> gleichen Exponenten normiert werden, d.h. die unteren Bits des
> reinlaufenden und des rausgeschmissenen Wertes gehen verloren, je mehr,
> um so länger dein Array ist. Was also versprichst du dir von der
> Float-Rechnerei?
Float wird halt im übergeordneten Programm verwendet. Nicht erst extra 
im Gleitenden Mittelwert.

[Hans W. (Firma: Wilhelm.Consulting) (hans-)]

Maxe schrieb:
> Allerdings möchte ich
> bei dem Filter wie es nur geht driftfrei sein.

Sorry, aber was soll da driften?

Wenn du das "richtig" implementierst, dann driftet da nix.

Nochmal, du hast eigentlich

NeuerWert=(AlterWert*N+Neues_Sample)/(N-1)

Das ist eine konvergente Folge!
Die kann dir nicht wesentlich über die 7. Stelle "driften" (bedingt 
durch 32-bit float; da spielen dann die tatsächlichen Werte und die 
größe von N mit rein...).
Da ist auch kein Ringbuffer notwendig.

Mit Ringbuffer ist das dann eine etwas andere Klasse... aber auch da 
kann nix driften.

Maxe schrieb:
> Ein Extrembeispiel wäre, wenn der Durchschnittswert beim "Reinstecken"
> des Wertes so groß ist, dass die kleine zugefügte Zahl komplett
> verschwindet. Also dass das kleinste Bit der Mantisse wertmäßig noch
> größer ist als der ganze zugefügte Wert. Zum Zeitpunkt, wo der Wert
> wieder aus dem Mittelwert fällt, könnte der Durchschnitt klein genug
> sein, dass sich der Wert doch wieder abbildet. Dann hätte man zwischen
> Reinstecken und Rausholen einen negativen Drift.

Das geht in der "richtigen" Implementierung nicht.
Im Code oben, wäre es entweder mit großem AlterWert oder großem N 
möglich, dass in der Addition das neueste Sample "unter geht".
Durch das "/(N-1)", gibt's aber wieder keine drift. "irgendwann" 
(Tiefpass) verringert sich der Wert soweit, dass alles wieder zum Wert 
von Neues_Sample konvergiert.

Mit beliebiger Genauigkeit könntest du natürlich
AlterWert - AlterWert/N + NeuerWert/N rechnen.
Aber genau hier kommt es zum beschriebenen Problem mit Runden.

Das ist aber auch kein "richtiger" Moving Average (siehe Wikipedia 
Link)!

Die mathematischen Hintergründe findet man, wenn man sich IIR und FIR 
Filter ansieht. Über alle Elemente vom Ringbuffer mitteln wäre ein FIR. 
Der ist immer stabil. Die Sache ohne Ringbuffer wäre IIR. Da muss man 
etwas aufpassen.

Wenn du das nicht mit "*N/(N-1)" machst, dann bekommst du einen Pol, der 
sich negativ bemerkbar machen kann.

73

[Purzel H. (hacky)]

Was waere denn der Vorteil eines Gleitenden Averages ueber die letzten 
256 Werte, gegen den exponentiellen Average ueber die letzten 256 Werte 
?

Dass der exponentiellen Average, schneller, einfacher und intuitiver 
ist. Er braucht nur eine speicherstelle.

[Maxe (maxemaxe)]

Hans W. schrieb:
> Nochmal, du hast eigentlich
> NeuerWert=(AlterWert*N+Neues_Sample)/(N-1)

Nein, das ist ein ganz anderer Filter, um den es hier überhaupt nicht 
geht.

> Die mathematischen Hintergründe findet man, wenn man sich IIR und FIR Filter 
ansieht.

Die Frage des Driftens ist m.E. nicht eine der Filtertheorie, sondern 
eine Frage der konkreten Implementierung. Einen Integrator könnte man 
rechnerisch (Assoziativregel) auf zwei parallele aufteilen, und hätte - 
rechnerisch - den gleichen Ausgang. Tatsächlich würden sie dann aber 
auseinanderlaufen und damit in der Praxis nicht funktionieren. (Das nur 
als Beispiel).

Ich würde mich auch überzeugen lassen, dass die Gleitkommazahlen sowieso 
nicht abdriften können, aber dann braucht es schon stichhaltige 
Argumente.

Purzel H. schrieb:
> Was waere denn der Vorteil eines Gleitenden Averages ueber die letzten
> 256 Werte, gegen den exponentiellen Average ueber die letzten 256 Werte
> ?

Was ist denn ein "exponentieller Average über 256 Werte"?

Ein Moving Average ist u.a. dann sinnvoll, wenn man eine einzige 
bekannte Frequenz herausfiltern möchte. Bspw. die Zacken eines 
PWM-Signals lassen sich mit einem passenden Moving Average komplett 
entfernen.

[Hans W. (Firma: Wilhelm.Consulting) (hans-)]

Maxe schrieb:
> Hans W. schrieb:
>> Nochmal, du hast eigentlich
>> NeuerWert=(AlterWert*N+Neues_Sample)/(N-1)
>
> Nein, das ist ein ganz anderer Filter, um den es hier überhaupt nicht
> geht.
>

Stimmt. Aber ein richtiger Moving Average ist eben so definiert!

Was ihr mach ist eine "optimierte" Variante.

Typischerweise machst du das so wie ihr wenn du mit 2^N Basis arbeitest. 
Dann brauchst dir Addition und Schiebeoperationen.

In float bekommst du damit aber ggf Probleme.

Wenn schon Moving Average in floating point, dann implementiere das 
bitte so, wie es eigentlich definiert ist. Damit ist das stabil und 
"driftet" nicht und ist alles in allem ein gutartiger Tiefpass.

Die taktmäßig optimierte Version ist eben nicht unter allen Bedingungen 
optimal.

Der Wikipedia Artikel ist schon recht informativ.

73

[Rainer W. (rawi)]

Hans W. schrieb:
> Nochmal, du hast eigentlich
>
> NeuerWert=(AlterWert*N+Neues_Sample)/(N-1)
>
> Das ist eine konvergente Folge!

Nein, das ist erstens kein gleitender Mittelwert, sondern ein IIR-Filter 
und zweites auch noch falsch. Wenn, dann muss es heißen .../(n+1), sonst 
konvergiert da gar nichts.

Hans W. schrieb:
> Aber ein richtiger Moving Average ist eben so definiert!

Nein, ein Moving Average ist ein FIR-Filter. Der gefilterte Wert ergibt 
sich als Mittelwert der letzten n Samples und alle Sample sind gleich 
gewichtet.
https://en.wikipedia.org/wiki/Moving_average#Simple_moving_average

[Rolf (rolf22)]

Rainer W. schrieb:
> Rolf schrieb:
>> Wenn der Wert des Divisionsergebnisses xxx,5 ist, dann rundet der von
>> Compiler generierte Code entweder immer auf xxx oder immer auf xxx+1.
>
> Float kennt gar keine xxx,5. Float Zahlen sind Binärzahlen, d.h.
> entweder ist ein Bit gesetzt oder nicht.

Ad1: Ich habe auf ein Posting geantwortet, in dem eben nicht mit Float, 
sondern mit Integer gerechnet wurde.

Ad2: Ich habe "Wert" gesagt. Selbstverständlich kann eine Float-Variable 
den Wert xxx,5 haben. Gerade WEIL Exponent und Mantisse Dualzahlen sind.
https://www.h-schmidt.net/FloatConverter/IEEE754de.html

[Hans W. (Firma: Wilhelm.Consulting) (hans-)]

Rainer W. schrieb:
> Wenn, dann muss es heißen .../(n+1)

Uppps. Asche auf mein Haupt!
Natürlich völlig richtig. Der Term mit der Historie muss mit kleiner 1 
gewichtet sein.

Rainer W. schrieb:
> Nein, ein Moving Average ist ein FIR-Filter.

Jain. Es gibt auch den IIR Ansatz... das ist zwar dann eigentlich eher 
ein klassischer Tiefpass, aber es ist auch ein Moving Average.

Aus der Einleitung des verlinkten Artikels:

Variations include: simple, cumulative, or weighted forms (described below).

Jedenfalls kann sowas weder driften, noch oszillieren noch sonst was.

Das was da oben beschrieben wurde sind auch IIR Filter, weil da die 
Historie mit verrechnet wird.

Siehe hier:

Maxe schrieb:
> Michael schrieb:
>> Ich habe das schon mal so gelöst:
>> out = summe_mittelwert + in[0] - in[-4000]
>> summe_neu = summe_neu + in[0]
>> alle 4000 Werte:
>> summe_mittelwert = summe_neu
>> summe_neu = 0
>> mittelwert = summe_mittelwert / 4000
>> Damit wird alle 4000 Werte der Rundungsfehler zurückgesetzt.
>> Die Zykluszeit ist immer nahezu gleich.
>
> Perfekt! (Deshalb auch das Vollzitat.)

Ich will da jetzt nicht irgendwie klugscheißen... Ich wollte nur darauf 
hinweisen, dass ein richtig ausgesuchter Moving Average einfach nicht 
driftet. Wenn er es tut, dann ist es entweder eine ungeeignete Variante, 
oder man hat sich vertippt (wie ich mit dem dummen -/+ Fehler)!

Wie dem auch sei, wenn es die Rechenleistung und der Speicher erlaubt, 
würde ich einfach über alle Elemente in einem Ringbuffer drübermitteln. 
Das macht am wenigsten numerische Fehler/Probleme.

Nebenbei (wenn das Rauschen mitspielt) kann man damit direkt seine 
Auflösung erhöhen... Einfach 2^N Samples addieren und aufs schieben 
vergessen. Das ist aber wieder ein ganz anderes Thema mit seinen eigenen 
Problemen.

73

[Rainer W. (rawi)]

Hans W. schrieb:
> Wie dem auch sei, wenn es die Rechenleistung und der Speicher erlaubt,
> würde ich einfach über alle Elemente in einem Ringbuffer drübermitteln.

Für die nummerischen Fehler ist es hier völlig egal, in welcher 
Reihenfolge man rechnet, solange man die Fehler nicht provoziert.
Jedes Mal n-1 Element noch einmal stumpf aufzusummieren ist so ziemlich 
der dämlichste Ansatz, wenn es doch die Möglichkeit gibt, die Summe 
rekursiv aus der vorherigen zu berechnen. Die letzten n Samples muss man 
bei einem FIR-Filter so oder so im Speicher bereit halten.

> Das macht am wenigsten numerische Fehler/Probleme.
Wieso macht das weniger Fehler als eine rekursive Berechnung der Summe, 
die letztendlich mit den selben Operationen bestimmt wird, nur in 
unterschiedlicher Reihenfolge?

[Hans W. (Firma: Wilhelm.Consulting) (hans-)]

Rainer W. schrieb:
> Wieso macht das weniger Fehler als eine rekursive Berechnung der Summe

Weil bei floats die Reihenfolge der Rechenoperationen eben nicht (ganz) 
egal ist.

Ich bezweifle aber, dass man oft den Fall hat, dass die Input Samples 
sich um  4-5 Größenordnungen unterscheiden... Erst dann wirst du beim 
float, mit seinen 7 signifikanten Stellen, relevante Unterschiede sehen.

73

[Steve van de Grens (roehrmond)]

Rainer W. schrieb:
> Jedes Mal n-1 Element noch einmal stumpf aufzusummieren ist so ziemlich
> der dämlichste Ansatz, wenn es doch die Möglichkeit gibt, die Summe
> rekursiv aus der vorherigen zu berechnen.

Sicher?

Ein paar Integer zu addieren kann durchaus schneller gehen, als die 
alternative (mit Division und so). Vor allem ist es ein geradliniger 
Ansatz, den jeder Programmierer versteht, ohne einen Knoten im Kopf zu 
bekommen. Man soll schließlich wartbaren Code schreiben, damit nach uns 
nicht die Sintflut kommt.

[Hans W. (Firma: Wilhelm.Consulting) (hans-)]

Siehe auch den letzten Teil von
https://www.analog.com/media/en/technical-documentation/dsp-book/dsp_book_ch15.pdf

Wenn du nur ein paar extrem unterschiedliche samples entfernen willst, 
bist du aber ohnehin mit einem Median Filter besser beraten.
Aber auch wieder ein anderes Thema...

73

[Marci W. (marci_w)]

Hallo Steve,

Steve van de Grens schrieb:
> Sicher?
>
> Ein paar Integer zu addieren kann durchaus schneller gehen, als die
> alternative (mit Division und so).

Er meinte vermutlich meine bereits in
Beitrag "Re: Gleitender Mittelwert -> Algorithmus"
dargestellte Variante.

Dividieren musst du an irgend einer Stelle auf jeden Fall, sonst wird es 
schwierig, einen Mittelwert zu berechnen.
Und 4000 Integer zu addieren dauert auch so seine Zeit, bzw. ist eine 
Division schnell erledigt, je nach CPU-Architektur.

ciao

Marci

[Purzel H. (hacky)]

Wenn man die Division so ansetzt, dass duch 2^N geteilt wird, geht's ne 
ganze Ecke schneller.

[Steve van de Grens (roehrmond)]

Marci W. schrieb:
> Und 4000 Integer zu addieren dauert auch so seine Zeit

Das stimmt wohl, das ist weit mehr als "ein paar".

[Rolf (rolf22)]

Marci W. schrieb:
> Und 4000 Integer zu addieren dauert auch so seine Zeit

Früher®, als die Prozessrechner noch viel größer aber viel langsamer 
waren (Minis à la DEC PDP), musste man viel mehr tricksen als heute.

Im konkreten Fall hier hätte man die Addition vielleicht schon begonnen, 
bevor der letzte Summand bekannt war. Dann bekäme man den Mittelwert 
fast sofort nach der letzten Messung.

[Klaus K. (Gast)]

Rolf schrieb:
> Marci W. schrieb:
>> Und 4000 Integer zu addieren dauert auch so seine Zeit
>
> Früher®, als die Prozessrechner noch viel größer aber viel langsamer
> waren (Minis à la DEC PDP), musste man viel mehr tricksen als heute.

Naja, ist halt die Frage wie schnell resp. "echtzeitig" die SPS sein 
muss.
Bei der angesprochenen ca. 50 Jahren alten DEC PDP kann man von 0.3 MIPS 
ausgehen, da sollte immer noch time slots  in Zehntelsekunden 
ermöglichen. Ggf. speichert man Teilsummen als Zwischenwerte.

[Jobst Q. (joquis)]

Maxe schrieb:
> Hallo zusammen, ich brauch für die SPS einen
> Moving-Average-Filter/Gleitenden Mittelwert. Das ist ja eigentlich auch
> leicht und effizient zu programmieren, indem man dem Durchschnitt
> jeweils den neu hinzukommenden Wert anteilig addiert und den
> herausfallenden Wert anteilig subtrahiert.

Kann mir jemand erklären, wozu das gut sein soll, für den gleitenden 
Mittelwert eine Menge Daten zu sammeln, die dann alle gleich gewichtet 
werden?

Ich benutze zum Glätten von Kurven einen einfachen Tiefpass, der außer 
dem Mittelwert selbst nur noch den neuen Wert braucht. Die Differenz des 
neuen Wertes gegenüber dem Mittelwert wird durch eine Konstante 
dividiert und zum  Mittelwert addiert. Fertig ist der Algorithmus.

Dadurch ergibt sich automatisch eine exponentiell abnehmende Gewichtung 
für die älteren Werte, was wiederum zu einer besseren Glättung und 
besserer Aktualität führt.

[Udo S. (urschmitt)]

Maxe schrieb:
> Eine Regelung in der Verfahrenstechnik, das Fluid ist 3 Minuten in der
> Leitung unterwegs, also 180 Sek. Bei 10 Zyklen pro Sekunde hat man 1800
> Werte.

Du schreibst von einer Regelung und einer Totzeit von 180 Sekunden.
Warum nimmst du dann eine so kurze Abtastzeit?
Ich habe vor langer Zeit in der Vorlesung Digitale Regelungstechnik 
gehört dass man die Abtastzeit auf die Geschwindigkeit der Regelstrecke 
und des Reglers anpassen soll. Mit einer so kurzen Abtastzeit bekommst 
du doch auch extrem unhandlich kleine Parameter.
Wenn ich mich recht erinnere war eine gute Faustregel die Abtastzeit 10 
mal kleiner wie die relevanten Zeitkonstanten zu wählen.
Und nicht fast 2000-fach wie bei dir.

[J. S. (engineer)]

Michael D. schrieb:
> Wahrscheinlich ist es dasselbe was J. S. vorgeschlagen hat.

eher nicht, das hier z.B. bringt wieder das Problem zurück:

> sum = sum - ((sum + count/2) / count) + next;

und auch das hier :

> is the '+count/2' for 5/10 rounding correct here?

löst das Problem nicht. Wenn man den Fehler integriert, bekommt man 
irgend ein niederfrequentes Spektrum und das widerspricht dem gesamten 
Ansatz, per Mittelung eine Filterung vorzunehmen und nur die Höhen 
abzuschneiden: ("high cut"):  Die Wahl der Anzahl der Werte, mit denen 
man mittelt, bemisst sich ja anhand einer Frequenzbetrachtung d.h. die 
Grenzfrequenz ist mit der Zahl verknüpft. Das darf man nicht 
konterkarieren. Der Rundungsfehler darf nur lokal auftreten, nach der 
jeweiligen Division.

[Purzel H. (hacky)]

Divisionen soll man vermeiden. Das ist zum Glueck relativ leicht 
moeglich.
- indem man als divisor etwas mit 2^N, zB 2^16 = 65536 waehlt.

[Abdul K. (ehydra)]

Warum sollte man beim Shiften nicht runden müssen?

[Hans W. (Firma: Wilhelm.Consulting) (hans-)]

Nochmal: Das Problem sind die Floats!

Dort sind durch das prinzipbedingte Skalieren die Additionen nicht 
Assoziativ. Damit macht es einen (kleinen) Unterschied, ob du das 1. 
Sample zuerst mit dem Rest addierst und dann abziehst, oder es garnicht 
mitzählst.

Hast du die Möglichkeit die Integer-Samples vom ADC zu verwenden und 
erst nach dem Mittelwert zum Float zu machen, dann hast du diesen Fehler 
nicht.

Das ist in der Analog-Devices Appnote oben beschrieben.

Ansonsten würde ich z.B. auf dspguide.com verweisen... so ein gleitender 
Mittelwert (FIR Filter mit konstanter Gewichtung) hat eigentlich 
vergleichsweise schlechte Eigenschaften... Filterwirkung naja, 
Berechnung naja, Speicherbedarf naja

73

[J. S. (engineer)]

Hans W. schrieb:
> Berechnung naja, Speicherbedarf naja
Bei "Filterwirkung NAJA" gehe ich voll mit, aber ich wüsste nicht, 
welcher gleitende Filter weniger Rechen- und Speicherbedarf hat. Musst 
ihn natürlich nicht als FIR mit ausdrücklichen Werten - sondern als CIC2 
laufen lassen. 2 Akkumulatoren sind das Minimale. Kleiner geht nur ein 
IIR mit Ganzzahlkoeffizienten.