Curve fit algorihtmuss

[BorisProgger]

Habe ein Array mit x und Y Koordinaten einer grafischen Kurve (sinus)
nun will ich anhand der Koordnaten rausfinden was für ne Kurve das ist?.. gibt es da spezielle Curve-fit algorithmen ??

[Phoemuex]

Ich weiß nicht genau, was du suchst, aber wenn (sinus) bedeutet, dass es eine Summe aus Sinus- und Kosinustermen ist, dann hilft dir das hier wahrscheinlich weiter:

http://de.wikipedia.org/wiki/Trigonometrische_Interpolation

Dazu habe ich auch mal ein Programm geschrieben, dass aus der Schwingung einer Wave-Datei die Schwingung berechnet. Wenn es das richtige ist, kann ich dir den Quellcode mal posten.

Felix

[BorisProgger]

ja das wäre schon mal super... wenn mir den code und vll. paar kommentakre zukommen lassen könntest;) danke

Das ist dann quasei nichte anderes wie ne Foruier analyse die du da hast oder? aber her damit..

[Phoemuex]

#include <iostream>
#include <vector>
#include <string>
#include <cmath>
#include <fstream>

#include "wav.h"

using namespace std;

/**Integral(anfang, ende, abstand) berechnet die Obersumme des Integrals
  *der Funktion, deren Punkte Äquidistant mit dem Abstand abstand
  *gemessen wurden und in dem Speicher, der durch anfang und ende
  *begrenzt wird, gespeichert sind*/
template<class T>
double Integral (T anfang, T ende, double abstand)
{
	//Es werden erst alle Werte aufsummiert
    double summe = 0;
    ++anfang;
    for(; anfang != ende; ++anfang)
    {
        summe += *anfang;
    }
	//und mit dem Abstand multipliziert
    return summe * abstand;
}

/** FourierCoefficients(anfang, ende, abstand, n, xAnfang) berechnet die Fourierkoeffizienten
  *(a_0,b_0)..(a_n,b_n) der Funktion, deren Punkte Äquidistant mit dem Abstand
  *abstand gemessen wurden und in dem Speicher, der durch anfang und ende begrenzt
  *wird, gespeichert sind. Dabei wird angenommen, dass die Schwingungsdauer der
  *Funktion gleich der Länge des gegebenen Intervalls ist.
*/
vector<pair<double,double> > FourierCoefficients (vector<int>::const_iterator anfang,
 vector<int>::const_iterator ende, const double abstand, const int n, const double xAnfang = 0.0)
{
	/*Zuerst wird die Länge des Intervalls und damit die Schwingungsdauer bestimmt
	  Dazu wird der Abstand mit der Anzahl der Werte multipliziert*/
    double laenge = 0;
    {
        int i = 0;
        for(vector<int>::const_iterator anf = anfang; anf != ende; ++anf)
        {
            ++i;
        }
        laenge = abstand * i;
    }
	//Der Vektor (Array mit dynamischer Länge) mit den Rückgabewerten
    vector<pair<double,double> > returnVector;
	/*Zuerst werden a_0 und b_0 hinzugefügt, wobei b_0 = 0 ist.
	  a_0 wird über die normale Formel mit Hilfe der Funktion Integral berechnet*/
    returnVector.push_back(pair<double,double>((1.0/laenge)*Integral(anfang, ende, abstand), 0.0));
    for(int i = 1; i < n; ++i)
    {
		//Berechnung von a_i
        double a = 0;
		//Vector von neuen Funktionswerten der Funktion cos(2*Pi/T*i*x)*f(x)
        vector<double> werte;
        double k = xAnfang;
		//Füllen des Vektors mit den zugehörigen Werten
        for(vector<int>::const_iterator anf = anfang; anf != ende; ++anf)
        {
            werte.push_back(*anf * cos((2*M_PI/laenge)*i*k));
            k += abstand;
        }
		//Berechnen des Koeffizienten a_i mit Hilfe der üblichen Formel
        a = (2/laenge)*Integral(werte.begin(), werte.end(), abstand);
        //Berechnung von b_i
        double b = 0;
		//Löschen der alten Werte des Vektors
        werte.clear();
        k = xAnfang;
		//Füllen des Vektors mit den zugehörigen Werten der Funktion sin(2*Pi/T*i*x)*f(x)
        for(vector<int>::const_iterator anf = anfang; anf != ende; ++anf)
        {
            werte.push_back(*anf * sin((2*M_PI/laenge)*i*k));
            k += abstand;
        }
		//Berechnen des Koeffizienten b_i mit Hilfe der üblichen Formel
        b = (2/laenge)*Integral(werte.begin(), werte.end(), abstand);
		//Zufügen von a_i und b_i zum Vektor, der am Ende zurückgegeben wird
        returnVector.push_back(pair<double,double>(a,b));
    }
    return returnVector;
}

