Schnelle Fourier-Transformation

Die schnelle Fourier-Transformation (englisch fast Fourier transform, daher meist FFT abgekürzt) ist ein Algorithmus zur schnellen Berechnung der Werte einer diskreten Fourier-Transformation (DFT). Die Beschleunigung gegenüber der direkten Berechnung beruht auf der Vermeidung mehrfacher Berechnung sich gegenseitig aufhebender Terme. Der Algorithmus wird James W. Cooley und John W. Tukey zugeschrieben, die ihn 1965 veröffentlichten. Genaugenommen wurde eine Form des Algorithmus jedoch bereits 1805 von Carl Friedrich Gauss entworfen, der ihn zur Berechnung der Flugbahnen der Asteroiden Pallas und Juno verwendete. Darüber hinaus wurden eingeschränkte Formen des Algorithmus noch mehrfach vor Cooley und Tukey entwickelt, so z.B. von Good (1960). Nach Cooley und Tukey hat es darüber hinaus zahlreiche Verbesserungsvorschläge und Variationen gegeben, so etwa von Georg Bruun, C. M. Rader und Leo I. Bluestein.

Analog kennt man auch für die diskrete inverse Fourier-Transformation einen schnellen Algorithmus (inverse FFT - iFFT) Es kommen bei der iFFT idente Algorithmen mit anderen Koeffizienten in der Berechnung zur Anwendung.

Informelle Beschreibung des Algorithmus (Cooley und Tukey)

Der Algorithmus von Cooley und Tukey ist ein klassisches Teile-und-herrsche-Verfahren. Voraussetzung für seine Anwendung ist, dass die Anzahl der Stützstellen bzw. Abtastpunkte eine Zweierpotenz ist. Da die Anzahl solcher Punkte im Rahmen von Messverfahren jedoch im allgemeinen frei gewählt werden kann, handelt es sich dabei nicht um eine gravierende Einschränkung. Das Problem der Berechnung einer DFT der Größe n wird nun zunächst in zwei Berechnungen der DFT der Größe n/2 aufgeteilt - der Vektor der Messwerte wird in die Teilvektoren der Einträge zu geraden bzw. ungeraden Indizes aufgeteilt und von beiden die DFT bestimmt. Die beiden Teilergebnisse werden dann wieder zu einem Gesamtergebnis zusammengefügt. Zur Berechnung der beiden Hälften werden die Eigenschaften der Einheitswurzeln aus der Fouriermatrix ausgenutzt. Eine rekursive Anwendung dieser Grundidee ermöglicht schließlich eine Berechnung in O(n log n) Zeit.

Formelle Beschreibung des Algorithmus (Cooley und Tukey)

Zunächst sei in Erinnerung gerufen, dass die DFT durch folgende Formel definiert ist:

f_j =  \sum_{k=0}^{n-1} x_k \;e^{-\frac{2\pi i}{n} jk }

\qquad j = 0,\dots,n-1.

Setzen wir n' = n/2 und notieren die Einträge zu geraden Indizes als

x'₀ = x₀, x'₁ = x₂, ..., x'_n'-1 = x_2n-2

und bezeichnen deren Transformierte nach DFT der Größe n' mit

f'₀, ..., f'_n'-1;

die Einträge zu ungeraden Indizes notieren wir entsprechend als

x"₀ = x₁, x"₁ = x₃, ..., x"_n'-1 = x_2n-1

und bezeichnen deren Transformierte nach DFT der Größe n' mit

f"₀, ..., f"_n'-1.

Dann folgt:

