Forum: Digitale Signalverarbeitung / DSP / Machine Learning Wie funktioniert dieser FFT-Algorithmus?


von Tobias P. (hubertus)


Lesenswert?

Hallo,
ich wollte mich (wieder einmal mehr....) mit dem Thema FFT auseinander 
setzen.
Anstatt zu versuchen, wie das ganze funktioniert, habe ich es mir 
diesmal etwas einfacher gemacht und unter

http://www.jjj.de/fft/

Die fix_fft heruntergeladen, eine FFT-Funktion mit Festpunktarithmetik.
Erstmal habe ich mit Visual Studio ein bisschen damit herumgespielt, das 
Programm scheint ansich zu funktionieren. Auch auf meinem Zielsystem, 
einem ARM7-Rechner, läuft es ohne Probleme recht fix durch (Performance 
habe ich allerdings noch nicht analysiert).
Aber, was mir nicht ganz klar ist:

Im Testprogramm der FFT wird eine komische "Umsortierung" der Samples 
vorgenommen:
1
for(i=0; i<N; i++)
2
{
3
  x[i] = AMPLITUDE*cos(i*FREQUENCY*(2*3.1415926535)/N);
4
  if (i & 0x01) /* wozu braucht es diese if-Abfrage?! */
5
    fx[(N+i)>>1] = x[i];
6
  else
7
    fx[i>>1] = x[i];
8
}

Mir ist klar, dass man bei der FFT erst diesen Bitreversal-Algorithmus 
durchführen muss; aber in der eigentlichen FFT-Funktion wird auch etwas 
derartiges gemacht. Im Kommentar heisst es lediglich "re-order data", 
was auch immer das bedeuten mag.

Meine Frage jetzt:

1. Wozu ist diese if-Abfrage in obigem Codeschnipsel nötig? Kann ich der 
FFT-Funktion nicht einfach die Rohdaten übergeben, welche ich von meinem 
AD-Wandler kriege?

2. Wie kann ich testen, ob das, was aus der FFT raus kommt, 
einigermassen korrekt ist? Ich habe mal versucht, ein DC-Signal 
einzuspeisen; das funktioniert, aber nur wenn ich die Umsortierung der 
Samples wie in oben gezeigtem Codeschnipsel mache. Sonst kommt Murks 
raus.
Testweise habe ich auch alle Samples auf 0 gesetzt, bis auf das erste, 
das den Wert 1000 erhielt. Nach durchführen der FFT über 256 Samples 
hatten dann die ersten 128 Ausgangswerte den Wert 78, alle anderen 0. 
Ist das jetzt richtig oder nicht?
Alles in allem zwar eine gute FFT-Routine, aber schlecht dokumentiert.

Hier noch der Link zum Sourcecode: http://www.jjj.de/fft/fix_fft.tar.gz

Wie gesagt - lässt sich problemlos compilieren und scheint auch zu 
funktionieren; nur weiss ich nicht ob es richtig rechnet.

Gruss Tobias


PS: was mir noch einfällt: lässt sich diese riesen grosse Sinustabelle 
zwecks Speicherersparnis noch etwas verkleinern? Wie ungenau wird die 
FFT, wenn ich weniger Samples in der Sinustabelle habe?
Und lässt sich das ganze noch etwas optimieren?

von Detlef _. (detlef_a)


Lesenswert?

Tobias Plüss schrieb:

> 1. Wozu ist diese if-Abfrage in obigem Codeschnipsel nötig? Kann ich der
> FFT-Funktion nicht einfach die Rohdaten übergeben, welche ich von meinem
> AD-Wandler kriege?

Ich kann mir den link zum sourcecode nicht anschauen, deswegen 
Vermutung: Normalerweise bildet eine N-Punkte FFT N komplexe 
Eingangswerte auf N komplexe Spektralwerte ab. Möglicherweise handelt es 
sich bei Deinem Beispielsourcecode um eine FFT, die mittels einer N/2 
FFT eine Anzahl von N reellen Eingangswerte (AD- Wandler) auf N 
komplexe Spektralwerte transformiert. Das Spektrum reeller Eingangswerte 
hat Symmetrieeigenschaften, die das erlauben. Dieses Vorgehen ist 
Standard.

