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

[sirius]

Zitat:

Zitat von gammatester

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

Nee, da war ich blind :stupid:

Das ist ja interessant. Zugegeben, ich hab mich noch nie intensiver mit FFT-Algos beschäftigt,... ich habs schon mal runtergeladen. Das werde ich mir mal reinziehen.


Edit: Sehr schön der zweite Link
Aber Hanning machen die alle nicht (und in den Kurzfassungen beim ersten Link stehts auch nicht drin). Die Bass.dll tuts laut Beschreibung. Sollte man vielleicht noch ergänzen.

[3_of_8]

OK, ich hab hier noch ein Problem.

Ich suche mir momentan den Spitzenwert und rechne aus dem die Frequenz aus.

Meine Idee war jetzt, alle Werte der ersten Hälfte des fouriertransformierten Vektors mit der Frequenz, für die sie stehen, zu multiplizieren und dann am Ende das ganze durch alle Werte dieser Vektorhälfte zu teilen. Also in etwa so wie im Anhang.

Da bekomme ich Werte die ziemlich (aber nicht aber nicht ganz) gar nichts mit den wirklichen Frequenzen zu tun haben.

Also habe ich etwas eingebaut, was ich als NOISE_THRESHOLD bezeichne. NOISE_THRESHOLD ist eine Single-Konstante mit dem Wert 0.001. Werte aus dem Vektor gehen nur dann in die Berechnung ein, wenn sie >= NOISE_THRESHOLD*a[max] sind, wobei a mein Vektor ist und max der Index des Spitzenwertes. Das heißt, alle Werte, die kleiner sind als 1/1000 des Maximalwertes, werden vernachlässigt.

Dadurch erreiche ich eine Genauigkeit von +/- 1Hz, also statt 1420Hz erhalte ich 1419Hz.

Das klingt sehr gut, funktioniert aber nur bei Sinustönen.

Bei Square- oder Sägezahntönen stimmen die Werte überhaupt nicht.

[sirius]

Ich habs jetzt mehrmals durchgelesen und verstehe nicht, was du willst.

Deine Formel kann man übrigens noch kürzen, dann bleibt nur noch die Summe aller Frequenzen übrig, die Amplituden (a) verschwinden.

[3_of_8]

Bei meiner Formel hab ich mich verschrieben, ich hab die richtige Version wieder im Anhang. Bin mir aber nicht ganz sicher, ob man die nicht auch kürzen kann.

Ich will sozusagen den Durchschnitt aller Frequenzen haben, gewichtet nach Lautstärke und hoffe, damit die Frequenz meines Tons genau herauszufinden.

Soweit ich weiß, bekommt man doch die Frequenz des vorherrschenden Tons genau dadurch, dass man den Durchschnitt nimmt, oder?

EDIT: Die Formel lässt sich nicht weiter kürzen.

[sirius]

Hmm, da fehlen mir wohl ein paar Akustik vorlesungen. Keine ahnung. Ich dachte man hat ein eindeutiges Maximum, dass sich von allen anderen abhebt. Dem ist anscheinend nicht so. Kannst du mal so ein Spektrum zeichnen. Das würde mich mal persönlich interessieren.

Was du da mathematisch machst, da bin ich noch nicht ganz durchgestiegen. Aber du wichtest irgendwie höhere Frequenzen stärker als niedrige.

[3_of_8]

Nein.

Ich habe eine eindeutige Spitze, aber nur bei Sinustönen. Bei Sägezahntönen zum Beispiel habe ich die nicht, da habe ich etwa 3 solche Spitzen, die von der Amplitude her abnehmen, also eine sehr große, eine mittlere und eine kleinere.

Mathematisch gesehen mache ich eine einfache, gewichtete Durchschnittsbildung. Höhere Frequenzen werden natürlich nicht höher gewichtet, sie sind nur einfach größer und ziehen damit den Durchschnitt "in die Höhe". Man könnte es vergleichen mit einem Vektor für Notenverteilungen bei einem Test in der Schule:

Man hat einen Vektor (a1,...,a6), bei dem a1 die Anzahl der 1en und a6 die Anzahl der 6en ist. Die Durchschnittsnote des Tests bildet sich aus (a1*1+a2*2+a3*3+a4*4+a5*5+a6*6)/(a1+a2+a3+a4+a5+a6).

Das ist mathematisch gesehen das gleiche wie bei meiner Frequenzberechnung.

[sirius]

Mit der Mathematik hast du Recht. Ist für mich jetzt zu spät.

Aus Wiki: "Die Schwingung mit der niedrigsten Frequenz (gleichartige Schwingung der gesamten Saite) ist hier die Grundfrequenz, die anderen Schwingungen Oberschwingungen."
Soweit

Was ist den Sägezahn für ein Ton? Harmonisch ist der bestimmt nicht. Und er besteht auch nicht aus einem Grundton.
Was ist den z.B., wenn ich am Klavier drei Tasten anschlage, oder auf der Gitarre alle 6 Saiten schlage?
Für was entscheidet sich dann so ein Gitarrenstimmgerät?

[3_of_8]

*kratz*

