Seltsame Ergebnisse bei DFT

[inherited]

Hi,
Ich, bzw. hauptsächlich der Borg hat versucht, ein DFT zu implementieren um aus einer wav-datei die Frequenzen zu bekommen. Doch offenbar mag das alles nicht so wie es soll. Manchmal kommen ganz authentische Ergebnisse raus und dann wieder nicht. es kommen Ergebnisse raus die nicht rauskommen können.
Sie liefert richtige Ergebnisse bei zB 440 Herz, bei 300 Herz funktioniert es bei zB 1000 Abtastwerten, aber bei zB 2000 nicht, da gibt es seltsame Ergebnisse (42 wär ja logisch, aber nicht mal das kommt raus :stupid:) (Anzahl der Abtastwerte=Länge des Arrays)
Anbei die DFT-Funktion.
Es wäre nett wenn sich das mal jemand anschauen könnte und uns zu einem *In-die-Holzkannte-Beiss-natürlich-da-liegt-das-Problem*-Erlebnis verhilft^^
Danke.

Delphi-Quellcode:

procedure DoDFT(var a: TComplexArray; n, lo: Integer; w: TComplex);  
var b, f: array of TComplex;  
    I, J: Integer;  
begin  
  setlength(b, n);  
  for I:=0 to n-1 do  
  begin  
    b[I]:=a[I];  
    a[I]:=MakeC(0, 0);  
  end;  
  setlength(f, n);  
  f[0]:=MakeC(1, 0);  
  f[1]:=w;  
  for I:=2 to n-1 do  
    f[I]:=MulC(f[I-1], w);  
  for I:=0 to n-1 do  
    for J:=0 to n-1 do  
    begin  
      a[I]:=AddC(a[I], MulC(b[J], f[I*J mod n]));  
    end;  
end;

[3_of_8]

Nachtrag:

AddC und MulC sind die Addition bzw. Multiplikation von komplexen Zahlen, MakeC erstellt eine komplexe Zahl aus den angegebenen Werten. Diese Funktionen sind garantiert richtig implementiert.

a ist der "Eingabevektor" (a1,...,an), der verändert wird.
n ist die Anzahl der Elemente in a.
lo ist hier nicht von Bedeutung, ich werds aus der Parameterliste rausstreichen.
w ist die n-te primitive Einheitswurzel.

f[I] ist das Array, das den Wert von w^I speichert.
b ist der Vektor, der die "ursprünglichen" Werte von a speichert, die bei der Matrixmultiplikation benötigt werden, weil die Werte in a ja verändert werden und die Matrixmultiplikation sonst falsch wäre.

EDIT: Nebenbei bemerkt ist dieser Algorithmus trotz einiger Optimierungen, die ich mir ausgedacht habe, sehr, sehr langsam. Nicht verwunderlich, wenn man bedenkt, dass bei n=2000 einige Millionen (Milliarden?) Additionen und Multiplikationen ausgeführt werden müssen.

[sirius]

Funktioniert doch :gruebel:

Ich hab alles so verwendet wie bei deiner FFT, mit folgenden Änderungen:

Delphi-Quellcode:

rocedure TForm1.Button1Click(Sender: TObject);
const Abtastwerte=3000;
      Abtastfrequenz=1000;
var l:Tlineseries;
    a:array[0..Abtastwerte-1] of single;
    i:integer;
begin
  for i:=0 to Abtastwerte-1 do
    a[i]:=2             //Gleichanteil
          +5*sin(2*pi*50*i/abtastfrequenz)  //50Hz
          +3*sin(2*pi*440*i/Abtastfrequenz) //440Hz
          +1*sin(2*pi*300*i/Abtastfrequenz);//300Hz
  fft(a);

  chart1.FreeAllSeries;
  l:=tlineseries.create(chart1);
  chart1.AddSeries(l);
  for i:=0 to Abtastwerte-1 do l.AddXY(i*abtastfrequenz/Abtastwerte,a[i]);
end;

Delphi-Quellcode:

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);
  DoDFT(b, length(a), 0, MakeC(cos((2*Pi/length(a))), sin((2*Pi/length(a)))));
  for I:=0 to high(a) do a[I]:=sqrt(sqr(b[I].re)+sqr(b[i].im))/length(a);  //nur den Betrag
end;

+Exakt der angegebenen DoDFT
-->Bild ist im Anhang

Zitat:

EDIT: Nebenbei bemerkt ist dieser Algorithmus trotz einiger Optimierungen, die ich mir ausgedacht habe, sehr, sehr langsam. Nicht verwunderlich, wenn man bedenkt, dass bei n=2000 einige Millionen (Milliarden?) Additionen und Multiplikationen ausgeführt werden müssen.