\begin{matrix} f_j & = & \sum_{k=0}^{\frac{n}{2}-1} x_{2k} e^{-\frac{2\pi i}{n} j(2k)} + \sum_{k=0}^{\frac{n}{2}-1} x_{2k+1} e^{-\frac{2\pi i}{n} j(2k+1)} \\ \\ & = & \sum_{k=0}^{n'-1} x'_{k} e^{-\frac{2\pi i}{n'} jk} + e^{-\frac{2\pi i}{n}j} \sum_{k=0}^{n'-1} x''_k e^{-\frac{2\pi i}{n'} jk} \\ \\ & = & \left\{ \begin{matrix} f'_j + e^{-\frac{2\pi i}{n}j} f''_j & \mbox{falls } j

Mathematische Beschreibung (allgemeiner Fall)

In der Mathematik wird die schnelle diskrete Fouriertransformation in einem wesentlich allgemeineren Kontext behandelt:

Sei $R$ ein kommutativer unitärer Ring. In $R$ sei die Zahl $2=1+1$ eine Einheit (d.h. invertierbar); ferner sei $w\in R$ eine $2^n$ -te Einheitswurzel mit $w^{2^{n-1}}=-1$ .

Dann lässt sich im Modul $R^{2^n}$ zu $a\in R^{2^n}$ die diskrete Fouriertransformierte $\hat a$ mit

\hat a_k = \sum_{j=0}^{2^n-1} a_j\cdot w^{j\cdot k}\qquad (k=0,\ldots, 2^n-1)

wie folgt optimiert berechnen:

Zunächst stellen wir die Indizes $j$ und $k$ wie folgt dual dar:

k = \sum_{\nu=0}^{n-1} k_\nu 2^\nu,\quad j = \sum_{\nu=0}^{n-1} j_\nu 2^{n-1-\nu}

Damit haben wir folgende Rekursion:

A_0(j_0,\ldots,j_{n-1}) = a_j\qquad (j=0,\ldots,2^n-1)

A_{r+1}(k_0,\ldots, k_r, j_{r+1},\ldots,j_{n-1}) = \sum_{j_r=0}^1 A_r(k_0,\ldots,k_{r-1},j_r,\ldots,

j_{n-1})\cdot w^{j_r\cdot 2^{n-1-r}\cdot (k_r 2^r+\ldots+k_0 2^0)}.

Wegen

A_n(k_0,\ldots,k_{n-1}) = \hat a_k

erhalten wir hieraus die diskrete Fouriertransformierte

\hat a\in R^{2^n}

Komplexität

Die FFT ist im Gegensatz zur DFT nur möglich, wenn die Länge des Eingangsvektors einer Zweierpotenz entspricht (zumindest gilt dies für die sogenannte klassische Variante der FFT, wie sie oben definiert wurde). Die Anzahl der Abtastpunkte kann also beispielsweise 1, 2, 4, 8, 16, 32 usw. betragen. Man spricht hier von einer Radix-2-FFT. Wenn man die Längen auf Potenzen von vier beschränkt, kann man analog eine Radix-4-FFT einsetzen usw.. Der Rechenaufwand reduziert sich für eine Radix-2-FFT von

T \left( N \right)=O(N^2)

für die DFT auf

T \left( N \right)=O(N \cdot \log(N))

Bei einer Faltung zweier Signale reduziert sich der Rechenaufwand gar von

T \left( N \right)=O(N^4)

auf

T \left( N \right)=O(N^2 \cdot \log(N))

Die Struktur des Datenflusses kann durch einen Schmetterlingsgraphen beschrieben werden, der die Reihenfolge der Berechnung festlegt. Diese Struktur ist identisch zur Sortiermethode Quicksort.

Zur Abschätzung der Rechenleistung für die Berechnung eines Algorithmus siehe auch Komplexitätstheorie.

Alternative Formen der FFT

Neben den oben dargestellten FFT-Algorithmus von Cooley und Tukey, auch Radix-2-Algorithmus genannt, existieren noch eine Reihe weiterer Algorithmen zur schnellen Fourier-Transformation. Die Unterschiede liegen in der benötigten Anzahl der Rechenoperationen, so können die Anzahl der Multiplikationen in bestimmten Umfang durch schnellere Addition getauscht werden, und in den notwendigen Speicherplatz für das Bereithalten der zu transformierenden Daten.

Im folgenden sind übersichtsartig einige weitere Algorithmen dargestellt. Details und genaue mathematische Beschreibungen samt Herleitungen finden sich in der unten angegeben Literatur.

Radix-4-Algorithmus

Der Radix-4-Algorithmus ist, analog dazu der Radix-8-Algorithmus oder allgemein Radix-2^N-Algorithmus, eine Weiterentwicklung des obigen Radix-2-Algorithmus. Der Hauptunterschied besteht darin, dass die Anzahl der zu verarbeitenden Datenpunkte eine Potenz von 4 bzw. 2^N darstellen muss. Die Verarbeitungstruktur bleibt dabei gleich, nur dass in dem Schmetterlingsgraph pro Element statt 2 Datenpfade 4 Datenpfade miteinander verknüpft werden müssen. Der Vorteil besteht in einem weiter reduzierten Rechenaufwand und damit Geschwindigkeitsvorteil. So sind, verglichen mit dem obigen Algorithmus von Cooley und Tukey, bei dem Radix-4-Algorithmus ca. 25% weniger Multiplikationen notwendig. Bei dem Radix-8-Algorithmus reduziert sich die Anzahl der Multiplikationen um ca. 40%.

Nachteil dieser Verfahren ist die gröbere Struktur und ein aufwendiger Programmcode. So lassen sich mit Radix-4-Algorithmus nur Blöck der Längen 4, 16, 64, 256, 1024, 4096, ... verarbeiten. Bei dem Radix-8-Algorithmus ist dies analog zu sehen. Wäre in einer bestimmten Anwendung aber eine Blocklänge von beispielsweise 2048 Stützstellen ideal und damit auch Speicherplatz gering zu halten, können diese schnelleren Algorithmen daher nicht eingesetzt werden. Es müsste dann ein grösserer Speicher eingesetzt werden und der Rechenaufwand gesteigert werden.

Winograd-Algorithmus

Bei diesem Algorithmus ist nur eine bestimmte, endliche Anzahl von Stützstellen der Anzahl N möglich, nämlich:

N = \prod^{m}_{j=1}N_j \qquad N_j \in \{2,3,4,5,7,8,9,16\}

wobei alle Kombinationen von N_j zulässig sind, bei denen die verwendeten N_j relativ prim sind. Dadurch ist nur eine maximale Blöcklänge von 5040 möglich. Die möglichen Werte für N liegen aber in dem Bereich bis 5040 dichter auf der Zahlengeraden als die Zweierpotenzen. Es ist damit eine bessere Feinabstimmung der Blocklänge möglich. Aufgebaut wird der Algorithmus aus Basisblöcken der DFT deren Längen mit N_j korrespondieren. Bei diesem Verfahren wird zwar die Anzahl der Multiplikationen gegenüber dem Radix-2-Algorithmus reduziert, gleichzeitig steigt aber die Anzahl der notwendigen Additionen. Ausserdem ist am Eingang und Ausgang jeder DFT eine aufwendige Permutation der Daten notwendig, die nach den Regeln des chinesischen Restsatzes gebildet wird.

Diese Art der schnellen Fourier-Transformation besitzt in praktischen Implementierungen dann Vorteile gegenüber der Radix-2 Methode, wenn der für die FFT verwendete Mikrocontroller, wie beispielsweise ein Intel 8051 ohne eigene Multipliziereinheit, für die notwendigen Multiplikationen sehr viel Rechenzeit in Form von Unterprogrammen benötigen würde. In heutigen Signalprozessoren mit eigenen Multipliziereinheiten hat dieser Algorithmus keine wesentliche Bedeutung.

Primfaktor-Algorithmus

Dieser FFT-Algorithmus basiert auf ähnlichen Ideen wie der Winograd-Algorithmus, allerdings ist die Struktur einfacher und damit der Aufwand an Multiplikationen höher als wie beim Winograd-Algorithmus. Der wesentliche Vorteil bei der Implementierung liegt in der effizienten Ausnützung des zur Verfügung stehenden Speichers durch optimale Anpassung der Blocklänge. Wenn in einer bestimmten Anwendung zwar eine schnelle Multipliziereinheit verfügbar ist und gleichzeitig der verfügbare Speicherplatz sehr beschränkt ist, kann dieser Algorithmus das Optimum bieten. Die Ausführungszeit ist bei in etwa gleichen Blöcklänge in etwa mit FFT-Algorithmus von Cooley und Tukey vergleichbar.

Anwendung

Die Kombination aus FFT und iFFT kann zur Kodierung und Dekodierung von Signalen auf Frequenzebene eingesetzt werden. Kompressionsalgorithmen wie der des MP3-Formats basieren hierauf. Ein weiteres wichtiges Anwendungsbeispiel ist die Breitbanddatenübertragung per OFDM, die Grundlage für ADSL und WLAN (Internet), DVB-T (Fernsehen), DRM und DAB (Radio) ist.

Basis für die Signalanalyse
Schwingungsmesstechnik - Umrechnung Zeit in Frequenzendarstellung
Akustik (Audiomessungen) - Berechnung von Spektrogrammen (Diagramme mit der Darstellung der Amplituden von den jeweiligen Frequenzanteilen)
Längstwellenempfang mit dem PC
Teil schneller, Polynome verarbeitender Algorithmen (z.B. Polynommultiplikation in $O(n \log n)$ ) in der Computeralgebra.
Zur Reduzierung des Berechnungsaufwandes bei der zirkularen Faltung im Zeitbereich von FIR-Filtern und Ersatz durch die schnelle Fouriertransformation mit einfachen Multiplikationen im Frequenzbereich. (siehe auch Schnelle Faltung)

FFT-Software zur Signalanalyse mit dem PC mit Hilfe der PC-Soundkarte

prettyFastFFT (OpenGL, Freeware)
Sonogram Visible Sound (OpenSource, Java-basiert)
SpecLab (Freeware)
Spectran (Freeware)
Spectrogram (Shareware)
Audiotester (Shareware)
SIGVIEW (Shareware)

Literatur

James W. Cooley und John W. Tukey: An algorithm for the machine calculation of complex Fourier series, Math. Comput. 19, 297–301 (1965).
Georg Bruun: z-Transform DFT filters and FFTs, IEEE Trans. on Acoustics, Speech and Signal Processing (ASSP) 26 (1), 56-63 (1978).
C. M. Rader: Discrete Fourier transforms when the number of data samples is prime, Proc. IEEE 56, 1107–1108 (1968).
Leo I. Bluestein: A linear filtering approach to the computation of the discrete Fourier transform, Northeast Electronics Research and Engineering Meeting Record 10, 218-219 (1968).
Heideman, Johnson, Burrus: Gauss and the History of the Fast Fourier Transform, Arch. Hist. Sc. 34(3) (1985).
Oppenheim: Zeitdiskrete Signalverarbeitung, Oldenbourg Verlag 1999, ISBN 3-486-24145-1

Weblinks

http://www.fftw.org
http://www.nr.com/ (Buch Numerical Recipes, engl., auch online verfügbar)
http://www.unet.univie.ac.at/~a0503736/php/floxwiki/pmwiki.php/Main/ProjekteAmpUni Bakkalaureatsarbeit (pdf)

Numerische Mathematik | Digitale Signalverarbeitung | Algorithmus

Gourt Home::Gourt :: Schnelle Fourier-Transformation