~ XV: Die schnelle Fourier-Transformation

John W. Tukey
John W. Tukey (1915-2000), Quelle: Wikimedia

Die schnelle Fourier-Transformation ist die Grundlage dafür, dass, neben der Eleganz der Formeln selbst, sich die Fourier-Transformation am Computer so häufig wiederfindet und von einigen sogar als einer der wichtigsten Algorithmen des 20. Jahrhunderts angesehen wird. Wie der Name schon verrät, ist die schnelle Fourier-Transformation eine Berechnung der diskreten Fourier-Koeffizienten, die eine schnellere Berechnung ermöglicht, als die Betrachtungen aus dem letzten Beitrag. Ihre Bekanntheit geht zurück auf einen Artikel von James W. Cooley (*1926) und John W. Tuches (1915-2000) aus dem Jahre 1965 CT65. Die Ideen zu diesem Algorithmus reichen aber viel weiter zurück.

Wie in den vorigen Beitrag ist das Ziel, die Summe \[ \hat y_j = \sum_{k=0}^{N-1} y_k \text{e}^{-2\pi\textrm{i}\frac{kj}{N}},\quad j=0,\ldots,N-1, \] für gegebene Werte \(y_j\), beispielsweise Abtastwerte einer Funktion, zu berechnenWenn man es genau nimmt, ist diese Formel eine Mischform: Die diskrete Fourier-Transformation besitzt noch den Vorfaktor \(\frac{1}{N}\), die inverse diskrete anstelle des \(-2\pi\) im Exponenten ein \(2\pi\), aber im Prinzip deckt diese Formel beide Fälle ab.. Dabei konzentriert sich dieser Beitrag auf den Fall, dass \(N\) eine Zweierpotenz ist, also als \(N = 2^n\), \(n\in\mathbb N\), geschrieben werden kann. Beispiele sind \(N=2,4,8,16,32,64\) zu den Exponenten \(n=1,2,3,4,5,6\). Nützlich ist hier, dass der Logarithmus (zur Basis 2) in diesem Fall direkt angegeben werden kann als \(\log N = n\).

Die Grundlegende Idee der schnellen Fourier-Transformation ist das „Teile und Herrsche“, einem Prinzip für Algorithmen, bei dem das eigentliche Problem in kleinere Probleme zerlegt wird, die sich effizient lösen lassen und ebenso effizient zur Gesamtlösung zusammengefasst werden können. Hier ist das Teilen ein Zerlegen der obigen Summe in zwei Teile, eine Summe, die alle geraden \(j\) (inklusive der \(j=0\)) enthält und eine, die \(y_j\) mit ungeradem Index als Summanden enthält. Die Formel lautet dann \[ \hat y_j = \sum_{k=0}^{N/2-1} y_{2k} \text{e}^{-2\pi\textrm{i}\frac{2kj}{N}} + \sum_{k=0}^{N/2-1} y_{2k+1} \text{e}^{-2\pi\textrm{i}\frac{(2k+1)j}{N}} ,\quad j=0,\ldots,N-1. \] Da das etwas monströs daherkommt, sei als Beispiel \(N=8\) betrachtet. Für ein \(j\in\{0,1,\ldots,7\}\) wird die gesamte Summe aufgeteilt in \[ \begin{align*} \hat y_j &= y_0 \text{e}^{-2\pi\textrm{i}\frac{0}{16}} + y_1 \text{e}^{-2\pi\textrm{i}\frac{1j}{16}} + y_2 \text{e}^{-2\pi\textrm{i}\frac{2j}{16}} + y_3 \text{e}^{-2\pi\textrm{i}\frac{3j}{16}} + y_4 \text{e}^{-2\pi\textrm{i}\frac{4j}{16}} + y_5 \text{e}^{-2\pi\textrm{i}\frac{5j}{16}} + y_6 \text{e}^{-2\pi\textrm{i}\frac{6j}{16}} + y_7 \text{e}^{-2\pi\textrm{i}\frac{7j}{16}}\\ &=\bigl( y_0 \text{e}^{-2\pi\textrm{i}\frac{0}{16}} + y_2 \text{e}^{-2\pi\textrm{i}\frac{2j}{16}} + y_4 \text{e}^{-2\pi\textrm{i}\frac{4j}{16}} + y_6 \text{e}^{-2\pi\textrm{i}\frac{6j}{16}} \bigr)\\&\qquad +\bigl( y_1 \text{e}^{-2\pi\textrm{i}\frac{1j}{16}} + y_3 \text{e}^{-2\pi\textrm{i}\frac{3j}{16}} + y_5 \text{e}^{-2\pi\textrm{i}\frac{5j}{16}} + y_7 \text{e}^{-2\pi\textrm{i}\frac{7j}{16}} \bigr). \end{align*} \] Um auf die obigen Darstellungen der zwei Summen zu kommen, muss man sich für das \(k\), die Summenlaufvariable, zwei Dinge überlegen: In welcher Schrittweite läuft der neue Exponent (also welchen Vorvaktor benötigt das \(k\)) und wo hört das \(k\) auf. Für die erste Summe der Zerlegung gilt: Die Summe beginnt mit \(y_0\), \(y_2\), und läuft somit in Zweierschritten, daher der Wert \(2k\) in den Indizes der \(y_{2k}\). Sie endet schließlich mit \(y_{6}\), was sich in diesem Beispiel schreiben lässt als \(6=8-2 = N-2\). Klammert man dort eine zwei aus (denn für den letzten Summanden gilt \(6=2k\) für ein bestimmtes \(k\), genauer \(k=3\)), so ergibt sich als letzter Index \(\frac{N}{2}-1 = \frac{8}{2}-1 = 3\). Der Nachteil der Summenschreibweise ist also, dass dieses \(k\) so genau bestimmt werden muss. Der Vorteil ist allerdings, dass sie sehr platzsparend ist und man schnell die allgemeine Regel in den Summanden wiedererkennt.

