Problem bei FFT

[sirius]

Ein sinus ist für den Anfang genau Richtig, später macht man dann einfach mal die Addition von Schwingungen unterschiedlicher Frequenz.
Aber welche Frequenz hat dein Sinus.
Du musst zum Testen dir das auch alles genau festlegen (128 Werte sin etwas wenig).
Versuche mal mindestens 5 Perioden (10 ist besser) bei einer Abtastfrequenz von mindestens dem 3-fachen hineinzubekommen.
Dazu musst du dir auch genau festlegen, welche Zeitwerte du hast und kannst damit dann die diskreten Frequenzwerte nach der FFt ausrechnen und dann gilt eigentlich (egal welcher FFT-Algo oder eben nur DFT):
-Das erste Element ist der Gleichanteil
-Und alles ab dem zweiten Element ist der entsprechende Frequenzanteil; allerdings nur die Hälfte, da die andere Hälfte im hinteren Teil des Spektrums ist, gespiegelt an Fs/2 (<<-Aliasing); Der Winkel ist konjugiert Komplex.

[3_of_8]

Erklär mir mal, wie ich meine FFT mit Parameter eines Arrays aufrufen muss und wie ich die Werte interpretieren muss, Beispieldaten wähl bitte du, weil ich grad keine Ahnung mehr hab.

Übrigens hab ich oben geeditet.

[sirius]

Ich vermute immernoch, dass dein w=e^(-2*i*pi/n) ist. Bedenke, dass die n-te Wurzel aus x immer n Lösungen hat. Deine ist auch eine, aber für den Fall wahrscheinlich nicht die "Richtige".


Ok, Eingangsdaten:

Nimm mal ein array von 1024 Werten:

Delphi-Quellcode:

for n:=0 to 1023 do a[n]:=5*sin(2*pi*10*n/1000);

Das entspricht A*sin(w*t)
A..Amplitude
w..omega/Kreisfrequenz=2*pi*f=2*pi*10
t..entspricht n/1000;wobei 1000 die Abtastfrequenz ist

Du bekommst also eine Sinuskurve mit 10 Hz, wo aller 1ms ein fiktiver Messpunkt ist; also 100 Mespunkte pro Periode -->für 1024 Werte etwa 10 Perioden.

Nach der FFT entsprechen die Absolutwerte der Elemente folgenden Frequenzen: n*1000/1024 (n=0..1023); also an 11. Stelle liegen etwa unsere 10Hz und dort sollte eine Spitze im Spektrum auftauchen, wie auch etwa bei Stelle 1013 (Aliasing).

Die Winkel kannst du natürlich bei einer Schwingung vergessen, da ja keine Referenz existiert.

[3_of_8]

Zitat:

Ich vermute immernoch, dass dein w=e^(-2*i*pi/n) ist. Bedenke, dass die n-te Wurzel aus x immer n Lösungen hat. Deine ist auch eine, aber für den Fall wahrscheinlich nicht die "Richtige".

Eben gerade nicht.

Mein w ist keine Lösung von w^n=1, denn w^n=1+0.1i.

Seltsamerweise ist bei mir folgendes:

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

Für n=1 ist w^1=1
Für n=2 ist w^2=1
Für n=3 ist w^2=~(-0.96-1.21i)

[sirius]

Für welches w? Ich verstehe nicht ganz, was du rechnest.

Dir fehlt anscheinend noch ein bisschen Wissen zu komplexen Zahlen. Das kann ich dir aber hier und mal eben nebenbei nicht erklären.

Im Anhang findest du wie es die DFT macht, für n=8 Abtastwerte
(Der Rote Pfeil ist w)
Da der Betrag der Zahl immer 1 ist, werden beim potenzieren nur die Winkel verändert (und zwar werden die addiert)
Dasselbe Ergebnis kommt bei w=e^( +2*pi*i/8 ) heraus.

[3_of_8]

