Fast-Fourier-Transformation und Zeitkomplexität

Das Problem

$$ X(n) = \sum_{k=0}^{N-1}{x[k]e^{-j(2 \pi kn/N)}} $$

Wenn wir noch einmal auf die DFT-Formel schauen, können wir sehen, dass zur Berechnung eines Bandes N Additionen und N Multiplikationen erforderlich sind (wobei N die Größe des Fensters ist). Um N-Bänder zu bekommen, sind 2·N² Operationen notwendig, was sehr viel ist.

Angenommen, wir haben einen dreiminütigen Song mit 44,1 kHz und berechnen das Spektrogramm des Songs mit einem 4096-Sample-Fenster. Wir müssen 10,7 (44100/4096) DFT pro Sekunde berechnen, also 1938 DFTs für das vollständige Lied. Jede DFT benötigt 3,35·10⁷ Operationen (2·4096²). Um das Spektrogramm des Songs zu erhalten, müssen wir 6,5·10¹⁰ Operationen durchführen.

Nehmen wir an, wir haben eine Musiksammlung von 1000 drei Minuten langen Songs, dann benötigen wir 6,5·10¹³ Operationen, um die Spektrogramme Ihrer Songs zu erhalten. Selbst mit einem guten Prozessor würde es Tage bzw. Monate dauern, um das Ergebnis zu erhalten.

Zum Glück gibt es schnellere Implementierungen der DFT namens FFT (Fast Fourier Transform). Einige Implementierungen erfordern nur 1,5·N·log(N) Operationen. Für die gleiche Musiksammlung erfordert die Verwendung der FFT anstelle der DFT 340 mal weniger Additionen (1,43·10¹¹) und es würde nur Minuten bzw. Stunden dauern, um das Ergebnis zu erhalten.

Dieses Beispiel zeigt einen weiteren Kompromiss: Obwohl die Vergrößerung des Fensters die Frequenzauflösung verbessert, erhöht es auch die Rechenzeit. Wenn wir für die gleiche Musiksammlung das Spektrogramm mit einem 512-Sample-Fenster (Frequenzauflösung von 86 Hz) berechnen, erhalten wir das Ergebnis mit der FFT in 1,07·10¹¹ Operationen, etwa 1/4 mal schneller als mit einem 4096-Sample-Fenster (Frequenzauflösung von 10,77 Hz).

Diese Zeitkomplexität ist wichtig, denn wenn wir Shazam für einen Song verwenden, muss unser Smartphone das Spektrogramm des aufgezeichneten Audios berechnen und ein mobiler Prozessor ist (derzeit noch) weniger leistungsfähig als ein Desktop-Prozessor.

Downsampling

Zum Glück gibt es einen Trick, um die Frequenzauflösung beizubehalten und gleichzeitig die Fenstergröße zu reduzieren, er wird "Downsampling" genannt. Nehmen wir einen Standard-Song bei 44100 Hz, wenn wir ihn bei 11025 Hz (44100/4) umrechnen (resampeln), erhalten wir die gleiche Frequenzauflösung, egal ob wir eine FFT auf dem 44.1 kHz-Song mit einem 4096-Fenster oder eine FFT auf dem 11 kHz-Song, der umgerechnet wurde, mit einem 1024-Fenster machen. Der einzige Unterschied besteht darin, dass der neu abgetastete Song nur Frequenzen von 0 bis 5 kHz hat. Aber der wichtigste Teil eines Songs ist zwischen 0 und 5 kHz. Tatsächlich hören die meisten von uns keinen großen Unterschied zwischen einer Musik mit 11 kHz und einer Musik mit 44,1 kHz. Die wichtigsten Frequenzen befinden sich also immer noch im umgerechneten Song, was für einen Algorithmus wie Shazam wichtig ist.

Abbildung 19

Das Downsampling eines 44,1 kHz Songs auf einen 11 025 kHz Song ist nicht sehr schwierig: Ein einfacher Weg dies zu tun ist, die Samples jeweils als 4er-Gruppe zu nehmen und diese Gruppe zu nur einem Sample zu transformieren, indem man den Durchschnitt der 4 Samples bildet. Der einzige knifflige Teil ist, dass man vor dem Downsampling eines Signals die höheren Frequenzen im Sound filtern muss, um Aliasing zu vermeiden (siehe Nyquist-Shannon-Theorem). Dies kann durch Verwendung eines digitalen Tiefpassfilters („low pass filter") erfolgen.

FFT

Gehen wir zurück zur FFT. Die einfachste Implementierung der FFT ist der Radix-2-Cooley-Tukey-Algorithmus, der ein Teile-und-herrsche-Algorithmus ist. Die Idee ist, dass anstatt die Fourier-Transformation direkt auf dem N-Sample-Fenster zu berechnen, der Algorithmus wie folgt verfährt:

• Teile das N-Sample-Fenster in 2 N/2-Sample-Fenster.
• Berechne (rekursiv) die FFT für die 2 N/2-Sample-Fenster.
• Berechne effizient die FFT für die N-Sample-Fenster der 2 vorherigen FFT.

Der letzte Teil kostet nur N Operationen mit einem mathematischen Trick an den Wurzeln der Einheit (die exponentiellen Terme).

Hier ist eine lesbare Version der FFT (in der Programmiersprache Python):

from cmath import *
def fft(x):
        N=len(x)
        if N==1: return x

        even=fft([x[k] for k in range(0,N,2)])
        odd= fft([x[k] for k in range(1,N,2)])

        M=N/2
        l=[ even[k] + exp(-2j*pi*k/N)*odd[k] for k in range(M) ]
        r=[ even[k] - exp(-2j*pi*k/N)*odd[k] for k in range(M) ]

        return l+r

Weitere Informationen zur FFT findet man in diesem Wikipedia-Artikel: Cooley–Tukey FFT algorithm