Problem bei FFT

[3_of_8]

Morgen.

Ich versuche gerade, eine FFT zu implementieren, aber irgendwie kriege ich das nicht so ganz hin...

Delphi-Quellcode:

unit complex;

interface

uses Math;

type
  TComplex=record
    re, im: Single; //Real- und Imaginärteil
  end;

  TComplexArray=array of TComplex;

function AddC(a, b: TComplex): TComplex; overload;
function AddC(a: TComplex; b: Single): TComplex; overload;
function SubC(a, b: TComplex): TComplex; overload;
function SubC(a: TComplex; b: Single): TComplex; overload;
function MulC(a, b: TComplex): TComplex; overload;
function MulC(a: TComplex; b: Single): TComplex; overload;

function MakeC(re, im: Single): TComplex;

procedure FFT(var a: TComplexArray); overload;
procedure FFT(var a: array of Single); overload;
procedure FFT(var a: array of Integer); overload;

implementation

function AddC(a, b: TComplex): TComplex;
begin
  Result.re:=a.re+b.re;
  Result.im:=a.im+b.im;
end;

function AddC(a: TComplex; b: Single): TComplex;
begin
  Result.re:=a.re+b;
  Result.im:=a.im;
end;

function SubC(a, b: TComplex): TComplex;
begin
  Result.re:=a.re-b.re;
  Result.im:=a.im-b.im;
end;

function SubC(a: TComplex; b: Single): TComplex;
begin
  Result.re:=a.re-b;
  Result.im:=a.im;
end;

function MulC(a, b: TComplex): TComplex;
begin
  Result.re:=a.re*b.re-a.im*b.im;
  Result.im:=a.re*b.im+a.im*b.re;
end;

function MulC(a: TComplex; b: Single): TComplex;
begin
  Result.re:=a.re*b;
  Result.im:=a.im*b;
end;

function MakeC(re, im: Single): TComplex;
begin
  Result.re:=re;
  Result.im:=im;
end;

procedure shuffle(var a: TComplexArray; n, lo: Integer); //Mischt die Elemente des Arrays
var I, m: Integer;
    b: TComplexArray;
begin
  m:=n shr 1;
  setlength(b, m);
  for I:=0 to m-1 do
    b[i]:=a[lo+i];
  for I:=0 to m-1 do
    a[lo+i+i+1]:=a[lo+i+m];
  for I:=0 to m-1 do
    a[lo+i+i]:=b[i];
end;

procedure DoFFT(var a: TComplexArray; n, lo: Integer; w: TComplex); //Führt die FFT rekursiv aus
var I, m: Integer;
    z, v, h: TComplex;
begin
    //Hier ist w möglicherweise falsch
    z:=MulC(w, MulC(w, MulC(w, MulC(w, MulC(w, MulC(w, w))))));
    if n>1 then
    begin
        m:=n shr 1; //Teilt die Arraylänge auf
        z:=MakeC(1, 0); //Soviel wie z:=1, nur eben als komplexe Zahl
        for I:=lo to lo+m-1 do
        begin
            h:=SubC(a[i], a[i+m]);
            a[i]:=AddC(a[i],a[i+m]);
            a[i+m]:=MulC(h,z);
            z:=MulC(z,w);
        end;
        v:=MulC(w,w);
        DoFFT(a, m, lo, v);
        DoFFT(a, m, lo+m, v);
        shuffle(a, n, lo);
    end;
end;

//In den folgenden Prozeduren ist exp(2*Pi/length(a)) jeweils e^2*Pi*n, 
//die primitive Einheitswurzel

procedure FFT(var a: TComplexArray);
begin
  DoFFT(a, length(a), 0, MakeC(cos(DegToRad(2*Pi/length(a))), sin(DegToRad(2*Pi/length(a)))));
end;

procedure FFT(var a: array of Single);
var I: Integer;
    b: TComplexArray;
begin
  setlength(b, length(a));
  for I:=0 to high(a) do b[I]:=MakeC(a[I], 0);
  DoFFT(b, length(a), 0, MakeC(cos(DegToRad(2*Pi/length(a))), sin(DegToRad(2*Pi/length(a)))));
  for I:=0 to high(a) do a[I]:=b[I].re;
end;

procedure FFT(var a: array of Integer);
var I: Integer;
    b: TComplexArray;
begin
  setlength(b, length(a));
  for I:=0 to high(a) do b[I]:=MakeC(a[I], 0);
  DoFFT(b, length(a), 0, MakeC(cos(DegToRad(2*Pi/length(a))), sin(DegToRad(2*Pi/length(a)))));
  for I:=0 to high(a) do a[I]:=trunc(b[I].re);
end;

end.

EDIT: Soo, das sieht jetzt schon *etwas* besser aus.
Der Algorithmus kompiliert, läuft und terminiert ohne Exceptions und die Ergebnisse für w sind auch schon fast in Ordnung.

Wo ich mir aber nicht ganz sicher bin, das ist die Berechnung von w.
Wenn ich das richtig sehe, dann ist w=e^(2*Pi*i/n) und e^(i*x)=cos(x)+sin(x)*i. Dadurch lässt sich das ganze auflösen in w=cos(2*Pi/n)+sin(2*Pi/n)*i, was ich oben auch habe. Nur hätte ich eigentlich angenommen, dass mit sin und cos die Sinus- und Kosinusfunktionen im Bogenmaß gemeint sind, aber einigermaßen passende Ergebnisse bekomme ich nur mit Gradangaben, also vorher ein DegToRad...