Schnelle Fourier Transformation (FFT)

Die schnelle Fourier Transformation (engl. Fast Fourier Transform, FFT) ist in der Informatik (besonders in der Signalverarbeitung) ein sehr wichtiger Algorithmus zur schnellen Berechnung der diskreten Fourier Transformation (DFT). Um die FFT zu erklären muss vorher noch ein Wenig über die DFT gesagt werden. Ich setze in diesem Artikel außerdem grundlegende Kenntnisse der komplexen Zahlen $\mathbb{C}$ voraus.

Diskrete Fourier Transformation

Die diskrete Fourier Transformation ist sehr schnell erklärt (wenn man sie nicht von der allgemeinen Fourier Transformation herleitet). Sie ist eine Abbildung $\hat{f}: \mathbb{C}^n \to \mathbb{C}^n$. Sie bildet also komplexe Vektoren auf komplexe Vektoren ab. Wir definieren uns zunächst

\begin{align*} \omega &:= e^{\frac{2 i \pi}{n}} \\ &= \cos \frac{2\pi}{n} + i \sin \frac{2\pi}{n} \qquad \text{Nach der Eulerschen Identität} \end{align*}

Man nennt $\omega$ auch die primitive n'te Einheitswurzel. Zur Erinnerung: $n$ ist die Länge unseres "Eingabevektors". OK, jetzt definieren wir uns noch eine Matrix $F$. Dies geschieht ganz einfach:

\begin{gather*} F_{kl} := w^{k \cdot l} \qquad k,l = 0 \ldots n-1 \end{gather*}

Wir nennen $F$ die Fouriermatrix. Anders gesprochen ist dies gerade die Vandermonde-Matrix $V(\omega^0, \omega^1, \ldots, \omega^{n-1})$, aber das ist erstmal egal.

Schließlich definieren wir uns die diskrete Fourier Transformation. Haben wir einen Vektor $x \in \mathbb{C}^n$, so ist die Fourier Transformation $\hat{f}: \mathbb{C}^n \to \mathbb{C}^n$ definiert durch

$\displaystyle\hat{f}(x) := F \cdot x$

Sie ist also nichts weiter als eine Matrix-Vektor-Multiplikation der Fouriermatrix mit dem Vektor $x$. In der Literatur findet man auch oft die Formel

\begin{align*} \hat{f}(x)_p &= \sum_{k=0}^{n-1} e^{\frac{2 i \pi}{n} \cdot k \cdot p} \cdot x_p \\ &= \sum_{k=0}^{n-1} \omega^{kp} \cdot x_p \\ &= \sum_{k=0}^{n-1} F_{kp} \cdot x_p \end{align*}

Man sieht, dies ist nichts anderes als die ausgeschriebene Matrix-Vektor-Multiplikation.

Beispiel

Wie sieht sowas aus? Sagen wir, wir wollen den Vektor $x = \left(\begin{smallmatrix}1\\2\\3\\4\end{smallmatrix}\right)$ transformieren. Es ist also $n = 4$. Die 4'te primitive Einheitswurzel ist $\omega = e^{\frac{2 i \pi}{4}} = i$. Konstruktion der Fouriermatrix ergibt

\begin{gather*} F = \begin{pmatrix}1&1&1&1\\1&i&-1&-i\\1&-1&1&-1\\1&-i&-1&i\end{pmatrix} \end{gather*}

Die Rechnung $F \cdot x$ ist nun wirklich trivial und es ist $\hat{f}(x) = \left(\begin{smallmatrix}10\\-2 - 2i\\-2\\-2+2i\end{smallmatrix}\right)$. Dies ist schon unser transformierter Vektor!

Inverse Diskrete Fourier Transformation

Jetzt stellt sich die Frage wie wir den Vektor $\hat{f}$ wieder zurücktransformiert bekommen. Klar dürfte sein, dass eine Multiplikation mit der inversen Fouriermatrix $F^{-1}$ wieder den ursprünglichen Vektor liefern wird. Aufgrund der Eigenschaften der primitiven Einheitswurzeln und der Symmetrie von $F$ ist die Inverse Matrix $F^{-1}$ gegeben durch

