Decay-Filter

1. Definition

Ein Filter berechnet einen Ausgabewert y_s aus einem aktuellen Eingabewert x_s und einem internen Zustand z_s, der von den Eingabewerten der Vergangenheit abhängt. Mathematisch ist dies eine Faltung, allerdings rechnen wir nicht kontinuierlich, sondern mit diskreten Schritten s ∊ {0, 1, …}.

                 +------------------+
                 |                  |
                 |      Filter      |
Ausgabe ys <--<--+                  +-<--<--  Eingabe xs
                 |    Zustand zs    |
                 |                  |
                 +------------------+

Beim Decay-Filter beeinflusst ein Eingabewert die gesamte Zukunft, aber umso schwächer, je weiter ein Zeitpunkt entfernt ist. Der Filter wird durch einen Paramter k definiert mit 0 ≤ k < 1: bei einem großen k bleibt die Wirkung eines Eingabewertes lange erhalten, bei einem kleinen k klingt sie schneller ab. Es gibt eine weitere Konstante k'; diese ergibt sich aber, wie wir noch sehen werden, eindeutig aus k.

Die Ausgabe unseres Filter ist:

⟪y_s := k' * (1 * x_s + k x_(s-1) + k^2 x_(s-2) + k^3 x_(s-3) + ...)⟫
⟪y_s := k' * sum_(i=0)^oo k^i * x_(s-i)⟫

(1.1) ⟪y_s := sum_(i=0)^oo k^i * k' * x_(s-i)⟫

2. Mathematik

Eine Forderung an unseren Filter ist, dass bei konstanter Eingabe ⟪x_s := x⟫ nach hinreichender Zeit die Ausgabe gleich der Eingabe ist:

