Fourier Transform DFT, FFT

[Dipl Phys Ernst Winter]

Diskrete Fourier Transformation DFT
Eine periodische Funktion der normierten Periode 2Pi sei mit n = 2q äquidistanten Stützstellen f(xi) i = 0..n-1 im Intervall 0..2Pi gegeben. Stützstellenabstand dx = 2Pi/n, x0 = 0, xn-1 = 2Pi - dx.

Delphi-Quellcode:

Wir interpolieren sie mit:
            q - 1 
f(x) = a0 + Summe (ai*cos(mx) + bi*sin(mx)) + aq*cos(qx) 
             m =1
wobei wir b(q) = 0 setzen, um die Zahl der Parameter an die Zahl der Stützstellen anzupassen.
Brechen wir die Summe mit weniger Gliedern ab, so erhalten wir eine Approximation mit minimalem quadratischen Fehler. Für die n + 1 Koeffizienten ergibt sich

         n - 1 
a0 = 1/n Summe yk
         k = 0

         n - 1 
ai = 2/n Summe yk*cos(i*xk)         i = 1, 2,..., q-1
         k = 0

         n - 1 
bi = 2/n Summe yk*sin(i*xk)         i = 1, 2,..., q-1
         k = 0

         n - 1 
aq = 1/n Summe yk*cos(q*xk)         
         k = 0

FFT Fast Fourier Transform
Die Rechenzeiten der DFT wachsen mit Stützstellenzahl N quadratisch an: t ~ N^2. Es wurden verschiedene Verfahren zur schnellen Fourier-Transformation FFT entwickelt, deren die Rechenzeit nur mit t ~ Ln(N)*N anwächst.
Sie beruhen alle auf der sukzessiven Zerlegung einer Transformation mit n Stützstellen in zwei Transformationen mit n/2 Stützstellen.

Das Demo zur DFT zeigt den Aufbau der Näherungen aus den Oberwellen graphisch an.

FFT-Demo1 demonstriert den Zeitgewinn der FFT gegenüber der DFT. Vorsicht: Die DFT kann bei großem nMax sehr lange dauern. Es enthält eine Implementation der FFT nach dem Tukey-Algorithmus.

Das Demo2 zur FFT enthält eine Implementation der FFT die reelle Werte wahlweise in der Frequenzbereich bzw. von dort wieder zurück in die Zeitdarstellung transformatiert.

[EWeiss]

Ist zu hoch für mich ;)
Trotzdem ein :thumb: für deine Arbeit

OT:
Möchte damit meine Hochachtung ausdrücken.

gruss Emil

[jfheins]

Sieht sehr gut aus - insbesondere die erste Demo mit der grafischen Darstellung gefällt mir :)

[stoxx]

ich begrüße es und finde es sehr schön, dass Du Dich etwas um die wissenschaftliche Seite kümmerst :-)
Da fehlen oft Quelltexte in Pascal, wenn man danach sucht.

Was mich und vielleicht auch andere interessieren würde, wäre eine Wavelet Transformation.

Hier gibts zwar schon Pascal Quelltexte

http://www.basegroup.ru/download/fre...wavelet_utils/

aber da ist noch ein Fehler drin irgendwie. Die Approximation der letzten Datenpunkte einer Zeitreihe funktioniert dort überhaupt nicht.

Was dann zu groben Fehlern bei der Rückapproximation ergibt.

Im Anhang sieht man das mal. die Rote Wavelet Kurve ist die Approximation einer Wavelet Transformation auf Stufe 10
Dummerweise geht die rote Kurve im letzten Bereich nach oben, obwohl die grüne Originalkurve eher nach unten tendiert.

Habe mich zeitlich nicht so recht in den Quelltext einarbeiten können, ein einfaches Beispiel für Wavelets wäre ziemlich gut :-)
z.b. eines Daubechies-Wavelets ( nicht unbedingt für Haar- Wavelets)