int main(void)
{
    //PointNumber gibt die Anzahl der zu lesenden Werte an
    const int PointNumber = 1000;
	string name;
	cout << "Dieses Programm berechnet die Fourierkoeffizienten der in einer ";
	cout << "Wav-Datei gespeicherten Schwingung." << endl << endl;
	cout << "Bitte geben Sie den Dateiname der Wav-Datei ein: ";
	cin >> name;
	WavDatei File(name);
	//Datei zur Ausgabe der Fourierkoeffizienten
	fstream Datei("Koeffizienten.txt", ios::out);
	//Datei zur Ausgabe der aus der Wav-Datei gelesenen Werte
	fstream Datei2("Werte.txt", ios::out);
	//In dem Vector werte werden die gelesenen Werte gespeichert
	vector<int> werte;
	cout << endl << "Daten der Wav-Datei:" << endl;
	cout << "Channels: \t" << File.GetChannels() << endl;
	cout << "Sampling rate: \t"<< File.GetSamplingRate() << endl;
	cout << "Bits per sample:" << File.GetBpS() << endl << endl;
	/*Am Anfang der Wav-Datei befinden sich einige ungütlige Werte,
	diese werden hiermit übersprungen*/
	for(int i = 0; i < 60; ++i)
	{
	    File.ReadNextValue();
	}
	for(int i = 0; i < PointNumber; ++i)
	{
	    werte.push_back(File.ReadNextValue());
	}
	//Berechnung des Mittelwertes der Schwingung
    double mittelwert = 0.0;
	for(vector<int>::const_iterator it = werte.begin(); it != werte.end(); ++it)
	{
	    mittelwert += *it;
	}
	mittelwert /= werte.size();
	vector<int>::const_iterator anfangsIterator = werte.begin();
	/*Der folgende Block stellt sicher, dass der Anfangs- und End-Wert der
	gelesenen Schwingung identisch ist.*/
	/* Es werden solange Werte übersprungen, bis der Anfang des Intervalls so gewählt ist,
	   dass die Funktion am Anfang den Mittelwert annimmt und eine positive Steigung hat*/
	while(*anfangsIterator > mittelwert)
	{
	    ++anfangsIterator;
	}
	while(*anfangsIterator < mittelwert)
	{
	    ++anfangsIterator;
	}
	vector<int>::const_iterator endIterator = werte.end();
	--endIterator;
	/* Es werden solange Werte übersprungen, bis das Ende des Intervalls so gewählt ist,
	   dass die Funktion am Ende den Mittelwert annimmt und eine positive Steigung hat*/
	while(*endIterator < mittelwert)
	{
	    --endIterator;
	}
	while(*endIterator > mittelwert)
	{
	    --endIterator;
	}
	//Ende des Blocks, der Gleichheit des Anfangs- und Endwerts sicherstellt
	//Berechnung der Fourier-Koeffizienten
	vector<pair<double,double> > koeffizienten = FourierCoefficients(anfangsIterator, endIterator, 1.0/File.GetSamplingRate(), 400);
	//Schreiben der gelesenen Schwingung und der Koeffizienten ind verschiedene Dateien
	for(vector<pair<double,double> >::const_iterator it = koeffizienten.begin(); it != koeffizienten.end(); ++it)
	{
	    Datei << it->first << " ";
	    Datei << it->second << " ";
	}
	for(vector<int>::const_iterator it = anfangsIterator; it != endIterator; ++it)
	{
	    Datei2 << *it << " ";
	}
	cout << "Fuer eine graphische Auswerung starten Sie bitte das Maple-Worksheet \"Auswertung.mws\"";
	cout << ", passen gegebenenfalls die Samplingrate an und klicken auf die drei Ausrufungszeichen.";
	//Die beiden folgenden Zeilen verhindern, dass sich das Fenster des Programms direkt wieder schließt
	cin.get();
	cin.get();
	return 0;
}

So, das ist der Code, den ich geschreiben habe, 1 zu 1 kopiert.

