Problem bei FFT

[sirius]

Zitat:

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

Ich nehme an, dass ich nicht mehr direkt antworten brauche.

Trotzdem noch ein paar Bemerkungen zur FFT/DFT: Da du sie anscheinend das erste Mal verwendest.

Zitat:

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.

Mit dem Anhängen von Nullen verfälschst du das Ergebnis. Besonders der Knick an der Stelle, wo die Nullen anfangen, bringt dir Frequenzen ins Spektrum, die du so nicht hättest. Teste das mal an dem einfachen 10Hz-Sinus.
Dem kannst du mit einer DFT aus dem Weg gehen. Die ist ein bisschen einfacher zu programmieren. Und in Zeiten von 3GHz-Prozessoren dürfte der Zeitunterschied sehr gering sein.

Zitat:

2. Problem: Erhalten einer Spitze. Ich nehme hierfür einfach den Maximalwert aus dem Array.

Da weis ich nicht, warum und wofür.

Zitat:

3. Problem: Umrechnen des Wertes in eine Frequenz.

Hast du anscheinend hinbekommen. Ansonsten stehts ein paar Beiträge weiter oben.

Zitat:

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.

Ok, Anwendung der FFT/DFT:
Du hast einen bestimmten Zeitabschnitt eines Signals und willst wissen, welche Frequenzanteile mit welcher Amplitude und welcher Phasenlage in dem Signal sind. Erweiterung für unser obiges Bsp:
for n:=0 to 1023 do a[n]:=1+5*sin(2*pi*10*n/1000)+8*sin(2*pi*80/1000); So jetzt hast du neben den 10Hz auch noch einmal 80Hz mit Amplitude 8 und 0Hz (also Gleichanteil) von 1.
Das sollte nach der FFT auch wieder herauskommen.
In dieser Form kann man quasi jedes Signal aus seinen einzelnen Schwingungen zusammensetzen. Das kann auch ein Rechtecksignal sein, oder ein Dreieck, oder eine e-Funktion. Kannst ja mal alles testen, welches Spektrum da so herauskommt.
Zu dem Beispiel kann man noch ein Phasenlage hinzubringen (also den Sinus etwas verschieben):
for n:=0 to 1023 do a[n]:=1+5*sin(2*pi*10*n/1000+pi/2)+8*sin(2*pi*80/1000); Ich hab die 10Hz mal um 90° verschoben. Das bekommst du aus der FFT heraus, wenn du die Differenz der Werte von arctan(Im/Re) auswertest.

Dementsprechend kannst du deinen Auswertezeitraum so festlegen, dass er von der Länge her einer Zweierpotenz entspricht(z.B. 1024 Werte). Wenn da dein Sinuston mit f Hz drinn ist, solltest du den auch im Sprektrum sehen, neben vielen anderen Tönen.
Jetzt hast du zwei Möglichkeiten
-gleitende FFT: Du verschiebst dein FFT-Fenster nicht um 1024 Werte weiter für die nächste Auswertung, sondern nur um 1, 10 oder 100
-nicht gleitend: Du verschiebst es halt um 1024 Werte



Jetzt etwas Grundsätzliches:
Eine Fourieranalyse ist etwas ganz Tolles. Aber sie ist nicht der Weisheit letzter Schluss. Sie ist eigentlich nur der Einstieg in die Signalanalyse. Folgende Probleme bei der FFT fallen mir spontan ein:
1. Die FFT ist natürlich diskret (deswegen haben wir nur für bestimmte Frequenzen den richtigen Wert, alle anderen Frequenzen verteilen sich mit auf die Umliegenden (wie wir bei 10Hz und 1024 Abtastwerten ja gesehen haben)
Lösung: Du machst für Frequenzen, die du exakt haben möchtest eine eigene DFT, wo die Anzahl der Abtastwerte ein Vielfaches der zu untersuchenden Frequenz ist. (Geht natürlich nur, wenn du auf bestimmte Frequenzen testen willst/kannst.)
2. Wenn du dir aus deinem Signal die 1024 Werte nimmst, bekommst du nur gute Ergebnisse, wenn du dieses Fenster quasi x-mal aneinander hängst und danach keine Sprünge im Signal sind (an den Fenstergrenzen).
Lösung:Fensterfunktionen
3. Du wirst es nie schaffen, alle Frequenzen exakt aus dem Signal herauszubekommen (ausser es ist so sauber, wie in unserem Beispiel), dass muss dir bewusst sein. Es befassen sich x Bücher und dahinter Wissenschaftler mit diesem Thema der Signalanalyse. Die FFT ist ein Werkzeug, aber kein Perfektes. Prony-Methode ist z.B. ein anderes. die Methode geht dann über Differentialgleichungen. Und hat natürlich auch wieder Nachteile.


Edit:
-Grobe Rechtschreibfehler korrigiert
-Für die Analyse von Musik um bunte Bildchen darzustellen, à la Visualisierung in diversen Media-Playern reicht die FFT(DFT) aus.

Edit2: Hier ist auch noch ein schöner Link, um bildhaft zu zeigen, was Knicke im Signal bedeuten (deine Nullen oder falsche Fensterung)

[negaH]

Du möchtest nur wenige und vorher bekannnten Frequenzen aus deinem Signal auswerten ?

Dann gibt es noch andere Möglichkeiten effektiver zum Ziel zu kommmen.

- IIR, FIR Filter können Banpässe,High/Lowpass, Notch und andere Filter sein. Sie dämpfen oder verstärken quasi ein definierbares Frequenzband
- Digitale Lock-In-Amplifier. Das sind sogenannte PSD->Phase Sensitive Devices oder auch Synchrone Gleichrichter genannt. Es gibt viele Namen für diese Methode. Sie können aber eine exakte Frequenz/Codesequenz/Spreadspektrum herausfiltern bzw. den Gleichstromanteil und Phasenlage dieser gesuchten Frequenz im vorhandenen Signalgemisch ausrechnen. Damit kann man zb. aus so stark verrauschten/gestörten Signalen die Nutzdaten/Frequenzen noch rausfinden selbst wenn sie nur noch Bruchteile in der Amplitude betragen. Fast jede neuere Funktechnik basiert auf solchen Methoden, zb. Bluetooth, ZigBee, GSM, GPS, Radar wären ohne nicht möglich.

Hier im Forum habe ich schon mal erklärt wie man in Software einen solchen Lookin Amplifier baut, suche mal danach. Ich selber habe damit DTMF Tonsequenzen (töne beim Wählen am Telefon) dekoriert, in Echtzeit. Der Vorteil ist dabei das diese Methode präziser und performanter als eine FFT ist. Eine FFT als Spekrumanalyse ist immer ein Breitbandgeschoss.

Die Frage ist also ob du wirklich das komplette Spektrum benötigst oder nur ganz spezielle Frequenzen. Je nach Antwort dann FFT oder Signalprozessing mit Filtern.


Gruß Hagen

[sirius]

Hallo Hagen,

Der Lookin-Verstärker macht doch nix anderes als DFT auf einer Frequenz.
Du multiplizierst dein Eingangssignal mit einem Referenzsignal bekannter Frequenz. Genau dasselbe macht die DFT. Die macht es allerdings iterativ für alle Frequnezen die in das Zeitfenster passen.

Damit ist natürlich, der zeitliche Aufwand erheblich kürzer gegenüber einer DFT. Außerdem, und das ist das Wichtigste, lässt sich ein Look-In-Verstärker recht einfach hardwareseitig zusammenbauen (im Vergleich zu einer FFT) und ist, wie auch die Softwarevariante, sehr fix.

[negaH]

Hi Sirius,

ich habe nicht das Gegenteil behauptet Auch die IIR/FIR Filter lassen sich auf die Theorien der Fourier Transformation zurückführen.

Entscheidend ist aber was anderes: Wie versteht ein Laie/Mensch am schnellsten eine Technologie ? Und da muß ich sagen das mir persönlich die Theorien der Lookin-Amplifiers, PSDs, FIR, IIR viel leichter nachvollziehbar waren als die der FFT. Obwohl ich erst mit der FFT angefangen habe.

Es gibt aber auch wesentliche Unterschiede. Ein Lookin-Amplifier kann zb. mit jeder dynamischen Frequenzform arbeiten. Dh. statt einer Sinusfrequenz mit fester Frequenz im Zeitspektrum kann ein Lookin Amlifier/PSD mit einer beliebigen Signalform, zb. ein Code basierend auf einem Bakercode, arbeiten. Das benutzt zb. die Radartechnoogie oder das Speadspektrum im WLAN. Sowas per FFT machen zu wollen dürfte sehr schwierig sein, ich wüsste jetzt nicht wie das machbar ist.

Gruß Hagen

[3_of_8]

Ich möchte eine große Anzahl von Frequenzen auf einem Spektrum von ~5000-8000 Hz überprüfen. Also ist das wohl eher ungeeignet.

[sirius]

Zitat von negaH:

ich habe nicht das Gegenteil behauptet

Ich weis.

Zitat von negaH:

Auch die IIR/FIR Filter lassen sich auf die Theorien der Fourier Transformation zurückführen.

Wenn man es so nimmt, lässt sich auch Prony darauf zurückführen (ich habs selber noch nicht ausprobiert, wurde mir aber gesagt)
Allerdings ist ein (Dual-) Look-In-Verstärker exakt DFT. Ob man nun mit 2 sinussen (Was ist eigentlich die Mehrzahl von sinus?) mit einer Phasen-Verschiebung von 90° multipliziert oder mit einer Komplexen Zahl, die sich mit derselben Frequenz auf dem Einheitskreis dreht... Es ist wirklich haargenau dasselbe. Nur die DFT macht es als hintereinanderweg für alle Frequenzen.

Zitat von negaH:

Entscheidend ist aber was anderes: Wie versteht ein Laie/Mensch am schnellsten eine Technologie ? Und da muß ich sagen das mir persönlich die Theorien der Lookin-Amplifiers, PSDs, FIR, IIR viel leichter nachvollziehbar waren als die der FFT. Obwohl ich erst mit der FFT angefangen habe.

Und genau deswegen hab ich es nochmal ergänzt, dass der LookIn-Verstärker eben die DFT für eine Frequenz ist. Weiter oben hab ich ja auch geschrieben, dass zum Verständnis des ganzen, eine Implementierung einer DFT besser ist, als gleich eine FFT zu nehmen. (Was ja in einem anderen Beitrag/Thread auch beherzigt wurde)

Zitat von negaH:

Es gibt aber auch wesentliche Unterschiede. Ein Lookin-Amplifier kann zb. mit jeder dynamischen Frequenzform arbeiten. Dh. statt einer Sinusfrequenz mit fester Frequenz im Zeitspektrum kann ein Lookin Amlifier/PSD mit einer beliebigen Signalform, zb. ein Code basierend auf einem Bakercode, arbeiten. Das benutzt zb. die Radartechnoogie oder das Speadspektrum im WLAN. Sowas per FFT machen zu wollen dürfte sehr schwierig sein, ich wüsste jetzt nicht wie das machbar ist.

Natürlich kann ich auch eine komplexe Funktion erstellen, die im Zeitberaich jedes beliebige Signal nachbildet und das durch eine (dann Pseudo)-DFT schicken. Und wenn das Signal nicht ganz chaotisch ist, ergeben sich auch wieder Dopplungen, sodass man (à la FFT) abkürzen kann. Aber hier ist, wie du auch schon sagtest, der Unterschied, ob ich ein Spektrum haben will, oder nur eine Frequenz/Frequenzfolge. Und dann brauch ich auch keine (Pseudo-)FFT entwickeln. Und wahrscheinlich gibts sowas schon, nur unter einem ganz anderen Namen....

[3_of_8]

Das ist meine momentan (funktionierende) Implementation einer DFT. Wenn möglich wird eine FFT angewandt. Ich empfehle, das ganze nicht mit mehr als 2000 Werten aufzurufen, es sei denn, die Anzahl der Werte ist eine Zweierpotenz oder lässt sich sehr oft durch 2 teilen.

Aufgerufen wird das einfach mit DFT(a), wobei a ein Array of Integer, Single oder ein TComplexArray ist.

Nach dem Aufruf enthält das Array die Amplitudenwerte der einzelnen Frequenzen, in Index 1 steht der Wert der Frequenz (Samplingrate/Abtastwerte) Hz, in Index 2 2*(Samplingrate/Abtastwerte) Hz und in Index I steht der Wert für I*(Samplingrate/Abtastwerte) Hz. Man sollte nur Werte bis zum Index (Abtastwerte/2), also bis (Samplingrate/2) Hz, der sogenannten Nyquist-Frequenz abfragen, da einem ansonsten der Alias-Effekt auftritt.

(Vielleicht wäre das was für die CL?)

Delphi-Quellcode:

unit complex;

interface

uses Math;

type
  TComplex=record
    re, im: Extended; 
  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 DFT(var a: TComplexArray); overload;
procedure DFT(var a: array of Single); overload;
procedure DFT(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);
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 DoDFT(var a: TComplexArray; n, lo: Integer; w: TComplex);
var b, f: array of TComplex;
    I, J, index: Integer;
begin
  setlength(b, n);
  for I:=0 to n-1 do
  begin
    b[I]:=a[I+lo];
    a[I+lo]:=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
  begin
    index:=0;
    for J:=0 to n-1 do
    begin
      a[I+lo]:=AddC(a[I+lo], MulC(b[J], f[index]));
      inc(index, I);
      if index>n then dec(index, n);
    end;
  end;
end;

procedure DoFFT(var a: TComplexArray; n, lo: Integer; w: TComplex);
var I, m: Integer;
    z, v, h: TComplex;
begin
  if n and (n-1)=0 then
  begin
    if n>1 then
    begin
        m:=n shr 1; 
        z:=MakeC(1, 0);
        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 else if n>1 then
  begin
    DoDFT(a, n, lo, w);
  end;
end;

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

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

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

end.

EDIT: Hier noch eine Funktion zum Analysieren der vorherrschenden Frequenz in einem Tonabschnitt.

Delphi-Quellcode:

function GetFrequency(Data: array of Single; SampleRate, SampleIndex, 
    Samples: Integer): Integer;
var I, n, max: Integer;
    a: array of Single;
begin
  max:=0;
  if length(Data)-SampleIndex<Samples then
    raise EAccessViolation.Create('Sound index or length out of bounds.');
  setlength(a, Samples);
  for I:=0 to Samples-1 do
    a[I]:=Data[SampleIndex+I];
  dft(a);
  for I:=0 to (n shr 1)-1 do if a[I]>a[max] then
  begin
    max:=I;
  end;
  Result:=trunc(SampleRate/Samples*max);
end;

[sirius]

Klappt das? Die DFT so in die Rekursion der FFT zu basteln?

Edit: Ja, müsste Edit2: Man könnte es natürlich auch vor der Rekursion überprüfen.

Edit2: zu GetFrequency
Und immer daran denken, dass der absolute Fehler (durch digitalisierung) bei SampleRate/Samples/2 liegt
Von anderen Fehlern habe ich schon genug gesprochen, dürften aber für diese (weniger kritische) Analyse nicht interessant sein.

[SunBlack]

Hallo 3_of_8,

ich hab mir gerade deinen Code angeguckt um zu verstehen, wie die DFt genau funktioniert (mit Wiki werde ich da nicht schlau ). Allerdings habe ich da ein Problem. Wie soll man den real- und Imaginaryteil verstehen? Also kann man die Daten irgendwie weiter gebrauchen, oder bringt einem immer nur der Pythagoras da etwas?

Gruß
SunBlack

[sirius]

Es ist schon wichtig, dass es eine komplexe Zahl ist (kommt ja auch aus dem algorithmus).
Eine komplexe Zahl besteht immer aus Betrag und Phasenwinkel. Den Betrag bekommst du ja hier aus dem Pythagoras. Den Winkel aus dem Tangens (Arcustangens).
Letztenendes bringt die ja die DFT den Anteil des Sinus für jeden Frequenzanteil. Und der kann ja auch verschoben sein, dass ist dann der Winkel.

Jedes beliebige Signal f setzt sich ja aus der Summe:

Code:

f(t) = A_1 * sin (w_1*t + phi_1) 
     + A_2 * sin (w_2*t + phi_2) 
     + A_3 * sin (w_3*t + phi_3) 
     + ...

zusammen. dabei sind die w_x deine Frequenzen, A_x die dazugehörigen Amplituden (aus Pythagoras) und phi_x der Phasenwinkel (=arctan(Imag/Real))

Edit: Achtung: Den winkel musst du auch noch wie den Betrag bezüglich des Aliasing korrigieren.