Ich weiß meiner Meinung nach genug über komplexe Zahlen.

Mein Algorithmus für die FFT erwartet ein w, das die n-te primitive Einheitswurzel ist, wobei n die Länge des Arrays ist.

Die n-te primitive Einheitswurzel lässt sich berechnen durch e^(2*Pi*i/n).

e^(x*i) ist nach der 1. Eulerschen Formel gleich cos(x)+sin(x)*i.

Daher schließe ich w=e^(2*Pi*i/n)=cos(2*Pi/n)+sin(2*Pi/n)*i. Mit dieser Formel berechne ich w, allerdings stimmen die Werte seltsamerweise nicht.

EDIT: Meine Formeln sind ein eindeutiger Schluss aus deinem Kreis.

EDIT2: :wall: Mein Wert für w war die ganze Zeit richtig, das Problem war meine Testroutine, die einfach falsch potenziert hat. Und das DegToRad() muss ich, wie ich gleich vermutet hatte, weglassen.

[sirius]

So, jetzt habe ich mal ein wenig geanuer reingeschaut

Ich bin über die Variante (array of single) reingegangen
sieht jetzt so aus:

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

Das Minus ist egal, da wir nur den Betrag auswerten.
Dein Fehler ist, dass du nur den Realteil auswertest.

Getestet mit Tchart:

Delphi-Quellcode:

procedure TForm1.Button1Click(Sender: TObject);
var l:Tlineseries;
    a:array[0..1023] of single;
    i:integer;
begin
  for i:=0 to 1023 do a[i]:=5*sin(2*pi*10*i/1000);
  for i:=0 to chart1.SeriesCount-1 do chart1.Series[0].free;
  l:=tlineseries.create(chart1);
  chart1.AddSeries(l);

  fft(a);
  for i:=0 to 1023 do l.AddXY(i*1000/1024,a[i]);
end;

Wenn du aus den beiden 1000 auch eine 1024 machst, bekommst du ein "scharfes" Ergebnis.

[3_of_8]

Ich habe keine Chart gefunden und hab mir daher mit einem Image beholfen:

EDIT: Wie oben schon beschrieben, habe ich ja das Problem gefunden. ;)

Der FFT scheint zu funktionieren, jetzt die Frage, warum der Ausschlag so ist, wie er ist.

EDIT2: Ich weiß, warum er bei 10/11 so ist, aber warum ist er hinten auch so? Du hast vorher was von (Anti-)Aliasing gesagt, was meinst du damit? Und ist es normal, dass einige Werte nach dem FFT negativ sind?

[sirius]

Ich hab erstmal mein Ergebnis angefügt.

Hab übersehen, dass du den Fehler gefunden hast. Nimm aber bitte trotzdem statt dem Realteil den Absolutwert, ansonsten bekommst du bei komplexeren Signalen Probleme. Ausserdem werden dann alle deine Werte positiv, so wie es sein soll.
Und ob w=w1=e^(-2*pi*i/n) oder w=w2=(+2*pi*i/n) entscheidet am Ende das Vorzeichen der Winkel der einzelnen Frequenzanteile und ist somit schon wichtig. Und wie ich schon sagte, w1^n=w2^n=1.

Wenn dein (Ton-)signal mein beschriebenes 10Hz Signal ist, verstehe ich nicht, wo in deinem Bild die Frequenzanteile unter 10 Hz sind, die müssten wieder Richtung 0 gehen.

Und für das Aliasing müsste ich eigentlich ziemlich weit ausholen, und mir fehlt die Zeit dazu.

