Schnelle Fourier-Transformation

Zeit-basierte Darstellung (oben) und Frequenz-basierte Darstellung (unten) desselben Signals, wobei die untere Darstellung aus der oberen durch Fouriertransformation gewonnen werden kann.

Die schnelle Fourier-Transformation (englisch fast Fourier transform, daher meist FFT abgekürzt) ist ein Algorithmus zur effizienten Berechnung der diskreten Fourier-Transformation (DFT). Mit ihr kann ein digitales Signal in seine Frequenzanteile zerlegt und diese dann analysiert werden.

Analog gibt es für die diskrete inverse Fourier-Transformation die inverse schnelle Fourier-Transformation (IFFT). Es kommen bei der IFFT die gleichen Algorithmen, aber mit konjugierten Koeffizienten zur Anwendung.

Die FFT hat zahlreiche Anwendungen (siehe auch unten) im Bereich der Ingenieurwissenschaften, der Naturwissenschaften und der angewandten Mathematik. Außerdem kommt sie in Mobilfunktechnologien wie UMTS und LTE und bei der drahtlosen Datenübertragung zum Einsatz, etwa in der WLAN-Funknetztechnik.

Die FFT gehört zu den Teile-und-herrsche-Verfahren, sodass – im Gegensatz zur direkten Berechnung – zuvor berechnete Zwischenergebnisse wiederverwendet und dadurch arithmetische Rechenoperationen eingespart werden können. Das bekannteste Verfahren wird James Cooley und John W. Tukey zugeschrieben, die es 1965 veröffentlichten. Genau genommen wurde eine Form des Algorithmus bereits 1805 von Carl Friedrich Gauß entworfen, der ihn zur Berechnung der Flugbahnen der Asteroiden (2) Pallas und (3) Juno verwendete. Zum ersten Mal publiziert wurde eine Variante des Algorithmus von Carl Runge im Jahre 1903 und 1905. Darüber hinaus wurden eingeschränkte Formen des Algorithmus mehrfach vor Cooley und Tukey entwickelt, so z.B. von Irving John 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.

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.

Der Algorithmus basiert auf der Beobachtung, dass die Berechnung einer DFT der Größe 2n in zwei Berechnungen einer DFT der Größe n zerlegbar ist (über den Vektor mit den Einträgen der geraden bzw. der ungeraden Indizes), wobei die beiden Teilergebnisse nach der Transformation wieder zu einer Fouriertransformation der Größe 2n zusammenzufassen sind.

Da die Berechnung einer DFT der halben Länge nur ein Viertel der komplexen Multiplikationen und Additionen der originalen DFT benötigt, und je nach Länge des Ausgangsvektors diese Vorschrift mehrfach hintereinander anwendbar ist, erlaubt die rekursive Anwendung dieser Grundidee schließlich eine Berechnung in {\mathcal  {O}}(n\log n) Zeit; zur Einsparung von trigonometrischen Rechenoperationen können bei der FFT zusätzlich die Eigenschaften der Einheitswurzeln aus der Fouriermatrix ausgenutzt werden.

Algorithmus von Cooley und Tukey

Die diskrete Fouriertransformation (DFT) eines Vektors (x_{0},\dotsc ,x_{{2n-1}}) der Dimension 2n lautet:

f_{m}=\sum _{{k=0}}^{{2n-1}}x_{k}\,e^{{-{\frac  {2\pi i}{2n}}mk}}\qquad m=0,\dotsc ,2n-1.

Die Einträge mit geraden Indizes werden notiert als

x'_{k}=x_{{2k}},\quad k=0,\dotsc ,n-1