\begin{gather*} F^{-1}_{kl} = \frac{1}{n} \cdot \overline{F_{kl}} = \frac{1}{n} \omega^{-k \cdot l} \end{gather*}

Dabei bezeichnet $\overline{F_{kl}}$ die konjugiert komplexe Zahl von $F_{kl}$. Die Rücktransformation ist also nun auch wieder nur eine Matrix-Vektor-Multiplikation, nämlich

$\displaystyle\hat{f}^{-1}(\hat{x}) = F^{-1} \cdot \hat{x}$

That's it! In unserem Beispiel von oben wechseln wir in der Matrix $F$ einfach das Vorzeichen vor dem Imaginärteil jeder komplexen Zahl, und ziehen $\frac{1}{n} = \frac{1}{4}$ vor die Matrix und erhalten

\begin{gather*} F^{-1} = \frac{1}{4} \cdot \begin{pmatrix}1&1&1&1\\1&-i&-1&i\\1&-1&1&-1\\1&i&-1&-i\end{pmatrix} \end{gather*}

womit dann $F^{-1} \cdot \left(\begin{smallmatrix}10\\-2 - 2i\\-2\\-2+2i\end{smallmatrix}\right) = \left(\begin{smallmatrix}1\\2\\3\\4\end{smallmatrix}\right)$ folgt.

Schnelle Fourier Transformation

Die Schnelle Fourier Transformation ist ein Divide & Conquer Algorithmus zur Berechnung der DFT. Im O-Kalkül ergibt sich, dass die DFT wegen der Matrix-Vektor-Multiplikationen einen Aufwand von $\mathcal{O}(n^2)$ hat. Der Algorithmus der FFT schafft die gleichen Berechnungen in $\mathcal{O}(n \log n)$, was eine wesentliche Verbesserung darstellt! Kleiner Haken an der Sache ist jedoch, dass die Länge der zu transformierenden Vektoren eine Zweierpotenz sein muss. In den meisten Fällen ist das aber egal, da man einfach den Vektor mit Nullen auffüllen kann bis seine Länge gerade die einer Zweierpotenz entspricht.

Das Verfahren

Wie schon angesprochen handelt es sich bei der FFT um einen Divide & Conquer Algorithmus, welcher sich am einfachsten rekursiv formulieren lässt. Wir möchten also den transformierten Vektor $b = F \cdot a$ berechnen. Dabei sind $a = (a_1, a_2, \ldots, a_{2^m})$ und $b = (b_1, b_2, \ldots, b_{2^m})$ mit $m \in \mathbb{N}$. Wir können den Vektor $b$ also aufspalten, und zwar in die Vektoren $b^{(g)}$, der gerade die Komponenten von $b$ mit geradem Index und $b^{(u)}$ der die Komponenten von $b$ mit ungeradem Index enthält. Ausgeschrieben heißt das für $b^{(g)}$

\begin{align*} b_k^{(g)} &= b_{2k}\\ &= \sum_{l=0}^{n-1} \omega^{l(2k)} a_l \\ &= \sum_{l=0}^{\frac{n}{2} - 1} w^{l(2k)} a_l + \sum_{l=0}^{\frac{n}{2}-1} \omega^{(\frac{n}{2} + l) 2k} a_{\frac{n}{2} + l} \\ &= \sum_{l=0}^{\frac{n}{2} - 1} w^{l(2k)} (a_l + a_{\frac{n}{2} + l}) \end{align*}

Uff, was wurde hier gemacht? Zunächst einmal berechnen sich die geraden Komponenten von $b$ ja dadurch dass wir die geraden Zeilen der Matrix $F$ mit $a$ multiplizieren. Dies lässt sich als Summe schreiben. Die Summe können wir in der Mitte teilen und erhalten damit die dritte Zeile. Jetzt nutzen wir aus, dass $\omega$ eine primitive Einheitswurzel ist. Dann gilt nämlich

