Entwurf: Spektralanalyse mit der FFT

Baustelle

Dieser Artikel wird gerade erstellt oder es werden umfangreiche Änderungen daran vorgenommen. Die Informationen können noch unvollständig oder inkorrekt sein.

Bevor du größere Änderungen vornimmst, stimme Dich bitte über die Diskussionsseite mit den Autoren darüber ab.

Oft findet man im DSP-Forum Fragen wie: "Ich habe 1024 Signalpunkte gespeichert und darauf eine FFT ausgeführt, was mache ich mit dem Ergebnis?"

Da muss man zuerst mal die Gegenfrage stellen: was willst du überhaupt mit dem Ergebnis machen?

In der Regel geht es dem Frager darum, ein Leistungsdichtespektrum darzustellen, also einen Plot von der Signalleistung über der Frequenz, ähnlich wie es eine Audio-Spektrumanalysator anzeigt.

Was macht eine FFT?

FFT (Fast Fourier Transform) ist eigentlich nur eine Bezeichnung für einen Algorithmus zur schnellen Berechnung der DFT (Discrete Fourier Transform) mit bestimmten Randbedingungen. Da die DFT sehr oft mit dem FFT-Algorithmus berechnet wird, verwenden viele beide Begriffe gleichwertig.

Die DFT/FFT berechnet, vereinfacht ausgedrückt, welche Frequenzen im Originalsignal enthalten sind, welche Sinusschwingungen also mit welcher Phase addiert werden müssten, um auf das ursprüngliche Signal zu kommen. Konkret bekommt man für jeden Frequenzpunkt (0 bis Abtastfrequenz/2) eine komplexe Zahl, die Amplitude und Phase der entsprechenden Sinusschwingung repräsentiert.

Grafisch lässt sich die FFT so interpretieren, dass ein orthogonales Koordinatensystem bestehend aus Sin/Cos auf das Signal angewendet wird, um dessen Vektoranteile in diesem Koordinatensystem zu bestimmen, die dann zu Betrag (Länge) und Phase (Winkel) umgerechnet werden.

Eine FFT arbeitet daher grundsätzlich mit komplexen Signalen, das heißt man gibt ein Signal mit einem Real- und Imaginärteil herein, und heraus kommt ein Spektrum mit einem Real- und Imaginärteil. Arbeitet man mit rein reellen Signalen (der Normalfall), dann setzt man den Imaginärteil des Eingangssignals einfach auf 0. Ansonsten gibt es auch spezielle FFT-Implementierungen für rein reelle Signale, und auch ein paar Tricks um die Transformation von rein reellen Signalen mit komplexen FFTs beschleunigen kann. Darauf soll hier aber erst mal nicht eingegangen werden.

Leistungsdichte

Oft interessiert man sich nicht für die Phasen der Sinuskomponenten des Signals, sondern nur für die Beträge, bzw. nach Quadrierung die Leistung. Wenn man die Beträge der komplexen Zahlen quadriert erhält man das sogenannte Leistungsdichtespektrum des Signals, kurz PSD (Power Spectral Density).

(Plot Zeitsignal -> PSD)

Fensterung

(Plot 1. abgeschnittener Sinus, 2. mit Fensterung)

Realität: stochastische Signale

In realen Anwendungen hat man meist kein Signal mit einer definierten Form wie einen Sinus, sondern Signale von denen man nicht weiss wie sie aussehen, die sich ständig ändern, und von denen man nur Erwartungswerte kennt, z. B. den Mittelwert oder die Leistung. Solche Signale nennt man Zufallssignale oder stochastische Signale. Wenn wir ein Stück aus so einem Signal, hier als Beispiel ein Sprachsignal, ausschneiden, fenstern, FFT-transformieren, den Betrag bilden und quadrieren, dann sieht das Ergebnis zum Beispiel so aus:

(Plot gefenstertes Sprachsignal, Spektrum)

Mittelung: Die Welch-Methode

Diese Varianz zu reduzieren geht glücklicherweise sehr einfach: man bildet einen Mittelwert über die FFTs von mehreren Blöcken. Wie viele Blöcke man verwendet bestimmt die Reaktionszeit des Spektrums auf Änderungen. Ein weiterer Faktor ist die Überlappung zwischen den Blöcken, hierfür wird meistens 50% gewählt.

(Plot gemitteltes Spektrum)

Mittelung: Rekursiv

Ein Nachteil bei der Welch-Methode ist der hohe Speicherbedarf, da die Ergebnisse aller Blöcke über die gemittelt werden soll ja irgendwo gespeichert werden müssen. Alternativ kann man daher eine rekursive Mittelung verwenden.

[math]\displaystyle{ S_{neu} = S_{alt} \cdot \lambda + \mathrm{FFT}(aktueller Block) \cdot (1 - \lambda) }[/math]

Beispiel:

[math]\displaystyle{ S_{neu} = S_{alt} \cdot 0.99 + \mathrm{FFT}(aktueller Block) \cdot 0.01 }[/math]

Jetzt ist neben dem Speicher für das FFT-Ergebnis und das Spektrum kein zusätzlicher Speicher mehr nötig.

Auch hier stellt sich wieder die Frage welche Fensterung verwendet wird, und ob aufeinanderfolgende Datenblöcke überlappen sollen oder nicht. Hierfür gelten die selben Überlegungen wie für die Welch-Methode.

Außerdem gilt es natürlich den Parameter [math]\displaystyle{ \lambda }[/math] zu wählen. In der Praxis ist i.d.R. Ausprobieren angesagt. Wenn man schon eine Mittelungslänge für die Welch-Methode gefunden hat, dann lässt sich mit folgender Faustregel ein äquivalenter Lambda-Faktor für die rekursive Mittelung berechnen:

[math]\displaystyle{ R = (1 + \lambda)/(1 - \lambda) }[/math]

(Plot gefiltertes Spektrum)

Effiziente Implementierung

Da Multiplikationen auf vielen Mikrocontrollern aufwändig sind, ist es sinnvoll sich auf einen Faktor [math]\displaystyle{ \lambda = 1 - 2^{-x} }[/math] zu beschränken und die Berechnung mit Schiebeoperationen durchzuführen:

[math]\displaystyle{ S_{neu} = S_{alt} \cdot (1 - 2^{-x}) + \mathrm{FFT}(aktueller Block) \cdot 2^{-x} }[/math]

Mit ein bisschen Umstellen und Ersetzen der Multiplikation mit [math]\displaystyle{ 2^{-x} }[/math] durch x-faches Rechtsschieben kommt man auf:

[math]\displaystyle{ Differenz = \mathrm{FFT}(aktueller Block) - S_{alt} }[/math]

[math]\displaystyle{ Aenderung = Differenz \gt \gt x }[/math]

[math]\displaystyle{ S_{neu} = S_{alt} + Aenderung }[/math]

Ein Nachteil der rekursiven Mittelung ist, dass Lambda-Werte sehr nahe bei 1 (bzw. ein großer Wert x) dazu führen können, dass "Aenderung" durch die begrenzte Genauigkeit immer 0 ist, und sich [math]\displaystyle{ S_{neu} }[/math] deshalb nicht ändert.