> 2. Wie kann ich testen, ob das, was aus der FFT raus kommt,
> einigermassen korrekt ist? Ich habe mal versucht, ein DC-Signal
> einzuspeisen; das funktioniert, aber nur wenn ich die Umsortierung der
> Samples wie in oben gezeigtem Codeschnipsel mache. Sonst kommt Murks
> raus.

Dann ist doch alles gut. Die FFT von DC: alles 0 ausser dem ersten bin.

> Testweise habe ich auch alle Samples auf 0 gesetzt, bis auf das erste,
> das den Wert 1000 erhielt. Nach durchführen der FFT über 256 Samples
> hatten dann die ersten 128 Ausgangswerte den Wert 78, alle anderen 0.
> Ist das jetzt richtig oder nicht?

Das ist die FFT eines Impulses, da sind alle bins belegt und die 
Ergebnisse sind rein reel, also passen die Nullen gut. Die 78 könnten 
auch passen: Parseval: 1000^2~sqrt(2)*78^2*128.  Die Reihenfolge der 
Ergebnisse ist etwas ungewöhnlich.
>
> PS: was mir noch einfällt: lässt sich diese riesen grosse Sinustabelle
> zwecks Speicherersparnis noch etwas verkleinern? Wie ungenau wird die
> FFT, wenn ich weniger Samples in der Sinustabelle habe?

Die Größe der Sinustabelle hängt von der Länge der FFT ab. Ein 
Viertelbogen eines Sinus reicht aber, bei ner 256er rellen FFT 
256/2/4=32 Tabellenwerte.

Schau dir an: Brigham, The Fast Fourier Transform, gibts auch auf 
Deutsch.

math rulez!

Cheers
Detlef

von Tobias P. (hubertus)


Angehängte Dateien:

Lesenswert?

Hi Detlef,
anbei mal der Sourcecode, damit du dir ein Bild machen kannst, von was 
ich rede.
Diese komische Sample-Sortierung ist in der test_fix_fft.c zu sehen, und 
zwar ab Zeile 28. Schaus dir mal an, vielleicht kannst du mir dann 
sagen, wozu man das braucht?

> math rulez!
dito; leider ist mir das ganze FFT-Zeug ein bisschen zu hoch. Nirgends 
wird es verständlich erklärt; ich kann zwar Integrale und Differentiale 
berechnen, aber um eine FFT zu verstehen scheint das noch nicht wirklich 
zu reichen. :-(

Gruss Tobias

von Detlef _. (detlef_a)


Lesenswert?

>>Nirgends wird es verständlich erklärt;
Der Brigham erklärts verständlich, das hat bei mir den Knoten des 
Unverständnis durchschlagen.

>>Schaus dir mal an, vielleicht kannst du mir dann sagen, wozu man das
braucht?

Wie vermutet, eine FFT für rein reele Eingangsdaten. Das steht doch in 
fix_fft.c über der Funktion fix_fftr schön erklärt (das 'r' ist von 
'reel'). Dieses Zusammenwürfeln der Daten findest Du auch im Brigham.

BTW: besser erstmal mit Matlab/Octave bißchen FFT hin und her 
transformieren um nen Gefühl zu kriegen. Die integer-FFT hat noch nen 
paar Spezialitäten wie Overflowproblematik und Rundungsrauschen, diese 
Dinge vom reinen Verständnis besser separieren.

Die Riesentabelle in der Source ist nicht so besonders effizient. Die 
ist auf die maximale Länge von 1024 ausgelegt, Du brauchst bei 256er 
Länge nur jeden vierten Wert. Mal die Source analysieren, da läßt sich 
was abspecken wenns zu groß wird.

Aber es funktioniert ja, es ist alles gut.

Cheers
Detlef

Bitte melde dich an um einen Beitrag zu schreiben. Anmeldung ist kostenlos und dauert nur eine Minute.
Bestehender Account
Schon ein Account bei Google/GoogleMail? Keine Anmeldung erforderlich!
Mit Google-Account einloggen
Noch kein Account? Hier anmelden.