Standardabweichung auf embedded Systemen schnell berechnen

19.10.2021

Im Folgenden möchte ich eine Methode vorstellen, wie man die Standardabweichung näherungsweise mit sehr geringem Ressourcenaufwand kontinuierlich berechnen kann. Dabei möchte ich nicht auf die mathematischen Hintergründe eingehen, sondern das Prinzip erklären und die optimale Lösung präsentieren.

Die Standardabweichung „std()“ einer Messreihe ist definiert als Quadratwurzel „sqrt()“ der Varianz „var()“ der Messreihe. Die Varianz wiederum ist die mittlere quadratische Abweichung vom Mittelwert „avg()“ der Messreihe.

std( x ) := sqrt( var( x )) := sqrt( avg( (x - avg( x ) )^2 ) )

Da auf kleinen Mikrocontrollern der verfügbare Speicher und die verfügbare Rechenzeit eigentlich immer sehr knapp sind, bringt die Umsetzung nach dieser Definition hauptsächlich zwei Herausforderungen mit sich:

Für die Berechnung der Wurzel bietet sich die „sqrt()“ Funktion aus der Standardbibliothek an. Diese kostet allerdings relativ viel Laufzeit und Speicher, darum vereinfachen wir die Berechnung und lass sie lieber gleich ganz weg. Dadurch entfällt auch das Quadrieren und wird durch den Betrag „abs()“ ersetzt.

~std( x ) := avg( abs( x - avg( x ) ) )

Die zweite Herausforderung bringt die doppelte Mittelwertberechnung. Hierfür müssten alle Messwerte aus einem Messfenster im Speicher gehalten werden, und für jede Ausgabe müssten alle Differenzen neu berechnet, aufsummiert und am Ende dividiert werden. Welche Funktion bietet näherungsweise einen Mittelwert und lässt sich ohne diese Nachteile berechnen? Ein Tiefpassfilter, am einfachsten ein Tiefpassfilter erster Ordnung „pt1()“.

~std( x ) := pt1( abs( x – pt1( x ) ) )

Der Tiefpass erster Ordnung lässt sich für zeitdiskrete Abtastwerte mit nur einer Multiplikation, einer Addition und einer Subtraktion implementieren.

pt1( xn ) = xn-1 + a*(un - xn-1)

Mit der Konstante „a“ wird quasi die Fensterbreite eingestellt, mit einem Wertebereich von ]0:1[. Je kleiner der Wert, desto breiter das Fenster. Es ist auch die Einschwingzeit des Filters zu beachten, welche bei ca. 10/a Abtastwerten liegt.

Eine einfache Umsetzung in C könnte so aussehen:

#include math.h
float std(float x) {
  const float a = 0.1;
  static float xn = 0.0;
  static float xm = 0.0;
  float tmp;
  xn = xn + a * (x - xn);    // pt1
  tmp = abs(x - xn);
  xm = xm + a * (tmp - xm);  // pt1
  return xm;
}

Da nicht alle Mikrocontroller mit Fließkommazahlen (float) umgehen können und eine Fließkommaberechnung eigentlich immer langsamer ist als eine Rechnung mit ganzen Zahlen, lässt sich der Algorithmus auch super als Festkommaberechnung umsetzen. Durch geschickte Wahl von „a“ kann man dann auch gleich noch auf die Multiplikationen verzichten. Dabei gibt es natürlich Einschränkungen im Wertebereich und der Auflösung. Auf einem 32-bit System, mit einer Auflösung von 2^-10 und einem Wertebereich von [-4096:4096] mit a = 2^-10 ließe sich der Algorithmus wie folgt optimieren:

#include math.h
int  std(int  x /* 2^-10 */) {
  // const int a = 1;      // 2^-10
  static int  xn = 0;      // 2^-20
  static int  xm = 0;      // 2^-20
  int  tmp;                // 2^-10
  xn = xn + x - (xn >> 10);
  tmp = abs(x - (xn >> 10));
  xm = xm + x - (xm >> 10);
  return  xm >> 10;        // 2^-10
}

Auf einem ARM Cortex M3 benötigt die optimierte Umsetzung 8 Byte RAM, 38 Byte ROM und ca. 17 Taktzyklen pro Messwert.

std:                // r0 := x
movw  r1, #0
movt  r1, #8192     // base address
ldr   r3, [r1, #0]  // load xn
asr   r2, r3, #10   // xn >> 10
sub   r2, r0, r2    // tmp = a * (x - xn)
add   r3, r3, r2    // xn += tmp
str   r3, [r1, #0]  // store xn
asr   r2, r3, #10   // xn >> 10
subs  r0, r0, r2    // tmp = x - xn
rsbmi r0, r0, #0    // negate if negative
ldr   r3, [r1, #4]  // load xm
asr   r2, r3, #10   // xm >> 10
sub   r2, r0, r2    // tmp = a * (tmp - xm)
add   r3, r3, r2    // xm += tmp
str   r3, [r1, #4]  // store xm
asr   r0, r3, #10   // xm >> 10
mov   pc, lr        // jump back

Eventuell sind die Eingangswerte auf den gültigen Wertebereich zu überprüfen, falls dies nicht garantiert werden kann.

In meinen Simulationen wird die Standardabweichung durch den vorgestellten Algorithmus um konstant ca. 20% unterschätzt. Da der Offset konstant ist, kann er nachträglich korrigiert werden.

Back to top