und deren DFT der Dimension n als (f'_{k}).

Entsprechend seien die Einträge mit ungeraden Indizes notiert als

x''_{k}=x_{{2k+1}},\quad k=0,\dotsc ,n-1

mit DFT (f''_{k}).

Dann folgt:

{\begin{aligned}f_{m}&=\sum _{{k=0}}^{{n-1}}x_{{2k}}e^{{-{\frac  {2\pi i}{2n}}m(2k)}}+\sum _{{k=0}}^{{n-1}}x_{{2k+1}}e^{{-{\frac  {2\pi i}{2n}}m(2k+1)}}\\[0.5em]&=\sum _{{k=0}}^{{n-1}}x'_{k}e^{{-{\frac  {2\pi i}{n}}mk}}+e^{{-{\frac  {\pi i}{n}}m}}\sum _{{k=0}}^{{n-1}}x''_{k}e^{{-{\frac  {2\pi i}{n}}mk}}\\[0.5em]&={\begin{cases}f'_{m}+e^{{-{\frac  {\pi i}{n}}m}}f''_{m}&{\text{ falls }}m<n\\[0.5em]f'_{{m-n}}-e^{{-{\frac  {\pi i}{n}}(m-n)}}f''_{{m-n}}&{\text{ falls }}m\geq n\end{cases}}\end{aligned}}


Mit der Berechnung von {\displaystyle f'_{m}} und {\displaystyle f''_{m}} ({\displaystyle m=0,\dotsc ,n-1}) ist sowohl f_{m}, als auch f_{{m+n}} bestimmt. Der Rechenaufwand hat sich durch diese Zerlegung also praktisch halbiert.

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. Zum Beispiel ist im Restklassenring

R=\mathbb{Z } /N\mathbb{Z } mit N=2^{{d2^{n}}}+1, d,n\in \mathbb{N} , d ungerade (das ist gleichbedeutend mit der Forderung „teilerfremd zu 2^{n}“),

das Element w=2^{{2d}} eine solche Einheitswurzel, die entsprechende FFT wird im Schönhage-Strassen-Algorithmus verwendet.

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

Diese klassische Variante der FFT nach Cooley und Tukey ist im Gegensatz zur DFT nur durchführbar, wenn die Länge des Eingangsvektors einer Zweierpotenz entspricht. Die Anzahl der Abtastpunkte kann also beispielsweise 1, 2, 4, 8, 16, 32 usw. betragen. Man spricht hier von einer Radix-2-FFT. Andere Längen sind mit den unten angeführten alternativen Algorithmen möglich.

Aus obiger Rekursion ergibt sich folgende Rekursionsgleichung für die Laufzeit der FFT:

T\left(N\right)=2T\left(N/2\right)+f(N)

Hierbei beschreibt der Term f(N) den Aufwand, um die Ergebnisse mit einer Potenz der Einheitswurzel zu multiplizieren und die Ergebnisse zu addieren. Es werden N Paare von Zahlen addiert und N/2 Zahlen mit Einheitswurzeln multipliziert. Insgesamt ist f(N) also linear beschränkt:

T\left(N\right)=2T\left(N/2\right)+{\mathcal  {O}}\left(N\right)

Mit dem Master-Theorem ergibt sich eine Laufzeit von:

T\left(N\right)={\mathcal  {O}}(N\cdot \log(N))

Die Struktur des Datenflusses kann durch einen Schmetterlingsgraphen beschrieben werden, der die Reihenfolge der Berechnung festlegt.

Implementierung als rekursiver Algorithmus

Die direkte Implementierung der FFT in Pseudocode nach obiger Vorschrift besitzt die Form eines rekursiven Algorithmus:

Dies wird nun fortgeführt, bis das Argument eines Aufrufs der Funktion nur noch aus einem einzigen Element besteht (Rekursionsabbruch): Die FFT eines einzelnen Wertes ist (er besitzt sich selbst als Gleichanteil, und keine weiteren Frequenzen) er selbst. Die Funktion, die nur noch einen einzigen Wert als Parameter erhält, kann also ganz ohne Rechnung die FFT dieses Wertes zurückliefern – die Funktion, die sie aufgerufen hat, kombiniert die beiden jeweils 1 Punkt langen FFTs, die sie zurückerhält, die Funktion, die diese wiederum aufgerufen hat, die beiden 2-Punkte-FFTs, und so weiter.

{\mathrm  {function}}\ \operatorname {fft}(n,{\vec  f}):

{\mathrm  {if}}\ (n=1)
{\mathrm  {return}}\ {\vec  f}
{\mathrm  {else}}\
{\vec  g}=\operatorname {fft}\left({\tfrac  {n}{2}},(f_{0},f_{2},\ldots ,f_{{n-2}})\right)
{\vec  u}=\operatorname {fft}\left({\tfrac  {n}{2}},(f_{1},f_{3},\ldots ,f_{{n-1}})\right)
{\mathrm  {for}}\ k=0\ {\mathrm  {to}}\ {\tfrac  {n}{2}}-1
{\begin{aligned}c_{k}=g_{k}+u_{k}\cdot e^{{-2\pi ik/n}}\\c_{{k+n/2}}=g_{k}-u_{k}\cdot e^{{-2\pi ik/n}}\end{aligned}}
{\mathrm  {return}}\ {\vec  c}

Der Geschwindigkeitsvorteil der FFT gegenüber der DFT kann anhand dieses Algorithmus gut abgeschätzt werden:

Implementierung als nichtrekursiver Algorithmus

Die Implementierung eines rekursiven Algorithmus ist im Regelfall vom Ressourcenverbrauch her nicht ideal, da die vielen dabei notwendigen Funktionsaufrufe Rechenzeit und Speicher für das Merken der Rücksprungadressen benötigen. In der Praxis wird daher meist ein nichtrekursiver Algorithmus verwendet, der gegenüber der hier abgebildeten, auf einfaches Verständnis optimierten Form je nach Anwendung noch optimiert werden kann:

{\text{links}}=2^{{{\text{Ebene}}+1}}\cdot {\text{Abschnitt}}+{\text{Element des Abschnitts}}
{\text{rechts}}={\text{links}}+2^{{{\text{Ebene}}}}
über einen Schmetterlingsgraph kombiniert:
{\begin{array}{rcl}x_{{{\text{links,neu}}}}&=&x_{{{\text{links}}}}+e^{{-i\pi {\frac  {{\text{Element des Abschnitts}}}{2^{{{\text{Ebene}}}}}}}}\cdot x_{{{\text{rechts}}}}\\x_{{{\text{rechts,neu}}}}&=&x_{{{\text{links}}}}-e^{{-i\pi {\frac  {{\text{Element des Abschnitts}}}{2^{{{\mathrm  {Ebene}}}}}}}}\cdot x_{{{\text{rechts}}}}\end{array}}

Alternative Formen der FFT

Neben dem oben dargestellten FFT-Algorithmus von Cooley und Tukey, auch Radix-2-Algorithmus genannt, existieren noch eine Reihe weiterer Algorithmen zur schnellen Fourier-Transformation. Die Varianten unterscheiden sich darin, wie bestimmte Teile des „naiven“ Algorithmus so umgeformt werden, dass weniger (Hochpräzisions-)Multiplikationen notwendig sind. Dabei gilt meist, dass die Reduktion in der Anzahl der Multiplikationen eine erhöhte Anzahl von Additionen sowie von gleichzeitig im Speicher zu haltenden Zwischenergebnissen hervorruft.

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

Radix-4-Algorithmus

Der Radix-4-Algorithmus ist, analog dazu der Radix-8-Algorithmus oder allgemein Radix-2N-Algorithmus, eine Weiterentwicklung des obigen Radix-2-Algorithmus. Der Hauptunterschied besteht darin, dass die Anzahl der zu verarbeitenden Datenpunkte eine Potenz von 4 bzw. 2N darstellen muss. Die Verarbeitungstruktur bleibt dabei gleich, nur dass in dem Schmetterlingsgraph pro Element statt zwei Datenpfade vier bzw. acht und allgemein 2N 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öcke der Längen 4, 16, 64, 256, 1024, 4096, … verarbeiten. Bei dem Radix-8-Algorithmus sind die Einschränkungen analog zu sehen.

Winograd-Algorithmus

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

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

wobei alle Kombinationen von Nj zulässig sind, bei denen die verwendeten Nj teilerfremd sind. Dadurch ist nur eine maximale Blocklä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 Nj 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. Außerdem 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 keine eigene Multipliziereinheit besitzt und für die Multiplikationen sehr viel Rechenzeit aufgewendet werden muss. In heutigen Signalprozessoren mit eigenen Multipliziereinheiten hat dieser Algorithmus keine wesentliche Bedeutung mehr.

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 beim Winograd-Algorithmus. Der wesentliche Vorteil bei der Implementierung liegt in der effizienten Ausnutzung 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 Speicher knapp, kann dieser Algorithmus optimal sein. Die Ausführungszeit ist bei ähnlicher Blocklänge mit der des Algorithmus von Cooley und Tukey vergleichbar.

Goertzel-Algorithmus

Der Goertzel-Algorithmus stellt eine besondere Form zur effizienten Berechnung einzelner Spektralkomponenten dar und ist bei der Berechnung von nur einigen wenigen Spektralanteilen (engl. Bins) effizienter als alle blockbasierenden FFT-Algorithmen, welche immer das komplette diskrete Spektrum berechnen.

Chirp-z-Transformation

Bluestein-FFT-Algorithmus für Datenmengen beliebiger Größe (einschließlich Primzahlen).

Die inverse FFT

Die Inverse der diskreten Fourier-Transformation (DFT) stimmt bis auf den Normierungsfaktor und ein Vorzeichen mit der DFT überein. Da die schnelle Fourier-Transformation ein Algorithmus zur Berechnung der DFT ist, gilt dies dann natürlich auch für die FFT.

Anwendung

Die Anwendungsgebiete der FFT sind so vielseitig, dass hier nur eine Auswahl wiedergegeben werden kann:

Trenner
Basierend auf einem Artikel in: Wikipedia.de
Seitenende
Seite zurück
©  biancahoegel.de
Datum der letzten Änderung:  Jena, den: 28.05. 2023