Nun folgt der eine kleine Schritt, den man gar nicht so einfach sieht. Der Faktor \(\text{e}^{-2\pi\textrm{i}\frac{(2k+1)j}{N}}\), der in jedem zweiten Summanden vorkommt, lässt sich aufspalten in \(\text{e}^{-2\pi\textrm{i}\frac{(2k+1)j}{N}} = \text{e}^{-2\pi\textrm{i}\frac{2kj}{N}}\cdot\text{e}^{-2\pi\textrm{i}\frac{j}{N}}\). Der zweite Faktor kommt jedoch in jedem Summanden vor und kann somit ausgeklammert, sprich vor die Summe gezogen, werden. Die Fourier-Transformation lässt sich also umgeformt auch schreiben als \[ \hat y_j = \sum_{k=0}^{N/2-1} y_{2k} \text{e}^{-2\pi\textrm{i}\frac{2kj}{N}} + \text{e}^{-2\pi\textrm{i}\frac{j}{N}}\sum_{k=0}^{N/2-1} y_{2k+1} \text{e}^{-2\pi\textrm{i}\frac{2kj}{N}} ,\quad j=0,\ldots,N-1. \]

Schließlich lässt sich in beiden Summen der Bruch im Exponenten noch etwas anders schreiben \(\frac{2kj}{N} = \frac{kj}{N/2}\). Mit diesen Umformungen ist die erste Summe der \(j\)te diskrete Fourier-Koeffizient der diskreten Fourier-Transformation mit Länge \(N/2\) und den Eingabewerten \(y_0,y_2,\ldots,y_{N}\), denn dies sind genau \(N/2\) Werte. Für die zweite Summe muss man etwas genauer hinschauen: Die Summe enthält ebenso \(N/2\) Summanden, benötigt die Werte \(y_1,y_3,\ldots,y_{N-1}\). Es handelt sich also auch um den \(j\)ten Wert der dazugehörigen Fourier-Transformation der Länge \(N/2\), nur dass danach das Ergebnis der Summe noch mit dem Vorfaktor multipliziert werden muss. Falls \(j>N/2\) ist, so lässt sich die Eigenschaft verwenden, dass die Fourier-Koeffizienten periodisch sind und sich somit nach einer Periode wiederholen. Dies liegt an der Periodizität der Exponentialfunktion \(\textrm{e}^{\textrm{i}x}\), siehe Beitrag IV. Daher ist dann der Koeffizient mit \(j’=j-N/2\) derjenige, der in beiden Summen von Interesse ist. Mit anderen Worten gilt von der andern Seite betrachtet: Der jeweils \(j\)te Koeffizient aus einer der beiden Fourier-Transformationen halber Länge (also \(j\in\{0,1,\ldots,N/2-1\}\) wird sowohl für \(\hat{y}_j\) als auch für \(\hat{y}_{j+N/2}\) verwendet. Der Vorfaktor, der vor der zweiten Summe (also dem Fourier-Koeffizienten zu den ungeraden \(y_j\)) ist jedoch für diese beiden Koeffizienten verschieden.

In dieser Zerlegung steckt schon die wesentliche Idee der schnellen Fourier-Transformation. Aber wieso sollte es schneller gehen, wenn wir statt einer Fourier-Transformation mit \(N\) Daten zwei kleine lösen? Einen ersten Hinweis dazu gibt die Tabelle aus dem letzten Beitrag. Für \(N=8\) stehen dort schon \(120\) Rechenoperationen, die zur Berechnung notwendig sind. Die gerade geführten Überlegungen erlauben jedoch, stattdessen zweimal die Fourier-Transformation \(N=4\) zu nutzen, dann eine Multiplikation (der Faktor vor der zweiten Summe) und eine Addition (der beiden Transformationen), um das Ergebnis zu erhalten. Nimmt man dazu vorerst ebendiese Tabelle wieder zu Hilfe, so sind dies \(2\cdot28 + 2 = 58\) Rechenschritte. Schon für so eine kleine Zahl reduziert sich mit nur einer Zerlegung der Summe der Aufwand also enorm.

Doch diese Idee lässt sich noch weiter führen, indem die beiden Fourier-Transformationen, die aus der Zerlegung entsteht, weiter zerlegt werden. Das illustriert die Abbildung 1: Die \(N=8\) Eingabewerte werden aufgeteilt in gerade und ungerade Koeffizienten, um die entstehenden zwei Blöcke wieder genauso aufzuteilen. Die dritte Zerlegung führt schließlich zu einzelnen Zahlen, also Blöcken der Länge \(N=1\).

Darstellung der schnellen Fourier-Transformation mittels Schmetterlingsgraphen: Die Eingabe (1. Zeile) wird in Zeilen 2-4 zerlegt. Danach in Zeilen 5-7 wird jeweils das Ergebnis aus der Fourier-Transformation auf den ungeraden Indizes (jeweils innerhalb eines Block betrachtet) mit einem twiddle factor multipliziert (Halbkreise) und diejenige der geraden Indizes addiert (Pfeile). Die unterschiedlichen Faktoren werden beim Überfahren angezeigt. Zählen aller Pfeile und Kreise ist eine Möglichkeit, die Anzahl Operationen nachzuvollziehen.

Diese unterste Stufe der Zerlegung ist allerdings sehr einfach: Eine Fourier-Transformation der Länge \(N=1\) besteht aus einer Summe, die lediglich einen Koeffizienten \(y_j\) involviert, mehr noch, der entsprechende Faktor dahinter ist (wegen \(k=0\) im ersten Summanden) \(\textrm{e}^0 = 1\). Für diesen Schritt ist keine einzige Operation notwendig. Danach folgt ein zweiter Trick: Statt nur ein \(\hat y_j\) zu berechnen, können alle gleichzeitig in dem Diagramm berechnet werden, denn für jeden geschieht die gleiche Zerlegung, einzig beim Zusammenaddieren ändern sich die Vorfaktoren vor der zweiten Summe, die auch twiddle factors, zu deutsch Drehfaktoren, genannt werden.

Für die Anzahl Rechenschritte ist wichtig, dass diese Zerlegung und Brechnung, wie sie schematisch in Abbildung 1 gezeigt ist, gleichzeitig die Fourier-Koeffizienten \(\hat y_j,\ j=0,\ldots,N-1\) ausrechnet. Auf jeder Stufe, die kleinere Transformationen zusammenfasst, sind höchstens \(2N\) Additionen und ebenso viele Multiplikationen reeller Zahlen (die eigentliche Eingabe wieder getrennt in Real und Imaginärteil) notwendig. Die Anzahl Stufen ergibt sich daraus, wie oft die Eingabelänge halbiert werden kann. Hier ist die Einschränkung vom Anfang wichtig, dass \(N\) eine Zweierpotenz ist, denn so kann \(N\) genau \(n=\log N\) mal halbiert werden. Insgesamt sind \(4*N\log N\) Operationen notwendig, um die gesamte diskrete Fourier-Transformation zu berechnen. Mit einigen Tricks lässt sich dieser Vorfaktor auf \(\frac{34}{9} \approx 3,7777\) reduzieren. Dieses Resultat ist noch gar nicht so alt JF07. Die Anzahl Rechenschritte ist für die gleichen Eingabelängen des letzten Beitrages in Tabelle 1 aufgetragen für den Vorfaktor \(4\).

N248163264128
Rechenschritte der FFT832962566401 5363 584
Rechenschritte der DFT12562409924 03216 25665 280
Vergleich der Anzahl Rechenoperationen der FFT und DFT für verschiedene Eingabelängen.

Die Laufzeiten sind zu Beginn nicht viel besser als bei der ersten Umsetzung als Algorithmus, aber mit zunehmener Größe ist der Einfluss des Faktor \(\log N\) kleiner und kleiner. So nimmt insbesondere die Laufzeit beim Vergrößern von 2 auf 128 um den Faktor 448 zu, was bei der langsamen Implementierung im Vergleich ein Faktor 5 440 war. In der Laundau-Symbolik hat die schnelle Fourier-Transformation insgesamt also eine Komplexität von \[\mathcal O(N\log N),\] sie wächst also ein klein wenig mehr als linear.

Schreibe einen Kommentar

Deine E-Mail-Adresse wird nicht veröffentlicht. Erforderliche Felder sind mit * markiert.