Was fehlt sind die Header- und die Implementierungsdateien der Klasse WaveDatei, aber die brauchst du ja nicht unbedingt, weil du ja wahrscheinlich andere Punkte und keine WaveDateien hast.

So wie da oben funktioniert das aber nur mit äquidistant gemessenen Punkten, d. h. der Abstand zwischen zwei gemessenen Punkten muss immer gleich sein. Wenn das nicht der Fall ist, musst du es ein bisschen ändern.

Ich hoffe, es hilft dir.

Felix

[BorisProgger]

herlichen dank FElix.. ich werde es mal aubrobieren:)

[BorisProgger]

Problem an FElix::

///////////////////////// WERTE EINER EIFNACHE SINUS KURVE ///////////////////
	int i=0,r=50;
	POINT *pWerte= new POINT[NUMBER];
	for(double d= 0; d<=(M_PI*2);d+=(M_PI*2)/NUMBER,i++){
		pWerte[i].y=sin(d)*r;
		pWerte[i].x=i+1;
	}
	//////////////////////////////////////////////////////////////////////////////

	/// WERTE ÜBERNAME ////
	for(int i = 0; i < PointNumber; ++i) 
    { 
        werte.push_back(pWerte[i].y); 
    } 
	delete[] pWerte;

HAbe eine POINT Array gemacht .. welche amplituten einer sinukurve (eine Periode) in 400 Schritten enthält und diese dann in deinen verte vector übernomen...
Funktioniertr das dann so???

[Phoemuex]

Im Prinzip müsste es so gehen. Ich sehe nur zwei mögliche "Probleme":

Du schreibst

d+=(M_PI*2)/NUMBER

, d.h. du erhöhst d immer soweit, dass alle Änderungen zusammen 2Pi ergeben. Das ist auch normal so. Aber dann schreibst du

pWerte[i].x=i+1;

. Das heißt, du erhöhst den X-Wert des Punktes immer um 1 und nicht um 2Pi/NUMBER. Das kann aber auch gewollt sein.

Du musst dann zustäzlich noch in der Zeile

vector<pair<double,double> > koeffizienten = FourierCoefficients(anfangsIterator, endIterator, 1.0/File.GetSamplingRate(), 400);

Das 1.0/File.GetSamplingRate() durch 1 oder durch

(M_PI*2)/NUMBER

ersetzen.

Ansonsten müsste es gehen, wobei ich nicht weiß, ob in einer Variable vom Typ POINT double-, oder Int-Werte gespeichert werden...

Felix

EDIT: Evtl. musst du die FUnktionen noch so ändern, dass sie vector<double> anstatt vector<int> benutzen. So wie du es im Moment machst (mit r = 50) sind die Werte groß genug um auch mit int zu funktionieren, wenn r aber zum Beispiel = 1 ist, sind alle Werte -1<x<1 gleich 0 und dann funktioniert das mit den Funktionen nicht richtig. Ist aber relativ einfach, dass zu ändern.

[BorisProgger]

in POINT werden sie als INT abgespeichert!!

POINT Array (x) ist quasie die x achse von 0- 400, und die jeweilige Amplitude

Das Rrogramm läuft durch und ich habe auch 400 Koeffizienten das passt ja!!

Wenn ich nun alles 400 Koofizenten zusammen rechne, müsst doch wieder die Original Sinus kurve rauskommen oder nich???

/////////////////// KOEFFIZENTEN SUMMIEREN ////////////////////	
	double res=0.0;
	double *pResult = new double[NUMBER];
	for(int i=0; i<NUMBER;i++){
		//pResult[i].x=i;
		t=0;
		for(vector<pair<double,double> >::const_iterator it = koeffizienten.begin(); it != koeffizienten.end(); ++it,t++){
			if(t%2==0)
				res+=it->first*cos(it->second*i*((2*M_PI)/NUMBER));
			else 
				res+=it->first*sin(it->second*i*((2*M_PI)/NUMBER));
		}
		pResult[i]=res;
		res=0.0;
	}

[Flow_cplus]

wenn, um arraydaten von sins oder geht, mach doch einfach ne FFT, und schon hast frequenzen und phasen!

[BorisProgger]

der code oben ist ne Fourieranalyse.. oder gibts da schon irgendweche (verständliche/übersichtliche) FFT alogirthem!!
ich will ein Array von Kurven-werten analsysieren, und dann auch ein array mit den gefundenen Schwingungnen (amplitude und freq) zurückahaben aber sowas gibts nich.. bzw. hab ich nix gefunden