[Dipl Phys Ernst Winter]

"stoxx"

Zitat:

Da fehlen oft Quelltexte in Pascal, wenn man danach sucht.

Hinweis: G. Engeln-MNüllges, F. Reuter: Formelsammlung zur Numerischen Mathemetik mit Turbo PASCAL Programmen.
Verlag und Erscheinungsjahr habe ich nicht zu Hand. Müsste sich aber auch ohne finden lassen.
PS: Der Pascal-Stil ist jedoch grauenhaft.

Mit Wavelets habe ich mich noch nicht beschäftigt.

[DP-Maintenance]

Dieses Thema wurde von "Daniel G" von "Neuen Beitrag zur Code-Library hinzufügen" nach "Tutorials und Kurse" verschoben.
Aufgrund des Umfangs eher für die Tutorial-Sparte geeignet als für die Code-Lib.

[omata]

Hier mal eine überarbeitete Version, ohne globale Variablen und mit definierten Schnittstellen zwischen den einzelnen Units.
Die Bereichsfehlerprüfung darf jetzt auch aktiviert werden und es sind keine Speicherlecks mehr vorhanden.

[R2009]

Hi Ernst,

ich habe deine FFT ausprobiert und bekomme eine Fehlermeldung an dieser Stelle:

Delphi-Quellcode:

    IF j<=sigma THEN BEGIN
      j_2:= j SHL 1; sigma2:= sigma SHL 1;
      u.Re:= Y[j_2]; u.Im:= Y[j_2+1];
      Y[j_2]:= Y[sigma2]*faktor;           <<<<<<<< ungültige Gleitkommaoperation
      Y[j_2+1]:= Y[sigma2+1]*faktor;
      Y[sigma2]:= u.Re*faktor; Y[sigma2+1]:= u.Im*faktor END END;

Aufgerufen habe ich das Ganze:

Delphi-Quellcode:

procedure TForm1.Button1Click(Sender: TObject);
var s:array[0..500] of extended;
begin
  s[0]:=0;s[1]:=1;s[2]:=2;s[3]:=3;s[4]:=4;s[5]:=5;
  real_fft(6,s,false);
end;

Wo liegt mein Denkfehler? Woher kommt die Fehlermeldung? Habe beim debuggen nichts unkeusches gefunden.


OK habs schon gefunden. s muss initialisiert werden.

Delphi-Quellcode:

  for n:=0 to 500 do s[n]:=0;

Grüsse
Rainer

[R2009]

Hi Ernst,

die Rückgabewerte aus dieser Berechnung sehen etwas diffus aus.
Bitte erkläre uns was diese Werte physikalisch bedeuten.
Bei meiner Implementierung entsprechen die Werte den Amplituden der Frequenzen im Zeitbereich.

Mit dem oben genannten Beispiel habe ich folgende Ergebnisse bekommen:

0,109375
0,216908595604189
0,0267547521355106
0,211425330450812
0,0529725984625451
0,202422030471313
0,0781293830536248
0,190098371826692
0,101726122007646
0,174726923542575
0,123300778934621
0,156646391024367
0,142439092232398
0,136253238243791
0,15878417807239
0,113991908753758
0,172044665859018
0,0903439020560376
0,182001162251785
0,065815991232699

Grüsse
rainer

[sirius]

Zitat:

Zitat von R2009

Hi Ernst,

die Rückgabewerte aus dieser Berechnung sehen etwas diffus aus.
Bitte erkläre uns was diese Werte physikalisch bedeuten.
Bei meiner Implementierung entsprechen die Werte den Amplituden der Frequenzen im Zeitbereich.

Ich würde sagen, er macht eine DFT aus deiner Zeitfunktion mit den Werten:
[0; 1 ; 2; 3; 4; 5; 0; 0; 0; 0; 0; ...; 0];
Was erwartest du denn da für ein Spektrum?


Edit: OK, in Matlab kommt auch ein anderes Spektrum heraus.