3. Exponential Smoothing Filter in C++ effizient anwenden

Exponential Smoothing Filter in C++ effizient anwenden

  • Y

    Hallo!

    Es handelt sich beim "Exponential Smoothing Filter" um eine Filterroutine mit Exponentialfunktion (wie's der Name schon sagt...).

    Die mathematische Funktion dafür seht ihr hier:
    http://msdn.microsoft.com/en-us/library/jj131429.aspx#ID4EJEAE

    (es geht um die zweite Formel mit dem Summenzeichen)

    Ich habe diese Funktion bereits erfolgreich angewendet, jedoch hantiere ich hier mit etlichen for-Schleifen herum und ich denke, die Code-Übersicht und -Effizient lässt sich noch um einiges steigern... dies ist wichtig, um die "Quasi-" Echtzeit-Darstellung meiner Sensordaten zu gewährleisten.

    Habt ihr vielleicht Ideen zur Verbesserung dieses Codes?:

    #include <vector>
#include <iostream>

using namespace std;
vector<float> data;

float smoothing(float alpha, int n)
{
	float smoothed_data;
	float summe= 0;
	vector<float> dämpfungsfaktor;
	vector<float> zwischenergebnis;

	for (int i = n-1; i > -1; --i) // der exponent i wird mit jedem Schleifendurchgang kleiner, damit der Dämpfungsfaktor für "aktuellere Datenwerte" abnimmt
	{
		dämpfungsfaktor.push_back(pow((1 - alpha), i));
	}

	for (int i = 0; i < dämpfungsfaktor.size(); ++i)
	{
		zwischenergebnis.push_back(dämpfungsfaktor[i] * data[i]); // Elementweise Multiplikation der beiden Vektoren
		summe += zwischenergebnis[i]; // Summe aller bisherigen Elemente dieses Vektors
	}

	smoothed_data = summe * alpha; // Summe wird nochmal mit dem alpha-Faktor multipliziert
	return smoothed_data;
}

void main()
{
	vector<float> smooth_data;
	float alpha = 0.54; // verändert die Dämpfungseigenschaften der "smoothing"-Funktion

	for (int i = 0; i < 40; ++i) // erstelle eine "Sprungfunktion" (bestehend aus Nullen und Einsen) um die "smoothing" Funktion zu testen.
	{
		if (i < 15)
			data.push_back(0);
		else
			data.push_back(1);

		smooth_data.push_back(smoothing(alpha, data.size())); // alpha = 0.54, entspricht einem mittleren Dämpfungsfaktor

		// Zeige die Ergebnisse in der Konsole:
		cout << "Original Data: " << data[i];
		cout << "; Smoothed Data: " << smooth_data[i];
		cout << endl;
	}

	cout << endl;

	system("pause");
}

    Vielen herzlichen Dank im Voraus! 🙂

    Edit, Arcoth: Code-Tags

    0
  • Mod 

    float smoothing(float alpha, int n)
{
	float summe = 0;
	float const base = 1 - alpha;
	float power = 1;

	for (int i = 0; i != n; ++i)
	{
		summe += power * data[n-i-1];
		power *= base;
	}

	return summe * alpha;
}

    ?

    0
  • Mod

    Da kann man noch weitaus mehr optimieren.

    for (int i = 0; i != 40; ++i) // erstelle eine "Sprungfunktion" (bestehend aus Nullen und Einsen) um die "smoothing" Funktion zu testen.
{
	data.push_back( i >= 15 );

	smooth_data.push_back( smoothing(alpha, i + 1) );

	cout << "Original Data: "<< data[i] <<"; Smoothed Data: "<< smooth_data[i] <<'\n';
}

// Mit smooth_data arbeiten?

    Aber: Warum jedes mal smoothing(alpha, i + 1) aufs neue berechnen? Warum nicht die Summe und die Potenz merken und das neuste Element draufaddieren - dann multiplizieren?

    Außerdem: warum nutzt du vector wenn die Größe schon zur Compile-Zeit bekannt ist? Warum nicht einfach ein gewöhnliches Array?

    float data[40];

    Oder, falls nötig, boost::static_vector .

    P.S.: void main ist kein standardkonformes C++. int main .

    0
  • Y

    Wahnsinn! Das steigert die Code-Geschwindigkeit ums etliche... Starke Leistung, das schüttelst du dir einfach so aus den Ärmeln?

    Hut ab!

    warum nutzt du vector wenn die Größe schon zur Compile-Zeit bekannt ist?

    Das ist nur ein kleiner Auszug aus dem eigentlichen Programm. Normalerweise lese ich Sensor-Daten laufend ein und dadurch verändert sich die Größe des Arrays ständig.

    Warum jedes mal smoothing(alpha, i + 1) aufs neue berechnen? Warum nicht die Summe und die Potenz merken und das neuste Element draufaddieren - dann multiplizieren?

    Guter Punkt, das sehe ich mir noch an.

    Vielen Dank für die Hilfe! 👍

    0
  • Y

    Aber: Warum jedes mal smoothing(alpha, i + 1) aufs neue berechnen? Warum nicht die Summe und die Potenz merken und das neuste Element draufaddieren - dann multiplizieren?

    gut, hab ich mir nochmal überlegt... das hat doch wenig Sinn, sich die Summe zu merken? Diese ist ja immer noch den neuen Daten abhängig, da wird nicht einfach was dazuaddiert, sondern die Gewichtung der bereits bestehenden Daten ist dann immer ganz anders. (Der aktuellste Datensatz hat die Gewichtung "1" und wird mit zunehmenden Index immer geringer).

    Was man sich merken könnte wäre die Potenz. Die bleibt gleich, jedoch sehe ich hier keine Möglichkeit den Code zu verbessern?

    power *= base;

    benötigt ja nicht sehr lange.

    Folgender Code zeigt bislang die beste Effizienz (unter 7sek für 1Mill Schleifendurchgänge):

    #include <vector>
#include <iostream>
#include <numeric>

using namespace std;
vector<float> data;

float smoothing(static float alpha, int n)
{
	float summe = 0;
	float const base = 1 - alpha;
	float power = 1;

	for (int i = 0; i != n; ++i)
	{
		summe += power * data[data.size() - i - 1];
		power *= base;
	}

	return summe * alpha;
}

int main()
{
	vector<float> smooth_data;
	static float alpha = 0.54; // verändert die Dämpfungseigenschaften der "smoothing"-Funktion

	for (int i = 0; i < 1000000; ++i) // erstelle eine "Sprungfunktion" (bestehend aus Nullen und Einsen) um die "smoothing" Funktion zu testen.
	{
		data.push_back(i >= 15);

		if (i < 18)
			smooth_data.push_back(smoothing(alpha, i+1)); // alpha = 0.54, entspricht einem mittleren Dämpfungsfaktor
		else
			smooth_data.push_back(smoothing(alpha, 19));

		//cout << "Original Data: " << data[i] << "; Smoothed Data: " << smooth_data[i] << '\n';

}

	cout << "Fertig!" << '\n';

	system("pause");
	return 0;
}
    0
  • ?

    Das ist übrigens ein simpler 1st order IIR Filter, wobei man jedes Ausgabe-Sample in O(1) berechnen kann. Ausgehend von der Formel auf der verlinkten Seite kann man das nämlich "rekursiv" vereinfachen:

    X^_n=α∑_i=0n(1−α)iXn−i=α(X_n+∑_i=1n(1−α)iXn−i)=αX_n+α∑_i=1n(1−α)iXn−i=αX_n+(1−α)X^_n−1\begin{array}{lcl} \hat{X}\_n & = & \alpha \sum\_{i=0}^n \left( 1 - \alpha \right)^i X_{n-i} \\ & = & \alpha \left( X\_n + \sum\_{i=1}^n \left( 1 - \alpha \right)^i X_{n-i} \right) \\ & = & \alpha X\_n + \alpha \sum\_{i=1}^n \left( 1 - \alpha \right)^i X_{n-i} \\ & = & \alpha X\_n + \left( 1 - \alpha \right) \hat{X}\_{n-1} \end{array}

    Daraus kann man sich ein Funktionsobjekt mit Parameter alpha und Zustand bauen:

    struct exp_smoothing {
    float alpha; // smoothing parameter (1 = identity mapping)
    float last;  // state

    float operator()(float x) { 
        float curr = alpha * x + (1 - alpha) * last;
        last = curr;
        return curr;
    }
};

    Und falls euch noch die Übertragungsfunktion davon interressiert, hier ist sie:

    H(z)=α1−(1−α)z−1H(z) = \frac{\alpha}{1 - \left( 1 - \alpha \right) z^{-1} }

    Diese Funktion liefert Euch zu einer bestimmten Frequenz ff mit 0≤f≤π0 \leq f \leq \pi und z=ei⋅fz = e^{i \cdot f} die Amplitudenübertragung ∣H(z)∣|H(z)| und die Phasenantwort arg(H(z))\arg \left(H(z) \right) des obigen Filters.

    Herzlich willkommen in der Welt der digitalen Signalverarbeitung!
    😃

    Wenn ihr Euch da mehr einlesen wollt, schaut mal auf dspguide.com nach.

    0
  • Y

    Ja, diese digitale Signalverarbeitung interssiert mich... damit möchte ich mich beschäftigen!

    Kannst du mir noch einen Hinweis, Links bez. Tipps geben diese Funktion innerhalb des Structs auch anwenden zu können? Das habe ich noch nie gemacht und stehe hier etwas an.

    Fertige Lösungen sind gar nicht nötig, würde gern selbst draufkommen 😉

    Danke

    0
  • ?

    oder auch kürzer und mit einer Multiplikation weniger:

    float operator()(float x) {
     return last += (x - last) * alpha; 
 }
    0
  • ?

    YviXbit schrieb:

    Kannst du mir noch einen Hinweis, Links bez. Tipps geben diese Funktion innerhalb des Structs auch anwenden zu können? Das habe ich noch nie gemacht und stehe hier etwas an.

    In etwa so.

    0
  • Y

    kkaw schrieb:

    oder auch kürzer und mit einer Multiplikation weniger:

    float operator()(float x) {
     return last += (x - last) * alpha; 
 }

    ich denke, kürzer geht's nun wirklich nicht mehr 😉

    Hab's geschafft:

    #include <vector>
#include <iostream>

using namespace std;
vector<float> data;

struct exp_smoothing {
	float alpha; // smoothing parameter (1 = identity mapping) 
	float last;  // state 

	float operator()(float x) {
		return last += (x - last) * alpha;
	}
};

int main()
{
	vector<float> smooth_data;

	exp_smoothing DataToSmooth;

	DataToSmooth.alpha = 0.54; // alpha = 0.54, entspricht einem mittleren Dämpfungsfaktor
	DataToSmooth.last = 0; // Initialisierung

	for (int i = 0; i < 1000000; ++i) // erstelle eine "Sprungfunktion" (bestehend aus Nullen und Einsen) um die "smoothing" Funktion zu testen.
	{
		data.push_back(i >= 15);

		smooth_data.push_back(DataToSmooth.operator()(data[i]));	
		//cout << "Original Data: " << data[i] << "; Smoothed Data: " << smooth_data[i] << '\n';
	}

	cout << "Fertig!" << '\n';

	system("pause");
	return 0;
}

    Das verdoppelte die Geschwindigkeit nochmal... Unter 3 Sekunden für 1Mill Schleifendurchgänge! Sehr stark. Vielen, vielen Dank euch allen!!

    0
  • V

    Schätze, statt

    smooth_data.push_back(DataToSmooth.operator()(data[i]));

    war die andere Schreibweise

    smooth_data.push_back(DataToSmooth(data[i]));

    gemeint.

    0
  • ?

    Unter 3 Sekunden für 1Mill Schleifendurchgänge!

    Du misst aber dabei noch die Zeit fuer das Erzeugen Deiner Eingabedaten mit!

    Besser so:

    #include <vector>
#include <iostream>
#include <chrono>

using namespace std;

struct ExpSmoothing 
{
	typedef float dataType;

	dataType operator()(dataType x) 
	{
        return last += (x - last) * alpha;
    }

    dataType alpha; // smoothing parameter (1 = identity mapping)
    dataType last;  // state
};

void createInputData(size_t sz, vector<ExpSmoothing::dataType>& data)
{
    for (size_t i = 0; i < sz; ++i) // erstelle eine "Sprungfunktion" (bestehend aus Nullen und Einsen) um die "smoothing" Funktion zu testen.
    {
        data.emplace_back(i >= 15);
	}
}

void filterData(const vector<ExpSmoothing::dataType> input, 
				vector<ExpSmoothing::dataType>& output, 
				ExpSmoothing& filter)
{
    for (size_t i = 0; i < input.size(); ++i) 
    {
        output.push_back(filter(input[i]));   
    }
}

int main()
{
	size_t loopCount = 1000000;

	vector<ExpSmoothing::dataType> inputData;
    vector<ExpSmoothing::dataType> smoothedData;

    ExpSmoothing filter;

    filter.alpha = ExpSmoothing::dataType(0.54); // alpha = 0.54, entspricht einem mittleren Dämpfungsfaktor
    filter.last = ExpSmoothing::dataType(0); // Initialisierung

	createInputData(loopCount, inputData);

	auto start = chrono::high_resolution_clock::now();
	filterData(inputData, smoothedData, filter);
    auto exp = chrono::duration_cast<std::chrono::milliseconds>(chrono::high_resolution_clock::now()-start).count();
    cout << exp << " ms" << '\n';

    system("pause");
    return 0;
}
    0
Anmelden zum Antworten