1D und 2D FFT DFT rekursiv

[alleinherrscher]

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]

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]

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]

Korrekt :) Ich werde es oben einbauen.

[Michael II]

Bei Delphi 11 wird zur Laufzeit ein "Variant Fehler" ausgegeben.
Lösung: In FFT und in iFFT v als v:variant; definieren.

Auch wenn's ein paar Jahre her ist: Herzlichen Dank für den Code. Ich nutze diesen für Polynominterpolation in der Knotentheorie.