TransformationenSchnelle Fouriertransformation (FFT) | |
Die Fouriertransformation ist ein fundamentales Verfahren in der Signalverarbeitung. Durch die Fouriertransformation lassen sich Signale von der Darstellung {(Zeitpunkt, Abtastwert)} in die Darstellung {(Frequenzanteil, Amplitude, Phase)} überführen. Viele Operationen, z.B. Filter, lassen sich im Frequenzraum leichter durchführen. Anschließend wird das Signal mit der inversen Fouriertransformation wieder zurück transformiert. [Beispiel]
Außerdem hat die Fouriertransformation Anwendungen in der Numerik; z.B. lässt sich eine Polynommultiplikation schneller durchführen, wenn die Polynome in Stützstellendarstellung statt in Koeffizientendarstellung vorliegen.
Das Verfahren der schnellen Fouriertransformation (engl.: Fast Fourier Transform – FFT) hat eine Zeitkomplexität von O(n log(n)). Dadurch ist die Polynommultiplikation sogar einschließlich Transformation in Stützstellendarstellung und Rücktransformation noch schneller als die direkte Multiplikation in Koeffizientendarstellung.
Definition: Sei
der Körper der komplexen Zahlen. Ein Element w
heißt n-te Einheitswurzel, wenn wn = 1 ist, und w heißt primitive n-te Einheitswurzel, wenn wn = 1 ist, aber wk
1 für alle k
{1, ..., n-1}.
Beispiel: Sei n = 4. Dann ist i primitive 4-te Einheitswurzel1). Alle 4-ten Einheitswurzeln sind:
i0 = 1, i1 = i, i2 = -1, i3 = -i.
Für n = 6 ist cos(2π/6) + i sin(2π/6) primitive 6-te Einheitswurzel (Bild 1).
Allgemein gilt in
:
w = cos(k·2π/n) + i sin(k·2π/n) mit k
{0, ..., n-1}
ist n-te Einheitswurzel; ist k teilerfremd zu n, so ist w primitiv.
In der eulerschen Schreibweise lässt sich eine n-te Einheitswurzeln w auch ausdrücken als
w = eik2π/n mit k
{0, ..., n-1}.
| |
| Bild 1: Einheitskreis in der gaußschen Zahlenebene mit 4-ten und 6-ten Einheitswurzeln | |
genau n verschiedene n-te Einheitswurzeln, diese sind darstellbar als die Potenzen einer primitiven n-ten Einheitswurzel w: w0, w1, w2, ..., wn-1.
(wk)n = wk·n = (wn)k = 1k = 1.
Dies gilt auch für negative k.
wn/2 = -1,
denn (wn/2)2 = wn = 1, d.h. wn/2 ist 2-te Einheitswurzel, also 1 oder -1. Da aber wn/2
1 ist, da w primitiv ist, gilt wn/2 = -1.
{1, ..., n/2-1} mit (w2)k = 1. Dann ist aber w2k = 1 mit 2k < n, im Widerspruch dazu, dass w primitiv ist.
{1, ..., n-1} mit (w-1)k = w-k = 1. Dann ist aber wk = 1/w-k = 1/1 = 1, im Widerspruch dazu, dass w primitiv ist.
gilt für die zu einer n-ten Einheitswurzel w konjugierte Zahl w w = w-1
denn es ist
| w · w | = | (cos(k·2π/n) + i sin(k·2π/n)) · (cos(k·2π/n) – i sin(k·2π/n)) | ||
| = | cos(k·2π/n)2 + sin(k·2π/n)2 | |||
| = | 1. |
Definition: Sei n
und w primitive n-te Einheitswurzel in
. Eine n
n-Matrix F mit
Fi,j = wi·j
für alle i, j
{0, ..., n-1} heißt Fouriermatrix.2)
Die lineare Abbildung f :
n
n mit
f(a) = a·F
für alle (Zeilen-)vektoren a
n heißt diskrete Fouriertransformation (DFT).3)
Beispiel: Sei n = 4. Dann ist i primitive n-te Einheitswurzel. Die zugehörige Fouriermatrix ist
| F = | ![]() |
| ![]() |
So ist etwa die -1 in der letzten Zeile der Matrix das Element
F3,2 = w3·2 = w6 = (-1)3 = -1
(die Elemente der Matrix werden von 0 bis n-1 indiziert).
Die Fouriertransformation des Vektors a = [1 1 1 0] ergibt
y = a·F = [3 i 1 -i]
Offenbar ist die Fouriermatrix F die Vandermonde-Matrix des Vektors w0, ..., wn-1. Die Matrix-Vektor-Multiplikation y = a·F lässt sich somit als Polynomauswertung an den Stellen w0, ..., wn-1 auffassen, wobei der Vektor a die Koeffizienten des Polynoms enthält. Das Ergebnis y0, ..., yn-1 ist die Stützstellendarstellung des Polynoms.
Satz: Die inverse Fouriermatrix F -1 existiert und ist gleich
F -1i,j = 1/n · w-i·j
für alle i, j
{0, ..., n-1}. Die inverse Fouriermatrix enthält also die zu den Elementen der Fouriermatrix inversen Elemente, dividiert durch n.
Definition: Die lineare Abbildung f -1 :
n
n mit
f -1(a) = a·F -1
für alle a
n heißt inverse Fouriertransformation.
Beispiel: Sei n = 4. Die inverse Fouriermatrix ist
| F -1 = 1/4 · | ![]() |
| ![]() |
Die inverse Fouriertransformation des Vektors y = [3 i 1 -i] ergibt
a = y·F -1 = [1 1 1 0].
In dieser Form ist die Fouriertransformation eine Matrix-Vektor-Multiplikation mit der Komplexität O(n2). Durch Ausnutzung der Symmetrie der n-ten Einheitswurzeln lässt sich die Berechnung auf O(n log(n)) beschleunigen. Dieses Verfahren heißt schnelle Fouriertransformation (Fast Fourier Transform – FFT) [CT 65].
Die Idee des Verfahrens der schnellen Fouriertransformation ist, die einzelnen Berechnungen der Matrix-Vektor-Multiplikation y = a·F in einer speziellen Reihenfolge auszuführen, so dass jeweils auf schon berechnete Zwischenergebnisse zurückgegriffen werden kann. Dabei werden die o.a. Eigenschaften (3) und (4) ausgenutzt, nämlich dass wn/2 = -1 ist und dass das Quadrat w2 der primitiven n-ten Einheitswurzel w primitive n/2-te Einheitswurzel ist. Das Verfahren setzt voraus, dass n eine Zweierpotenz ist.
Zunächst werden die Komponenten von y mit geradem Index berechnet, indem der Vektor a mit den entsprechenden Spalten der Fouriermatrix multipliziert wird. Es gilt für alle k
{0, ..., n/2-1}:
yk' = y2k =
i = 0, ..., n-1 ai wi·2k.
Die Summe wird in zwei Hälften aufgespalten:
yk' =
i = 0, ..., n/2-1 ai wi·2k +
i = 0, ..., n/2-1 ai+n/2 w(i+n/2)·2k.
Nun ist aber
w(i+n/2)·2k = wi·2k + nk = wi·2k·wnk = wi·2k,
da wnk = 1 ist, und damit
yk' =
i = 0, ..., n/2-1 (ai + ai+n/2) wi·2k.
Mit m = n/2 sowie v = w2, v primitive m-te Einheitswurzel, gilt:
yk' =
i = 0, ..., m-1 (ai + ai+m) vi·k,
d.h. yk' ist nichts anderes als die k-te Komponente der Fouriertransformation des Vektors
(ai + ai+m) i = 0, ..., m-1
der Länge m.
Ähnlich werden die Komponenten von y mit ungeradem Index berechnet. Es gilt für alle k
{0, ..., n/2-1} :
yk'' = y2k+1 =
i = 0, ..., n-1 ai wi·(2k+1).
Wiederum wird die Summe in zwei Hälften aufgespalten:
yk'' =
i = 0, ..., n/2-1 ai wi·(2k+1) +
i = 0, ..., n/2-1 ai+n/2 w(i+n/2)·(2k+1).
Nun ist aber
w(i+n/2)·(2k+1) = wi·2k + nk + i + n/2 = - wi·wi·2k,
da wnk = 1 ist und wn/2 = -1 ist. Somit gilt:
yk'' =
i = 0, ..., n/2-1 ai wi·wi·2k +
i = 0, ..., n/2-1 – ai+n/2 wi·wi·2k
=
i = 0, ..., n/2-1 wi·(ai – ai+n/2) wi·2k.
Mit m = n/2 sowie v = w2 gilt wiederum:
yk'' =
i = 0, ..., m-1 wi·(ai – ai+m) vi·k,
d.h. yk'' ist nichts anderes als die k-te Komponente der Fouriertransformation des Vektors
wi·(ai – ai+m) i = 0, ..., m-1.
Durch rekursive Anwendung dieses Verfahrens auf Vektoren jeweils halber Länge wird im Ergebnis die Fouriertransformation berechnet. Bild 2 zeigt schematisch den Ablauf der Berechnung für n = 8.
| |
| Bild 2: Datenfluss der schnellen Fouriertransformation | |
Um die Fouriertransformation eines Vektors a zu berechnen, ist zunächst ein Vektor a' zu berechnen, dessen Komponenten ai' für i = 0, ..., m-1 wie oben gesehen folgende Werte haben:
ai' = ai + ai+m sowie
ai+m' = wi·(ai – ai+m).
Hierfür sind m Additionen, m Subtraktionen und m Multiplikationen erforderlich sowie nochmals m Multiplikationen, um jeweils wi aus wi-1 zu berechnen. Insgesamt ergibt dies also 2m = n Additionen und 2m = n Multiplikationen.
Anschließend ist rekursiv auf die beiden Hälften von a' die Fouriertransformation anzuwenden.
Die Ergebnisvektoren y' und y'' stellen die geraden und die ungeraden Komponenten von y dar; sie sind noch ineinander zu verschränken (perfect shuffle), um den gewünschten Vektor y zu erhalten. Die Permutation perfect shuffle lässt sich, etwa mit der unten angegebenen Prozedur, in 3/2 n Schritten durchführen.
Die Zeitkomplexität T(n) der FFT ergibt sich somit als
T(n) = 3.5n + 2·T(n/2) sowie
T(1) = 0.
Die Auflösung dieser Rekursion ergibt
T(n) = 3.5n·log(n), d.h.
T(n)
O(n log(n)).
Die folgende Prozedur berechnet die Fouriertransformation eines komplexen Vektors a, beginnend beim Index lo und der Länge n. Der Parameter w steht für die primitive n-te Einheitswurzel. Rechenoperationen mit komplexen Zahlen sind der Übersichtlichkeit halber mit normalen Rechenzeichen (+, -, *) dargestellt, obwohl diese in Java eigentlich nicht zur Verfügung stehen.
void fft(Complex[] a, int n, int lo, Complex w) { int i, m; Complex z, v, h; if (n>1) { m=n/2; z=1; for (i=lo; i<lo+m; i++) { h=a[i]-a[i+m]; a[i]=a[i]+a[i+m]; a[i+m]=h*z; z=z*w; } v=w*w; fft(a, m, lo, v); fft(a, m, lo+m, v); shuffle (a, n, lo); } } |
Die Prozedur shuffle verschränkt die beiden, von den rekursiven Aufrufen von fft erzeugten Hälften reißverschlussartig ineinander. Die entsprechende Permutation für n = 8 lautet
| 0 1 2 3 4 5 6 7 |
| 0 4 1 5 2 6 3 7 |
Zur Ausführung der Permutation wird ein Hilfsarray b verwendet, in das zunächst die eine Hälfte der Folge ausgelagert wird.
void shuffle(Complex[] a, int n, int lo) { int i, m=n/2; Complex[] b=new Complex[m]; for (i=0; i<m; i++) b[i]=a[lo+i]; for (i=0; i<m; i++) a[lo+i+i+1]=a[lo+i+m]; for (i=0; i<m; i++) a[lo+i+i]=b[i]; } |
Die inverse Fouriertransformation lässt sich mit demselben Verfahren durchführen. Aufgrund der Definition der inversen Fouriermatrix F -1 wird jedoch statt mit der primitiven n-ten Einheitswurzel w mit der inversen n-ten Einheitswurzel w-1 gearbeitet. In
ist dies die konjugiert komplexe n-te Einheitswurzel w (vgl. o.a. Eigenschaften (5) und (6)). Ferner werden die Elemente der invers zu transformierenden Folge zunächst durch n geteilt.
Die diskrete Fouriertransformation lässt sich interpretieren als Transformation eines Polynoms vom Grad n-1 von der Koeffizientendarstellung in die Stützstellendarstellung. Als Stützstellen werden die n-ten Einheitswurzeln w0, w1, ..., wn-1 verwendet. Im Körper
der komplexen Zahlen sind die Werte
wk = cos(k·2π/n) + i sin(k·2π/n) = eik2π/n mit k
{0, ..., n-1}
die n-ten Einheitswurzeln.
Durch Ausnutzung der Symmetrien der n-ten Einheitswurzeln lässt sich die Fouriertransformation beschleunigen; das Verfahren ist die schnelle Fouriertransformation FFT.
Mit Hilfe der FFT lässt sich die Zeitkomplexität der Polynommultiplikation von Θ(n2) auf Θ(n log(n)) verbessern.
| [AHU 74] | A.V. Aho, J.E. Hopcroft, J.D. Ullman: The Design and Analysis of Computer Algorithms. Addison-Wesley (1974) | |
| [CT 65] | J.M. Cooley, J.W. Tukey: An Algorithm for the Machine Calculation of Complex Fourier Series. Math. Comp. 19, 297-301 (1965) | |
| [Sed 88] | R. Sedgewick: Algorithms. 2. Auflage, Addison-Wesley (1988) | |
| [Lan 06] | H.W. Lang: Algorithmen in Java. 2. Auflage, Oldenbourg (2006)
Den FFT-Algorithmus finden Sie auch in meinem Buch über Algorithmen. |
|
1) Die Zahl i ist die imaginäre Einheit i =
-1
2) Da es für n>2 stets mehrere primitive n-te Einheitswurzeln gibt, ist die Fouriermatrix insofern nicht eindeutig festgelegt.
3) Da die Fouriermatrix symmetrisch ist, lässt sich die Fouriertransformation auch als f(a) = F·a für Spaltenvektoren a definieren.
Weiter mit: [Diskrete Kosinus-Transformation] oder ![]() |
|
|