⟪x = k' * (1 * x + k x + k^2 x + k^3 x + ...)⟫
⟪x = k' * sum_(i=0)^oo k^i x⟫
⟪1 = k' * sum_(i=0)^oo k^i ⟫
⟪1/(k')= sum_(i=0)^oo k^i ⟫

Für ⟪|k| lt 1⟫ ergibt sich mit (M.5):

⟪1/(k')= 1 / (1-k)⟫

(2.1) ⟪ k' = (1-k) ⟫

Die Gleichung enthält alle Eingangswerte der Vergangenheit x_s-i. Wir können sie aber in eine inkrementelle Darstellung umformen, die zur Berechnung von y_s nur x_s und y_s-1 braucht:

⟪y_s = sum_(i=0)^oo k^ i * k' * x_( s - i ) ⟫
⟪y_s = k^0 * k' * x_s + sum_(i=1)^oo k^ i * k' * x_( s - i ) ⟫
⟪y_s = k^0 * k' * x_s + sum_(i=0)^oo k^(i+1) * k' * x_( s - (i+1) ) ⟫
⟪y_s = k^0 * k' * x_s + k * ubrace(sum_(i=0)^oo k^ i * k' * x_( (s-1) - i ))_(y_(s-1))⟫

Mit (1.1) und (2.1) erhalten wir die inkrementelle Darstellung:

(2.2) ⟪y_s = (1-k) * x_s + k * y_(s-1) = k * y_(s-1) + (1-k) * x_s⟫

3. Implementierung mit Ganzzahl-Arithmetik

Die Filtergleichung nutzt die Werte k und 1-k, beide zwischen 0 und 1. Wir skalieren mit 2³¹ und erhalten den Wertebereich von 0 bis 2³¹. Durch Verzicht auf die beiden Extremwerte erhalten wir für k_s und 1-k_s den Bereich 1 bis 2³¹-1, dessen Werte sich als int32_t darstellen lassen.

static volatile	int32_t	ks;

Wir formen (2.2) um:

⟪y_t = k * y_(s-1) + (1-k) * c_t⟫

⟪y_t = (2^31 k)/2^31 * y_(s-1) + (2^31-2^31 k)/2^31 * c_t⟫

⟪y_t = 1/2^31 * ( 2^31 k * y_(s-1) + (2^31-2^31 k) * c_t )⟫

(3.1) ⟪y_t = 1/2^31 * ( k_s * y_(s-1) + (2^31-k_s) * c_t )⟫

Für Eingabe x_s und Ausgabe y_s nutzen wir den Wertebereich von int32_t:

static int32_t out;

Die beiden Multiplikationen haben int32_t Argumente, aber ein int64_t Ergebnis:

static inline int64_t
mul32_64( int32_t a, int32_t b ){

	return (int64_t) a * b;
}

Die Summe der beiden Produkte passt in ein int64_t.

Die Division duch 2³¹ kann durch ein Shift implementiert werden. Bei negativer Summe muss aber der Absolutbetrag geshifted werden, weil ein Shift negativer Werte von 0 weg rundet und das einen Überlauf verursachen kann.

void
decay_filter( int32_t in, int32_t ks, int32_t *outp ){

	// if( ks < 1 ){ ks = 1; }

	const int32_t minus_ks = INT32_MAX - ks + 1;

	const int64_t a = mul32_64( *outp, ks       );
	const int64_t b = mul32_64(  in  , minus_ks );

	const int64_t sum = a + b;

	if( sum >= 0 ){

		const uint64_t abs_sum = (uint64_t) +sum;
		const uint64_t shifted = abs_sum >> 31;

		*outp = +( (int32_t) shifted );
	} else {
		const uint64_t abs_sum = (uint64_t) -sum;
		const uint64_t shifted = abs_sum >> 31;

		*outp = -( (int32_t) shifted );
	}
}

4. Reale Implementierung

void
decay_filter( int32_t in, int32_t ks, int32_t *outp ){

	const register int64_t sum =
		mul32_64( *outp,             ks ) +
		mul32_64(  in  , INT32_MIN - ks );

	if( sum >= 0 ){

		*outp = +(int32_t)( (+sum) >> 31 );
	} else {
		*outp = -(int32_t)( (-sum) >> 31 );
	}
}

decay_filter( in, ks, &out );

Download: decay-filter.cc

Die Funktion benutzt:

• 1 32-Bit Signed Subtraktion,
• 2 Signed 32-Bit Multiplikationen mit 64-Bit-Ergebnis,
• 1 Signed 64-Bit Addition,
• 1 64-Bit Negation,
• 1 Unsigned 64-Bit Shift, und
• 1 32-Bit Negation.

5. Berechnung von k aus der Zeitkonstante

Wir wollen k aus der Halbswertzeit t_h bestimmen. Nach s_h Schritten soll die Wirkung einer Eingabe auf 0.5 abgeklungen sein:

⟪ k^(s_h) = 0.5 ⟫
⟪ k = 0.5 ^(1/s_h) ⟫
⟪ k = (1//2)^(1/s_h) ⟫

(5.1) ⟪ k = 1//2 ^(1/s_h) ⟫

Wenn die Schritte mit einer konstanten Frequenz f erfolgen, werden in der Zeit t_h die Anzahl f * t_h Schritte durchlaufen:

⟪ s_h = f * t_h ⟫

Zusammen mit (5.1) erhalten wir:

(5.2) ⟪ k = 1//2 ^(1/(f * t_h)) ⟫

Zur Vorbereitung der Implementierung wird die Formel in Schritte zerteilt:

(5.3) ⟪ c := 1/(f * t_h) ⟫

(5.4) ⟪ d := 2^c ⟫

(5.5) ⟪ k := 1 // d ⟫

• Die Frequenz f, mit der die Filterfunktion aufgerufen wird, wird in Hertz (1/Sekunde) angegeben.
• Um mit ganzen Zahlen rechnen zu können, geben wir die Halbwertszeit t_h in Mikrosekunden an.

Damit erhalten wir:

(5.3a) ⟪ c := 10^6 / (f^([Hz]) * t_h^([mu s])) ⟫

6. Implementierung mit Ganzzahl-Arithmetik

Zur Berechnung der Potenz nutzen wir die Funktion exp2m1(x) := 2^x-1, eine Verschmelzung aus exp2 und expm1.

Die Ganzzahl-Implementierung u32p32exp2m1(x), wir nennen sie ab hier p(x), hat als Definitionsbereich und Wertebereich 0 ... 1-2^-32, wobei sowohl Argument als auch Ergebnis mit 2³² skaliert werden:

(6.1) ⟪ p( 2^32 x ) = 2^32 * (2^x-1), \ 2^32 x in tt"uint32_t" wedge p( 2^32 x ) in tt"uint32_t" ⟫

Für die Nutzung dieser Funktion führen wir die skalierten ganzzahligen Variablen c_s, d_s und k_s ein.

Aus (5.3a) wird c_s definiert:

(6.2) ⟪ c_s := 2^32 c = min{ (2^32 * 10^6) / (f^([Hz]) * t_h^([mu s])), 2^32-1} in tt"uint32_t" ⟫

Der verfügbare Wertebereich für c_s erfordert s_h ≥ 1+1/(2³²-1), Die Filterroutine muss also etwas häufiger als einmal je Sekunde aufgerufen werden. Diese Bedingung wird trivial erfüllt.

Aus (5.4) und (6.1) wird d_s definiert:

⟪ d-1 = 2^c - 1 ⟫
⟪ 2^32 (d-1) = 2^32 ( 2^c - 1) = p( 2^32 c ) = p( c_s ) ⟫

(6.3) ⟪ d_s := 2^32 (d-1) = p( c_s ) in tt"uint32_t" ⟫

Aus (5.5) wird k_s definiert:

⟪ 2^31 k = 2^31/d = 2^63 / ( 2^32 d ) = 2^63 / ( underbrace(2^32 (d-1))_(d_s) + 2^32 ) = 2^63/(d_s + 2^32) ⟫

(6.4) ⟪ k_s := 2^31 k = min{ 2^63/(d_s + 2^32), 2^31-1} in tt"int32_t" ⟫

Der verfügbare Wertebereich für k_s erfordert d_s ≥ 2, mit (6.3) also d ≥ 1+2^-31, mit (5.4) also 2^c ≥ 1+2^-31 oder umgeformt c ≥ log₂(1+2^-31). Wegen (6.2) ist das erfüllt für c_s = 2³²c ≥ 2³² log₂(1+2^-31) = 2.89, also c_s ≥ 3, mit c_s = 2³² / s_h demnach für s_h ≤ 2³² / 3 ~ 1.4·10⁹. Die Filterroutine darf also maximal 1.4 Milliarden mal je Sekunde aufgerufen werden. Auch diese Bedingung wird trivial erfüllt.

Damit haben wir unsere Funktion:

#include <stdint.h>

// extern uint32_t u32p32exp2m1( uint32_t );
#include "u32p32exp2m1.h"

int32_t
compute_ks( uint32_t f_hz, uint32_t th_us ){

	const uint64_t nominator = (1000ul * 1000ul) << 32;
	const uint64_t cs64 = nominator / f_hz / th_us;
	const uint32_t cs   = cs64 <= UINT32_MAX ? (uint32_t) cs64 : UINT32_MAX;
	const uint32_t ps   = u32p32exp2m1( cs );

	const uint64_t nom  = 1ul << 63;
	const uint64_t den  = ps + ( 1ul << 32 );

	const uint64_t ks64 = nom / den;
	const  int32_t ks   = ks64 <=  INT32_MAX ? ( int32_t) ks64 :  INT32_MAX;

	return ks;
}

Download: decay-compute-ks.cc

Wir nutzen (die Konversionen erzeugen keinen Code und die konstanten Ausdrücke optimiert der Compiler weg):

• 1 Addition von oder Oderieren mit einer konstanten ganzzahligen Potenz von 2;
• 2 Triviale Vergleiche mit Auswahl;
• 3 Unsigned 64-Bit-Divisionen;
• 1 Aufruf von u32p32exp2m1().

Divisionen und besonders der Aufruf von u32p32exp2m1() dominieren den Aufwand.

Wenn garantiert ist, dass das Produkt aus f_hz und th_us kleiner als UINT32_MAX ist, kann in der zweiten Zeile eine 64-Bit-Division eingespart werden:

	const uint64_t temp = nominator / ( f_hz * th_us );