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...

[3_of_8]

*push*

[dino]

da du ja schon selber am pushen bist, kann ich vielleicht mal sagen, dass ich mich mal mit fft beschäftigt habe und es auch mal programmieren wollte, aber schnell demotiviert wurde

und tatsächlich: es ist kompliziert(wenn ich dein code mal so sehe)

aber ich möchte auch, dass jemand, der mehr erfahrung hat dir sagt, was du wissen willst, also: wenn jemand die antwort kennt, so soll er dies alles hier ignorieren

[sirius]

Zitat von 3_of_8:

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...

Ich habe bisher nur einfache DFT-Algos auf eine spezielle Anzahl von Abtastwerten bzw. Samplingraten "optimiert"/angepasst. Deswegen kann ich nicht genau sagen, was an deiner FFT flasch, oder fehlerhaft ist. Was komisch aussieht, hast du ja selber bemängelt:
Das DEGtoRad muss falsch sein.
Und was ist dieses w^7, was du am Anfang der Funktion machst.

Und nebenbei: Ist dir klar, was du am Ende einer DFT/FFT für Ergebnisse bekommst? Sagt dir Aliasing/Antialiasing etwas?


PS:
Warum denn gleich FFT? Reicht nicht erstmal eine einfache DFT? Da versteht man noch eher was passiert.


Edit: Wenn ich von der DFT ausgehe, müsstest du für w am Anfang die konjugiert Komplexe nehmen.
du tust die dir anscheinend mit deinem w^7 irgendwie zurechtzubasteln. Das dürfte aber nur bei einer bestimmten Zahl an Abtastwerten klappen.

[IngoD7]

Zitat von sirius:

Zitat von 3_of_8:

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...

Was komisch aussieht, hast du ja selber bemängelt:
Das DEGtoRad muss falsch sein.

Ich habe keine Ahnung, was FFT und dergleichen ist, aaaaber:

Die Winkelangaben dieser Funktionen müssen in Delphi im Bogenmaß angegeben werden. Daher muss in Delphi eine Winkelangabe von z.B. 30° zuvor mit DegToRad umgewandelt werden, damit der Sinus von 30° auch wirklich 0,5 ergibt.

Wenn also in cos(DegToRad(2*Pi/length(a))) das 2*Pi/length(a) schon Bogenmaß ist, dann ist DegToRad falsch und deine Ergebnisse wohl nur zufällig "einigermaßen passend". Wenn das 2*Pi/length(a) eine Grad-Angabe ist, dass ist das DegToRad korrekt, um in Delphi auf den korrekten Cos-Wert zu kommen.


Kann aber natürlich auch sein, dass ich euch überhaupt nicht verstanden habe ...

[EWeiss]

Einen eigenen FFT Algo zu programmieren setzt ein fundiertes wissen an Mathematik vorraus.
Wenn du nicht darüber verfügst dann lass es lieber und benutze vorhandene algos.

Irgendwo macht das auch keinen sinn etwas zu erfinden was schon vorhanden ist.
Es sein denn !!!
Du könntest ihn optimieren, damit haben sich aber schon einige Professoren beschäftigt
und viele sind gescheitert.

Information darüber was FFT (fast Fourier transform) bedeutet bitte hier ein link.

http://de.wikipedia.org/wiki/Schnell...Transformation

gruß

[sirius]

Zitat von EWeiss:

Information darüber was FFT (fast Fourier transform) bedeutet bitte hier ein link.
http://de.wikipedia.org/wiki/Schnell...Transformation
gruß

Ich glaube, den Link und den weiterfürhenden Link hat er schon selber gefunden.

Ich bin zwar der Meinung, dass man durch Abschreiben (und gleichzeitiges darüber Nachdenken) auch ein Menge lernen kann. Aber bei so komplexen Sachen...
Hut ab, vor dem, der es schafft.

[3_of_8]

Mein Wissen über Mathematik ist immerhin etwa auf Erstsemesterniveau, die Gleichungen der FFT verstehe ich größtenteils.

Zitat:

Und was ist dieses w^7, was du am Anfang der Funktion machst.

Das hätte eigentlich schon längst draußen sein sollen, das habe ich nur gemacht, um zu testen, ob w^n=1.

@Dino: Du wurdest nicht demotiviert, dir wurde nur gesagt, dass eine FFT garantiert nicht das tut, was du tun willst, weil das, was du tun wolltest, nicht möglich ist.

[sirius]

Was gibst du den zum Testen für eine Zeitreihe ein, und was für Ergebnisse bekommst du?

[3_of_8]

Ich weiß nicht so genau, was ich erwarte.

Deswegen probier ichs einfach mal aus mit einem array[0..127] of Single, die die Werte einer Sinuskurve enthalten.

EDIT: Ich schätze mal, das ganze würde funktionieren, wenn ich wenigstens w ausrechnen könnte.

Also, ich habe w^n=1

w=e^(2*Pi*i/n)

1. Eulersche Formel:
e^(Phi*i)=cos(Phi)+sin(Phi)*i

Also Umformung:
w=cos(2*Pi/n)+sin(2*Pi/n)

Nur erhalte ich dann Ergebnisse, bei denen wenn w die 7. primitive Einheitswurzel ist, Re(w^7) bei ~1 liegt, aber Im(w^7) bei 0.1.