12.2.2 Varianten der schnellen Fouriertransformation

Die im Abschnitt 12.2.1 beschriebene Grundform der schnellen Fouriertransformation verfährt folgendermaßen:

• Die Eingabe der Länge n wird auf zwei Unterfelder der Länge m = n/2 verteilt, wobei das erste nur mit Summen im Indexabstand m, das zweite mit “gedrehten” Differenzen aufgefüllt wird. Das kann man einen “Verteilungs”-Schritt nennen.

• Die rekursive Anwendung der FFT auf diese Teilfelder liefert zwei Felder der Länge m = n/2, die dann umsortiert werden müssen, um die Einträge mit geraden bzw. ungeraden Indizes im Ergebnisfeld zu bilden (“Sortierschritt”).

Eine andere Darstellungsweise bildet eine FFT der Länge n als Summe zweier FFT’s der Länge m = n/2, die einmal über die Eingaben mit ungeradem und einmal über die Eingaben mit geradem Index laufen. Man macht die Aufteilung gerade/ungerade also nicht in der Ausgabe, sondern in der Eingabe. Das führt zu einer im Prinzip äquivalenten Methode, die aber algorithmisch anders ist. Sie soll hier nicht weiter verfolgt werden.

Stattdessen geht es jetzt um Implementationsfragen. Es ist relativ einfach (siehe das “naive” FFT–Programm auf dem Web), den durch (12.2.27) beschriebenen Verteilungs–Schritt und die in (12.2.25) und (2.12.26) gegebene Sortierung rekursiv zu implementieren, wenn man keine Rücksicht auf die Speichertechnik nimmt. Es bedarf aber einiger Zusatzüberlegungen, um das Ganze ohne irgendwelche Zwischenspeicherungen “in place” auszuführen.

Man kann den rekursiven Prozeß schreiben als

F_{n}^{1} = S_{n}^{1}F_{n/2}^{2}V_{n}^{1},

wobei alle Abbildungen auf Vektoren des CI ^{n} wirken. Dabei ist V_{n}^{1} der Verteilungsschritt, S_{n}^{1} der Sortierschritt und F_{n/2}^{2} die Ausführung zweier FFT der Länge n/2. Wendet man das in naheliegender Weise rekursiv an, so folgt

F_{n}^{1} = S_{n}^{1}S_{n/2}^{2}F_{n/4}^{4}V_{n/2}^{2}V_{n}^{1} = S_{n}^{1}S_{n/2}^{2} … S_{1}^{n} ∘ V_{1}^{n} … V_{n/2}^{2}V_{n}^{1},

und man sieht, daß sich die gesamte FFT als “erst verteilen, dann sortieren” schreiben läßt. An dieser Darstellung kann man ablesen, daß die inverse FFT die Form “erst sortieren, dann verteilen” hat. Weil aber die inverse FFT bis auf Kleinigkeiten (negative statt positive Potenzen der Einheitswurzeln und Abschlußdivision durch n) mit der FFT identisch ist, gibt es auch für die normale FFT eine Formulierung als “erst sortieren, dann verteilen”. Auch diese Variante soll hier nicht verfolgt werden.

Die Grundidee zur Effektivierung der Basisversion aus Abschnitt 12.2.1 besteht nun darin, alle Verteilungsschritte zuerst und nichtrekursiv auszuführen, und die Sortierung erst ganz am Schluß und “en bloc” nachzuholen. Der erste Teil ist klar zu sehen: man hat für n = 2^{p} das gesamte Feld p mal zu durchlaufen, wobei man jeweils Verteilungsschritte der Länge n, n/2, n/4 … 4, 2 macht, von denen jeweils 1, 2, …, n/4, n/2 anfallen. Das läßt sich durch einfache Schleifen gut realisieren, wobei man keinerlei Zusatzspeicher braucht und Datenlokalität voll ausnutzen kann. Man sehe sich dazu die Funktion split in der beschleunigten Programmversion auf dem Netz an.

Problematischer ist das Verständnis der abschließenden Sortierung. Der Schritt S_{n}^{1} besteht in C–Notation darin, die Inverse der Indextransformation

\begin{matrix}k & ↦ & 2k & für 0 ≤ k < m = n/2 \\ k & ↦ & 2(k − m) + 1 & für m ≤ k < n \end{matrix}

auszuführen. Hat k eine Bitdarstellung aB mit führendem Bit a ∈ {0, 1} und einem Rest–Binärwort B, so hat das Ergebnis die Bitdarstellung Ba. Die Indextransformation besteht also darin, das erste Bit zu nehmen, abzuschneiden und nach hinten zu setzen. Der Schritt S_{n/2}^{2} bearbeitet die Zahlen mit Binärdarstellungen 0B und 1B unabhängig voneinander und in genau gleicher Weise, wobei er das mit B allein macht, was der vorige Schritt mit der gesamten Binärdarstellung machte. Beide Schritte zusammen transformieren also ein Bitmuster abC wie

abC ↦ aCb ↦ Cba,

weil erst S_{n/2}^{2} das Bit b auf die andere Seite von C stellt, woraufhin S_{n}^{1} noch a umstellt. Man sieht jetzt, wie die weiteren Schritte ablaufen; es wird einfach die Bitfolge des zu transformierenden Index “umgedreht”! Diese Abbildung von {0, 1, 2, …, n−1 = 2^{p}−1} auf sich wird Bitinversion genannt. Sie stimmt natürlich mit ihrer eigenen Inversen überein (Idempotenz), und das ist ein ganz wesentlicher Vorteil gegenüber den einzelnen Sortierschritten. Denn letztere sind zwar auch über Permutationen von {0, 1, 2, …, n − 1 = 2^{p} − 1} zu schreiben, aber diese Permutationen sind nicht idempotent und erfordern zur Realisierung entweder Hilfsfelder oder komplizierte Prozeduren.

Die Bitinversion läßt sich nämlich wegen ihrer Idempotenz als eine reine Folge von einzelnen Indexvertauschungen schreiben, und die zugehörigen Umsortierungen der FFT-Ergebnisdaten lassen sich dann ohne ein größeres Hilfsarray durchführen. Man kann das Ganze so implementieren, daß man alle Indizes k von 0 bis n − 1 durchläuft und eine Vertauschung ausführt, wenn die Bitinverse von k größer ist als k. Die Abfrage vermeidet eine zweifache und damit insgesamt unwirksame Vertauschung. Die Bitinversion ist in der beschleunigten FFT–Routine, die auf dem Netz verfügbar ist, als separate Funktion (etwas ineffizient realisiert, aber hoffentlich verständlich) bereitgestellt.

Ein weiterer Beschleunigungstrick, der sich in vielen Programmen findet, setzt voraus, daß der Aufruf der Bibliotheksfunktionen sin und cos deutlich aufwendiger ist als 2 Additionen und 2 Multiplikationen. Denn dann lohnt es sich, die Potenzen ζ_{2m}^{j} aus (12.2.26) rekursiv aus ζ_{2m}^{1} über die Additionstheoreme von Sinus und Cosinus zu berechnen. Diese Rechnung ist nichts anderes als die Zerlegung von

ζ_{2m}^{j+1} = ζ_{2m}^{1}ζ_{2m}^{j}

in Real– und Imaginärteil. Auch das ist in dem beschleunigten FFT–Programm implementiert, obwohl der Beschleunigungseffekt diskutabel ist. Hat man genug Speicher, so sollte man die Potenzen der n–ten Einheitswurzeln vorher berechnen und aus einem array abrufen. Wegen der Kürzungsregel sind darin auch alle Einheitswurzeln zu Teilern von n enthalten.