Na ein "sehr" kann man weglassen (zumindest bei 2GHz Dual Core) :wink:
Aber nicht umsonst hat die FFT ihren Namen bekommen. Sie hat aber eben den Nachteil, dass du immer 2^n Abtastwerte nehmen musst. Wenn das nicht geht, musst du eben auf die DFT zurückgreifen.
Nebenbei bemerkt: Man kann die For-Schleifen noch optimieren. Dazu müsste man wahrscheinlich zu Assembler gehen. Mathworks schaffts nämlich auch ohne merkliche Zeitverzögerung (eben ausprobiert). [Edit: Quatsch Mathworks, die benutzen auch nur das hier]

Du kannst mir ja mal sagen, was genau du vor hast. Wahrscheinlich kann man es für eine FFT hinbiegen.

[3_of_8]

Das wars!

Als du von Betrag gesprochen hast, dachte ich an Abs(Re(b[I])), allerdings war es wohl eher sqrt(Re(b[I])²+Im(b[I])²). Meine DFT war also korrekt, nur die Auswertung war falsch.

Ich geb dir nochmal einen an der DP-Bar aus. ;)

Ich habe schon überlegt, wie ich eine FFT machen kann. Die Anzahl der Abtastwerte wird exakt vorgegeben, ich hab überlegt, ob ich nicht die vorhandenen Werte auf die nächste 2-er-Potenz extrapolieren kann...

[sirius]

Zitat:

Als du von Betrag gesprochen hast

Das ist das, was ich meinte. Du kannst sicherlich hervorragend mit komplexen Zahlen umgehen. Dir fehlt aber einfach noch das Gefühl für diesen Zahlenbereich. Aber das kommt noch :wink:

Zitat:

Die Anzahl der Abtastwerte wird exakt vorgegeben, ich hab überlegt, ob ich nicht die vorhandenen Werte auf die nächste 2-er-Potenz extrapolieren kann...

Oder halt runterkürzen und gleiten.
Wenn du z.B. 1500 Werte vorgegeben hast. Dann fängst du halt bei dem ersten an und machst eine FFT über 1024 Werte, dan springst du z.B. 1 oder 10 Werte weiter und machst eine FFT über die Werte von 11 bis 1034. Ist die Frage, ob das zeitlich sich sehr nachteilig auswirkt und iwe man die Ergebisse am Ende wieder zusammenfasst (arithmetisch wäre der erste Ansatz).

[3_of_8]

Das wäre eine Idee. Ich teile den Zahlenbereich einfach auf. Das dürfte auf jeden Fall sehr viel schneller gehen als diese DFT, die hier auf meinem Notebook für 2000 Messwerte rund 1 Sekunde braucht. Ich arbeite gerade daran, die Ergebnisse bei der Multiplikation zu "puffern", was allerdings nicht viel bringen wird.

[gammatester]

Zitat:

Zitat von sirius

Zitat:

Als du von Betrag gesprochen hast

Das ist das, was ich meinte. Du kannst sicherlich hervorragend mit komplexen Zahlen umgehen. Dir fehlt aber einfach noch das Gefühl für diesen Zahlenbereich. Aber das kommt noch :wink:

Zitat:

Die Anzahl der Abtastwerte wird exakt vorgegeben, ich hab überlegt, ob ich nicht die vorhandenen Werte auf die nächste 2-er-Potenz extrapolieren kann...

Oder halt runterkürzen und gleiten.
Wenn du z.B. 1500 Werte vorgegeben hast. Dann fängst du halt bei dem ersten an und machst eine FFT über 1024 Werte, dan springst du z.B. 1 oder 10 Werte weiter und machst eine FFT über die Werte von 11 bis 1034. Ist die Frage, ob das zeitlich sich sehr nachteilig auswirkt und iwe man die Ergebisse am Ende wieder zusammenfasst (arithmetisch wäre der erste Ansatz).

Oder man nimmt eine FFT, die für fast beliebige Längen schnell ist und nur für Primzahllängen zur DFT entartet. Hier zwei Links für Delphi-Implementationen:

http://www-rab.larc.nasa.gov/nmp/nmp...tm#HerbsterFFT
http://www.simdesign.nl/fft.html

Wolfgang

[sirius]

Zitat:

Zitat von gammatester

Oder man nimmt eine FFT, die für fast beliebige Längen schnell ist und nur für Primzahllängen zur DFT entartet. Hier zwei Links für Delphi-Implementationen:
Wolfgang

Schau mal: Da waren wir bereits schon :wink:
DP - Problem bei FFT

[3_of_8]

Naja, mir ist grade aufgefallen, dass ich es durch ein bisschen tricksen auch mit einer FFT laufen lassen.

[gammatester]

Zitat:

Zitat von sirius

Zitat:

Zitat von gammatester

Oder man nimmt eine FFT, die für fast beliebige Längen schnell ist und nur für Primzahllängen zur DFT entartet. Hier zwei Links für Delphi-Implementationen:
Wolfgang

Schau mal: Da waren wir bereits schon :wink:
DP - Problem bei FFT

Ich weiß nicht ganz, vorauf sich dein Einwand bezieht. Mein Hinweis bezog sich auf FFTs, deren Längen keine 2-er-Potenz sind, und davon finde ich nichts (oder bin blind).

Wolfgang