Python: Fast Fourier Transform (FFT)
Fast and Fourier, Teil 1
Man stolpert als Nachrichtentechniker immer wieder über die FFT. Bei mir geschah es vor ein paar Abenden, als ich mir LoRa genauer angeschaut habe. Während der Demodulation eines LoRa-Signals, wird dieses zunächst mit einem sog. downchirp multipliziert. Als Ergebnis dieser Operation bleibt Signal mit konstanter Frequenz übrig. Mit dieser resultierenden Frequenz ist nun das eigentliche Symbol codiert. Wie misst bzw. decodiert man diese Frequenz? Mit Hilfe einer Fast Fourier Transformation.
Ich wusste einmal, wie man eine FFT programmiert. Das ist gut 30 Jahre her. Mein persönlicher Ehrgeiz verbietet mir die KI zu befragen. Also leite ich mir den Algorithmus noch ein letztes Mal selber her, schreibe eine eigene Software dazu und vergleiche das Ergbnis mit dem einer einer frei verfügbaren Library. Zuletzt schauen wir noch, welcher Code der KI dazu einfällt...
Man startet zunächst mit der Gleichung der Diskreten Fourier-Transformation (DFT), die sich samt Skalierungsfaktor 2/N auch leicht aus dem Fourier-Integral herleiten lässt.
Bei den weiteren Betrachtungen wird auch der Skalierungsfaktor 2/N zunächst vernachlässigt.
Bei einer DFT über N Stützpunkte sind N komplexwertige Multi-plikationen N-Mal auszuführen.
Der Ansatz ist nun die DFT zunächst in zwei DFT aufzutrennen: eine für die geraden Stützpunkte x(2n) und eine für die ungeraden Stützpunkte x(2n+1) im Zeitbereich.
Man zerlegt also die ursprüngliche DFT nach dem Prinzip
X(k)=Xg(k) + wN(k) * Xu(k) in zwei "kleinere" DFT, für die der Rechenaufwand nur noch 0,5 * N^2 + N beträgt. Die beiden Teil-DFT Xg(k) und Xu(k) lassen sich ihrerseits wiederum zerlegen, solange bis N/2 DFT über lediglich zwei Stützstellen übrigbleiben.
Die rekursive Umsetzung des FFT-Algorithmus ist nun offensichtlich. Bevor es jedoch ans Programmieren geht, sollte man sich unbedingt noch weitere Aspekte der FFT anschauen. Vielleicht lässt sich ja noch mehr Rechenzeit einsparen.
Da sind beispielsweise die sog. Twiggle-Faktoren wN(k), die nicht unbedingt bei jedem Durchlauf k neu berechnet werden sollten. Die komplexwertigen wN(k) werden vorab berechnet und lediglich aufgerufen.
Da N = 2^p, kann man auch auf zahlreiche Symmetrien komplexer Zeiger zurückgreifen. Als Beispiel sein hier genannt wN(N) = +1, das der Computer gerne als (0.999....-j0.123E-15) berechnet. Was für eine Verschwendung!
Nicht ganz so offensichtlich und dennoch nützlich ist die Eigenschaft, dass wN(-nk) = (wN(nk))*, also konjugiert komplex. Wer es selber nachrechnen möchte, verwende bitte die nebenstehenden Formeln und folgenden Ansatz: wN( k[N-n] ).
Bitte noch nicht mit dem Programmieren beginnen. Erst Denken, dann Programmieren!
Welchen Laufzeitvorteil der FFT gegenüber einer DFT kann man erwarten?
N lässt sich dazu als Potenz zur Basis 2 formulieren. Schreibt man nun den Aufwand für p = 1, 2, 3 auf, so lässt sich bald eine Abschätzung für beliebige p ableiten.
Die O-Notation war zwar noch nie meine Stärke, aber das Ergebnis ist an dieser Stelle zumindest gleich zu den Literaturwerten (und zu dem, was die KI nachplappert).
Endlich programmieren!
Als Referenz lasse ich eine DFT() über N Stützpunkte mitlaufen. Danach erarbeite ich eine einfache FFT() nach dem oben vorgestellten Muster. Weiters nutze ich die hier vorgestellten Optimierungspotenziale in einer advFFT() aus. Den Rechenzeitbedarf dieser Funktionen vergleiche ich mit der durch numpy bereit-gestellten FFT.
Die DFT() wird bereits ab 4096 Stützstellen sehr anstrengend. Die FFT() ist natürlich überlegen und belegt zudem die Zeitkomplexität von NlogN. Dennoch fördert die advFFT() weitere Rechenzeitvorteile, insbesondere durch die Ausnutzung der konjugiert komplexen Symmetrien der Twiggle-Faktoren. Dadurch wird die Dekomposition der Rekursion beschleunigt.
Bei bis zu 1048 Stützstellen ist meine advFFT() der von numpy gelieferten FFT stets überlegen. Doch schon ab 2048 Stützstellen gewinnt numpy immer. Den FFT-Code von numpy werde ich mir in einem nächsten Schritt anschauen müssen. Mal schauen, was da für ein Voodoo drin steckt...
Zuletzt frage ich die KI nach einem Python-Algorithmus und werde von ihrem ersten Entwurf enttäuscht. Dieser Code ist sogar noch ein paar Millisekunden langsamer als meine ursprüngliche FFT()! Ich habe jetzt auch keine Lust, der KI das beizubringen. Soll sie doch selber schauen, wo sie bleibt...
Wie könnte es von hier aus weitergehen? Man könnte noch einen Versuch mit einem 8-bit Arduino machen! Der FFT-Algorithmus sollte dazu in Integerarithmetik konvertiert werden. Ob damit die Echtzeit-FFT von Sprache möglich wäre?
Den Code zum Blog gibt es hier: