AGB  ·  Datenschutz  ·  Impressum  







Anmelden
Nützliche Links
Registrieren
Zurück Delphi-PRAXiS Tutorials Delphi Fourier Transform DFT, FFT

Fourier Transform DFT, FFT

Ein Tutorial von Dipl Phys Ernst Winter · begonnen am 16. Mai 2009 · letzter Beitrag vom 20. Okt 2009
Antwort Antwort
Seite 1 von 2  1 2   
Dipl Phys Ernst Winter
Registriert seit: 14. Apr 2009
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.
Angehängte Dateien
Dateityp: zip fft_demo2_143.zip (149,8 KB, 320x aufgerufen)
Dateityp: zip fft_demo1_393.zip (160,2 KB, 230x aufgerufen)
Dateityp: zip dft_demo_982.zip (167,9 KB, 230x aufgerufen)
Autor: DP Ernst Winter
 
EWeiss
 
#2
  Alt 16. Mai 2009, 17:43
Ist zu hoch für mich
Trotzdem ein für deine Arbeit

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

gruss Emil
  Mit Zitat antworten Zitat
Benutzerbild von jfheins
jfheins
 
#3
  Alt 16. Mai 2009, 18:05
Sieht sehr gut aus - insbesondere die erste Demo mit der grafischen Darstellung gefällt mir
  Mit Zitat antworten Zitat
Benutzerbild von stoxx
stoxx
 
#4
  Alt 16. Mai 2009, 18:09
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)
Miniaturansicht angehängter Grafiken
r_ckapproximation_919.png   wavelet_details_level_5_210.png  
  Mit Zitat antworten Zitat
Dipl Phys Ernst Winter

 
Delphi 3 Professional
 
#5
  Alt 16. Mai 2009, 20:44
"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.
  Mit Zitat antworten Zitat
19. Okt 2009, 08:52
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

 
Delphi 7 Enterprise
 
#7
  Alt 20. Okt 2009, 02:05
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.
Angehängte Dateien
Dateityp: zip dft_demo_134.zip (19,5 KB, 97x aufgerufen)
  Mit Zitat antworten Zitat
R2009

 
Delphi 2007 Professional
 
#8
  Alt 20. Okt 2009, 07:06
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.
  for n:=0 to 500 do s[n]:=0;
Grüsse
Rainer
Rainer Unger
  Mit Zitat antworten Zitat
R2009

 
Delphi 2007 Professional
 
#9
  Alt 20. Okt 2009, 07:27
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
Rainer Unger
  Mit Zitat antworten Zitat
Benutzerbild von sirius
sirius

 
Delphi 7 Enterprise
 
#10
  Alt 20. Okt 2009, 09:47
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.
  Mit Zitat antworten Zitat
Themen-Optionen Tutorial durchsuchen
Tutorial durchsuchen:

Erweiterte Suche
Ansicht

Forumregeln

Es ist dir nicht erlaubt, neue Themen zu verfassen.
Es ist dir nicht erlaubt, auf Beiträge zu antworten.
Es ist dir nicht erlaubt, Anhänge hochzuladen.
Es ist dir nicht erlaubt, deine Beiträge zu bearbeiten.

BB-Code ist an.
Smileys sind an.
[IMG] Code ist an.
HTML-Code ist aus.
Trackbacks are an
Pingbacks are an
Refbacks are aus

Gehe zu:

Impressum · AGB · Datenschutz · Nach oben
Alle Zeitangaben in WEZ +1. Es ist jetzt 10:42 Uhr.
Powered by vBulletin® Copyright ©2000 - 2024, Jelsoft Enterprises Ltd.
LinkBacks Enabled by vBSEO © 2011, Crawlability, Inc.
Delphi-PRAXiS (c) 2002 - 2023 by Daniel R. Wolf, 2024 by Thomas Breitkreuz