\begin{align*} \omega^{(\frac{n}{2} + l) 2k} &= \omega^{nk + l(2k)} \\ &= \omega^{l(2k)} \cdot \omega^{nk} \\ &= \omega^{l(2k)} \cdot e^{i \frac{2\pi nk}{n}} \\ &= \omega^{l(2k)} \cdot e^{i (2\pi k)} \\ &= \omega^{l(2k)} \cdot ( \underbrace{\cos (2k \pi)}_{= 1} + i \underbrace{\sin (2k \pi)}_{= 0})\\ &= \omega^{l(2k)} \end{align*}

Damit kann man also die Summe zusammenziehen (wie oben in der letzten Zeile gezeigt)! Substituiere nun $m = \frac{n}{2}$ und $\nu = \omega^2$ ($\nu$ ist also m-te primitive Einheitswurzel) und wir erhalten

\begin{gather*} b_k^{(g)} = \sum_{l=0}^{m - 1} \nu^{lk} (a_l + a_{m+l}) \end{gather*}

Das ist also nichts anderes als die Fouriertransformation des Vektors $\left(\begin{smallmatrix}a_0 + a_{m} \\ \vdots \\ a_{m-1} + a_{n-1} \end{smallmatrix}\right)$. Der Vektor hat gerade die Länge $\frac{n}{2} = m$.

Machen wir das gleiche Spiel mit den ungeraden Indizes von $b$:

\begin{align*} b_k^{(u)} &= b_{2k+1} \\ &= \sum_{l=0}^{n-1} \omega^{l(2k+1)} a_l \\ &= \sum_{l=0}^{\frac{n}{2} - 1} w^{l(2k+1)} a_l + \sum_{l=0}^{\frac{n}{2}-1} \omega^{(\frac{n}{2} + l)(2k+1)} a_{\frac{n}{2}+l} \\ &= \sum_{l=0}^{\frac{n}{2} - 1} \omega^l w^{l(2k)} a_l + \sum_{l=0}^{\frac{n}{2}-1} - \omega^l \omega^{l(2k)} a_{\frac{n}{2}+l} \\ &= \sum_{l=0}^{n-1} \omega^{l(2k)} \omega^l (a_l - a_{\frac{n}{2} + l}) \end{align*}

Die Umformungen von Zeile 3 nach 4 in der rechten Summe ausführlich:

\begin{align*} \omega^{(\frac{n}{2} + l)(2k+1)} a_{\frac{n}{2}+l} &= \omega^{nk + \frac{n}{2} + l(2k) + l}\\ &= \underbrace{\omega^{nk}}_{= 1} \omega^{\frac{n}{2}} \omega^{l(2k)}\omega^l\\ &= e^{i \pi} \omega^{l(2k)} \omega^l \\ &= (\underbrace{cos \pi}_{= -1} + i \underbrace{\sin \pi}_{=0}) \omega^{l(2k)} \omega^{l} \\ &= - \omega^{l(2k)} \omega^{l} \end{align*}

Ok, wir können also wieder $m = \frac{n}{2}$ und $\nu = \omega^2$ substituieren und erhalten

\begin{gather*} b_k^{(u)} = \sum_{l=0}^{m-1} \nu^{lk} \omega^l (a_l - a_{m+l}) \end{gather*}

Wieder ist dies nichts anderes als die Fouriertransformation des Vektors $\left(\begin{smallmatrix}\omega^0(a_0 - a_{m}) \\ \vdots \\ \omega^{m-1}(a_{m-1} - a_{n-1}) \end{smallmatrix}\right)$.

Diese Zerlegung lässt sich jetzt rekursiv weiter durchführen bis wir Vektoren der Länge 1 haben. Die Fouriertransformierte eines Vektors der Länge 1 ist jedoch der Vektor selbst (Es gilt ja bei $n = 1$ dass $F = (1)$ ist), somit haben wir auch eine rekursive "Abbruchbedingung".

