Numerische Faltung zweier Funktion



  • Hallo,

    ich würde gerne zwei Funktionen numerisch falten.
    Die Funktionen sind an n Punkten ausgewertet und in den Arrays f und g gespeichert. Beide Funktionen fallen zum Rand hin schnell auf Null ab, d.h. die Funktionen sind nur an sehr wenigen Stützstellen in der Mitte von Null verschieden.
    Für die Faltung muss ich das Integral über x' von minus bis plus unendlich über
    f(x') * g(x-x') ausrechnen.

    Leider komme ich beim programmieren irgendwie immer mit den Indizes durcheinander, mein Ansatz war:

    for(int i=0; i<n; i++)
    {
    	for(int j=0; j<n; j++)
    	{	
    		if( (i-j) > 0 )
    		{
    		      ergebnis[i] += f[j]*g[i-j]*weight[j];
    		}
    	}
    }
    

    In dem array weight sind die Gewichte für die Integration gespeichert (die Gewichte stimmen auch definitiv).

    Leider kommt bei mit immer Mist raus. Ich vermute, dass das ganze irgendwas damit zu tun hatm, dass ich nur für j=0..i aufsummiere.
    Habt ihr eine Idee wo mein Fehler liegt?



  • Ich weiss ehrlich gesagt nicht was numerisches Falten ist (bin erst im 1. Semester ^^) aber ich habe deinen Post trotzdem durchgelesen und mir ist was aufgefallen.

    Du schreibst:

    dass ich nur für j=0..i aufsummiere.

    In deinem Code summierst du allerdings für j = 0,...,n-1 auf.

    Keine Ahnug ob dir was hilft. Vielleicht hast du dich ja auch nur vertippt und wolltest n schreiben.



  • @icarus2: Da ist ein "if( (i-j) > 0 )". Warum man allerdings nicht "if(i > j)" oder gleich "j<=i" als Abbruchkriterium nimmt, ist mir allerdings auch schleierhaft.

    Der Fehler kann vieles sein. Stimmt n? Ist 'ergebnis' mit 0 initialisiert? usw.



  • Ich würde auch darauf Tippen, dass die Variable ergebnis nicht mit 0 initialisiert wurde.

    Zudem ist das irgendwie nicht alles ganz eindeutig mit deinen Indizes. Die Formel für die numerische Faltung schaut so aus:
    $y(kT) = i=0kx(iT)h(kTiT)\displaystyle\sum\limits_{i=0}^k x(iT)*h(kT-iT)$

    Also wieso nicht einfach C++ daraus machen:

    int y[N] = {0}; // Alles wird mit 0 initalisiert
    size_t K = ... // Faltung bis zu welchem Element durchführen
    
    for(size_t k = 0; k <= K; ++k)
    {
        for(size_t i = 0; i <= k; ++i)
            y[k] += x[i] * h[k-i];
    }
    

Anmelden zum Antworten