Delphi-PRAXiS

Delphi-PRAXiS (https://www.delphipraxis.net/forum.php)
-   Neuen Beitrag zur Code-Library hinzufügen (https://www.delphipraxis.net/33-neuen-beitrag-zur-code-library-hinzufuegen/)
-   -   Delphi 1D und 2D FFT DFT rekursiv (https://www.delphipraxis.net/175009-1d-und-2d-fft-dft-rekursiv.html)

alleinherrscher 24. Mai 2013 12:13

1D und 2D FFT DFT rekursiv
 
Im Gegensatz zu den hier im Forum bereits bekannten (diskreten) FFT Algorithmen, präsentiere ich hier einen rekursiven Algorithmus mit der Delphi nativen Unterstützung für komplexe Zahlen. Der Code wird dadurch sehr kurz und übersichtlich. Der Pseudocode kann dem Wikipedia Artikel zum Thema FFT entnommen werden.

Hinweis: Die Anzahl an Komponenten (n) des Eingangsvektors f mit einer Zahl 2^n entsprechen.

Delphi-Quellcode:
uses System.VarCmplx;
Delphi-Quellcode:
type
  TComplexVector = array of Variant;
  TFourierMatrix = array of array of Variant;
Delphi-Quellcode:
function FFT(n:integer;f:TComplexVector):TComplexVector;
var g,u,v,f1,f2:TComplexVector;
    k:integer;
begin
 if n=1 then
    result:=f
 else
   begin
     setlength(f1, n div 2);
     setlength(f2, n div 2);
     for k := 0 to (n div 2)-1 do
       begin
         f1[k]:=f[2*k+1];
         f2[k]:=f[2*k];
       end;

     g:=FFT(n div 2,f2);
     u:=FFT(n div 2,f1);

     setlength(result,n);

     for k := 0 to (n div 2)-1 do
       begin
         v := VarComplexExp(VarComplexCreate(0,-2*Pi*k/n));
         result[k]:=g[k]+u[k]*v;
         result[k+(n div 2)]:= g[k] - u[k]*v;
       end;
   end;
end;
Für die inverse FFT (iFFT) sind lediglich die Minuszeichen in der Exponentialfunktion zu entfernen:

Delphi-Quellcode:
function iFFT(n:integer;f:TComplexVector):TComplexVector;
var g,u,v,f1,f2:TComplexVector;
    k:integer;
begin
 if n=1 then
    result:=f
 else
   begin
     setlength(f1, n div 2);
     setlength(f2, n div 2);
     for k := 0 to (n div 2)-1 do
       begin
         f1[k]:=f[2*k+1];
         f2[k]:=f[2*k];
       end;

     g:=FFT(n div 2,f2);
     u:=FFT(n div 2,f1);

     setlength(result,n);

     for k := 0 to (n div 2)-1 do
       begin
         v := VarComplexExp(VarComplexCreate(0,2*Pi*k/n));        
         result[k]:=g[k]+u[k]*v;
         result[k+(n div 2)]:= g[k] - u[k]*v;
       end;
   end;
end;
Beide Funktionen können auch zweidimensional verwendet werden. Hierfür werden erst alle Zeilen einer Matrix und dann alle Spalten Fouriertransformiert:

Hinweis: Auch hier gilt, dass die Dimension der Matrix (w,h) in beiden Komponenten w=2^n bzw h=2^m geschrieben werden kann.

Delphi-Quellcode:
function FFT2D(Input:TFourierMatrix;w,h:integer):TFourierMatrix;
var F,R:TComplexVector;
  j: Integer;  
  i: Integer;  
begin

setlength(F,w);
setlength(result,w,h);

for j := 0 to h-1 do
   begin
      for i := 0 to w-1 do
          F[i]:=Input[i,j];
      R:=FFT(w,F);
      for i := 0 to w-1 do
        result[i,j]:=R[i];
   end;

setlength(F,h);
for i := 0 to w-1 do
   begin
      for j := 0 to h-1 do
          F[j]:=result[i,j];
      R:=FFT(h,F);
      for j := 0 to h-1 do
        result[i,j]:=R[j];
   end;

end;
iFFT2D folgt analog durch ersetzen der procedure FFT durch iFFT in die procedure FFT2D.

sx2008 24. Mai 2013 13:29

AW: Vorschlag für CodeLib: 1D und 2D FFT DFT rekursiv
 
Kann das sein, dass 2 mal der gleiche Ausdruck
Delphi-Quellcode:
VarComplexExp(VarComplexCreate(0,-2*Pi*k/n))
berechnet wird?
Den Faktor bräuchte man doch nur einmal in der Schleife zu berechnen:
Delphi-Quellcode:
for k := 0 to (n div 2)-1 do
begin
   v := VarComplexExp(VarComplexCreate(0,-2*Pi*k/n));
   result[k]:=g[k]+u[k]*v;
   result[k+(n div 2)]:= g[k] - u[k]*v;
end;

DP-Maintenance 24. Mai 2013 20:17

Dieses Thema wurde am "24. May 2013, 20:17 Uhr" von "TBx" aus dem Forum "Algorithmen, Datenstrukturen und Klassendesign" in das Forum "Neuen Beitrag zur Code-Library hinzufügen" verschoben.

alleinherrscher 25. Mai 2013 16:48

AW: 1D und 2D FFT DFT rekursiv
 
Korrekt :) Ich werde es oben einbauen.


Alle Zeitangaben in WEZ +2. Es ist jetzt 03:53 Uhr.

Powered by vBulletin® Copyright ©2000 - 2021, Jelsoft Enterprises Ltd.
LinkBacks Enabled by vBSEO © 2011, Crawlability, Inc.
Delphi-PRAXiS (c) 2002 - 2021 by Daniel R. Wolf