Also kurz (und damit alle Angaben ohne Gewähr):
Aliasing entsteht durch die Digitalisierung. Nehmen wir an du hast ein 10Hz-signal analog und tastest es Digital mit 50Hz ab, also aller 20ms ein Wert.
Und jetzt stehst du mit deinem abgetastetem Signal da, aller 20ms ein Wert und überlegst, wie das Ursprungssignal aussah. du kannst es mal ausprobieren, in dein abgetastetes Singnal passen folgende Frequenzen einzeln hervorragend rein:
10Hz
40Hz
60Hz
90Hz
110Hz
140Hz
160Hz
.
.
.
also alle 50Hz*k-10Hz und 50Hz*k+10Hz (k..ganze Zahl)

Nun musst du dich für eine Frequenz entscheiden.
Das geht am einfachsten, wenn du weist, was im Ursprungssignal überhaupt drinn war.
Wenn du z.B. wüsstest, dass im Ursprungssignal nur Werte zwischen 0<=f<25Hz waren, ist das Ergebnis klar. entsprechend kannst du das Fenster um 25Hz immer weiter verschieben.

Deswegen gibt es Antialiasingfilter, die machen nix weiter als das analoge Eingangsignal in einen Bandpass (oder Tiefpass) zu quetschen, damit beim A/D-wandeln nur signale (in unserem eispiel) zwischen 0 und 25 Hz auftauchen (oder eben ein anderes Frequenz-Band).
Man muss also immer mit der doppelten Frequenz abtasten, wie das Frequenzband breit ist (am besten mehr).

Nun zurück zu unserr FFT. Wir setzen ja beim bereits abgetasten Signal an und versuchen die Ursprungsfrequenzanteile vom analogen Signal zu bekommen.

In dem von mir angegebenen Beispiel haben wir ein mit 1000Hz abgetastetes 10Hz Signal. Damit gilt die Bedingung -> Unser Frequenzband ist 0 bis 10Hz =10Hz und wir tasten mit (deutlich) mehr als 2*10Hz=20Hz ab.
Wenn wir mit 1000Hz abtasten spiegeln sich, wie oben beschrieben, alle anderen (Aliasing)-Frequenzen um 1000Hz, also
10Hz
990Hz
1010Hz
1990Hz
.
.
.
Die FFt rechnet aber nur bis 1000Hz und teilt die Anteile gerecht und fair, gleichmäßig auf beide Anteile auf. Demnach haben wir bei einer Ursprungsamplitude von 5 bei 10Hz, nachher einen Frequenzanteil von 2,5 bei 10Hz und 2,5 bei 990Hz.

Soweit, was mir dazu einfällt.
Kannst das Alising mal damit testen, dass du statt den 10Hz, die anderen frequenzen testest. Es sollte das gleiche Ergebnis liefern. du kannst also nur bis Abtastrate/2 (=500Hz) die Ergebnisse verwenden.

Edit+Edit2: Aliasing-Bild angefügt

[3_of_8]

Folgendes Problem:

Ich habe eine WAV-Datei geparst.

Sie enthält einen Ton der länge t ms, die Abtastfrequenz ist a Hz. Damit ist die Anzahl der Abtastungen n=t*a/1000.
Die WAV-Datei enthält einen Sinuston von f Hz.

1. Problem: Eine FFT kann nur mit 2^m (0<=m<unendlich) Abtastungen arbeiten. Meine Lösung: Anhängen von Nullen, bis die Zahl der Abtastungen eine Zweierpotenz ist.
2. Problem: Erhalten einer Spitze. Ich nehme hierfür einfach den Maximalwert aus dem Array.
3. Problem: Umrechnen des Wertes in eine Frequenz.

Hier habe ich jetzt das Problem, dass ich nicht weiß, wie ich das umrechnen soll. Meine erste Lösung wäre gewesen Arrayindex*a/length(Array). Das funktioniert erwartungsgemäß überhaupt nicht und ich erhalte unsinnige Werte, die sich außerdem verändern, je nachdem, wie lange der Abschnitt des Sinustons ist, den ich analysiere.

Wo liegen meine Fehler?

EDIT: Heureka, es funktioniert. Ich geb euch allen ne Runde an der DP-Bar aus.