Zusammengefasst geht man also so vor:

  1. Berechne temporäre Vektoren $a' = \left(\begin{smallmatrix}a_0 + a_{m} \\ \vdots \\ a_{m-1} + a_{n-1} \end{smallmatrix}\right)$ und $ a'' = \left(\begin{smallmatrix}\omega^0(a_0 - a_{m}) \\ \vdots \\ \omega^{m-1}(a_{m-1} - a_{n-1}) \end{smallmatrix}\right)$ aus $a$. Diese haben jeweils nurnoch halbe Länge!
  2. Berechne rekursiv die Fouriertransformierten der Vektoren $a'$ und $a''$. Sie seien mit $\hat{a}'$ und $\hat{a}''$ bezeichnet. Beachte dass wegen der halben Länge der Vektoren die $\frac{n}{2} = \nu$-te primitive Einheitswurzel benutzt werden muss!
  3. Füge $\hat{a}'$ und $\hat{a}''$ reißverschlussartig zu $\hat{a}$ zusammen, so dass die geraden Indizes von $\hat{a}$ gerade $\hat{a}'$ und die ungeraden Indizes von $\hat{a}$ gerade $\hat{a}''$ entsprechen.

Dieses Verfahren hat nun eine Komplexität von nurnoch $\mathcal{O}(n \log n)$!

Inverse Schnelle Fourier Transformation

Die Umkehrung der FFT ist wieder eine FFT. Wie wir jedoch im Kapitel der DFT gesehen haben muss diesmal die Matrix $F^{-1}$ benutzt werden. Wir benutzen also die konjugiert komplexen Einheitswurzeln $\overline{\omega}$ statt $\omega$ und müssen die Endergebnisse noch mit $\frac{1}{n}$ multiplizieren. Aufwand: Wieder $\mathcal{O}(n \log n)$.

Das Butterfly Schema

Das oben angegebene Verfahren lässt sich algorithmisch sehr schön implementieren, aber will man die FFT auf dem Papier durchführen ist es etwas umständlich, da man viele temporäre Vektoren ausrechnen muss. Abhilfe schafft da das hochbedenkliche Butterfly Schema.

Als Butterfly bezeichnet man einen Graph bei dem jeder Knoten einen numerischen Wert enthält der sich als gewichtete Summe der eingehenden Kanten zusammensetzt. Beispiel:

Standard Butterfly Graph

Beim Butterfly Schema fängt man so an, dass man sich erst den Vektor jeweils in ungerade und gerade Komponenten aufspaltet bis man lauter Vektoren der Länge 1 erhält. Diese schreibt man dann untereinander. Nehmen wir zum Beispiel den Vektor

\begin{gather*} a = \begin{pmatrix} 0\\1\\2\\3 \end{pmatrix} \end{gather*}

Aufspalten jeweils in gerade und ungerade Teilvektoren ergibt diesen Baum

Beispielbaum

Wir schreiben nun einfach die Blätter des Baums untereinander und erhalten somit den Anfang.

Anfang des Schemas

Im $l$'ten Schritt des Schemas werden nun die Rechenoperationen wie folgt dargestellt durchgeführt:

Schritt

Wir haben auf der rechten Seite den zusammengesetzten Vektor aus den beiden Vektoren halber Länge die sich auf der linken Seite befinden. Man sieht hier schön die reißverschlussartige Verschränkung. In unserem Beispiel sähe das komplette Schema folgendermaßen aus:

Komplettes Beispielschema

Da unser Vektor die Länge 4 hatte, haben wir folgende Werte für die Potenzen von $\omega$ benutzt

\begin{gather*} (\omega^0, \omega^1, \omega^2, \omega^3) = (1, i, -1, -i) \end{gather*}

Das Ergebnis ist also

\begin{gather*} b = \begin{pmatrix}6\\-2-2i\\-2\\-2+2i\end{pmatrix}\end{gather*}

That's it!

Weiteres zur FFT

Website by: Thomas Pajor, Augartenstr. 56, 76137 Karlsruhe - www.darkviper.de - with help from: flowhase