Frag mich nicht.

Ein Sägezahnton hat halt statt der Sinusfunktion eine lineare Funktion.

Schwingung (rot) und FFT (grün) hab ich mal angehängt.

EDIT: http://de.wikipedia.org/wiki/Kippschwingung

[sirius]

Zitat:

Zitat von Wiki

Sägezahnschwingung als Überlagerung von Obertönen

wzbw (was zu beweisen war) :wink:
Der Sägezahn ist eine reine Überlagerung von Oberschwingungen.

Wenn also die Länge eines Zahns 0,2s dauert, also T=0,2s, dann ist f=5Hz deine Grundschwingung. Zu dieser Grundschwingung überlagern sich jetzt noch alle Harmonischen (bzw. Oberschwingungen), also alle n*5Hz (n..jede natürliche Zahl) Die Amplituden verkleinern sich allerdings mit jeder Oberschwingung (nur der Gleichanteil 0Hz kann etwas aus der Reihe tanzen, je nachdem wo der Sägezahn die Mitte hat).

Und das sieht man auch in deinem Bild: Eine Spitze immer im gleichen Abstand. die erste Spitze nach dem Gleichanteil ist die Grundschwingung.
Evtl. ist die Grundschwinung kein Vielfaches von Abtastfrequenz/(Anzahl der Abtastwerte), dann musst du die umliegenden Frequenzen irgendwie gruppieren.

btw: Wenn du ein Hanningfenster drüber legst, werden die Spitzen schärfer.

[3_of_8]

Ich schätze mal, ein Hanningfenster ist eine Art Filter. Wie funktioniert der?

EDIT: Aha. Das Ding heißt korrekt von Hann-Fenster und hat die Gleichung 0.5*(1+cos(2*Pi*n/M)) mit -M/2<=n<=M/2

Meine Implementation sieht so aus:

Delphi-Quellcode:

procedure vonHannWindow(var a: array of Single);
var I, M, n: Integer;
begin
  M:=length(a);
  for I:=0 to (M shr 1)-1 do
  begin
    n:=I-(M shr 2);
    a[I]:=a[I]*0.5*(1+cos(2*Pi*n/(M shr 2)));
  end;
end;

Ich erhalte klare Spitzen, kann aber immer noch keine Rechtecks/Sägezahntöne erkennen. Aber vielleicht geht das auch einfach nicht.

[sirius]

Nee, ein Filter ist es eigentlich nicht.

Es ist mehr so eine Funktion, die du überlagern (multplikativ) musst. Sagt dir "Faltung" etwas?

Die Idee dahinter ist, dass (wie ich schon mal erwähnte) du für deine Frequenzanalyse im Prinzip einen Ausschnitt verwendest.
dieser müsste allerdings problemlos unendlich mal hintereinenader passen ohne das Sprünge entstehen (so wie in fenster1 und fenster 2->Anhang)

Der Ausschnitt beträgt zufällig die Länge einer Farbe aus Fenster1. Wenn ich dieses Fenster wiedrum hintendrannhänge entsteht ein Sprung (Fenster2)

Dieser Sprung wirkt sich auf dein Frequenzspektrum aus.

Das Hanning-Fenster (ich sehe grad, das "bessere" heißt Hamming) aus Wiki sieht so aus:

Delphi-Quellcode:

for i:=0 to Abtastwerte-1 do //Abtastwerte ist die Länge des Arrays
    a[i]:=a[i]*(0.54+0.46*cos(2*pi*i/abtastwerte+pi));
  //und jetzt fft oder dft

Du kannst ja mal schauen, wie dein Zeitsignal vor und nach der Transformation aussieht.

Im übrigen musst du nach der FFT die Amplitudenwerte auf Grund der Fensterfunktion verdoppeln.



Aber ich sehe grade, dass Grouping dir eher hilft, nur da weis ich jetzt auch nicht viel...


Edit: zu deinem Edit: Ja fast, aber so sieht es bei mir zumindest richtig aus (bei vonHann Fenster steht da zweimal 0,5 statt 0,54 und 0,46)
Nur dein größeres Problem ist wahrscheinlich, dass dein Grundfrequenz kein Vielfaches von Abtastfrequenz/Abtastwerte ist. Und da hilft dieses Grouping. Vielleicht gibts auch noch andere Sachen. Aber du rechnest halt diskret (geht ja nicht anders) und dadurch bekommst du auch nur diskrete Frequenzen.

[3_of_8]

Aaah, erst das Fenster, dann die DFT... :lol: Ich habs andersrum

Delphi-Quellcode:

for i:=0 to Abtastwerte-1 do //Abtastwerte ist die Länge des Arrays
    a[i]:=a[i]*(0.54+0.46*cos(2*pi*i/abtastwerte+pi));
  //und jetzt fft oder dft

Ich dachte, man n (in dem fall i) müsste zwischen -M/2 und +M/2 liegen...

[sirius]

Zitat:

Ich dachte, man n (in dem fall i) müsste zwischen -M/2 und +M/2 liegen...

Das machte für mich keinen Sinn, da hab ichs einfach mal verschoben :wink: