Fast Fourier Transform (FFT)
Contents
Rechenaufwand von DFT bzw. IDFT
Ein Nachteil der direkten Berechnung der (im Allgemeinen komplexen) DFT–Zahlenfolgen
- ⟨D(μ)⟩∙−−(N)−−∘⟨d(ν)⟩
gemäß den in Kapitel Diskrete Fouriertransformation (DFT) angegebenen Gleichungen ist der große Rechenaufwand. Wir betrachten als Beispiel die DFT, also die Berechnung der D(μ) aus den d(ν):
- N⋅D(μ)=N−1∑ν=0d(ν)⋅wν⋅μ=d(0)⋅w0+d(1)⋅wμ+d(2)⋅w2μ+ ...+d(N−1)⋅w(N−1)⋅μ
Der hierfür erforderliche Rechenaufwand soll abgeschätzt werden, wobei wir davon ausgehen, dass die Potenzen des komplexen Drehfaktors w=e−j⋅2π/N bereits in Real– und Imaginärteilform in einer Lookup–Tabelle vorliegen. Zur Berechnung eines einzelnen Koeffizienten benötigt man dann N−1 komplexe Multiplikationen und ebenso viele komplexe Additionen, wobei zu beachten ist:
- Jede komplexe Addition erfordert zwei reelle Additionen:
- (R1+j⋅I1)+(R2+j⋅I2)=(R1+R2)+j⋅(I1+I2).
- Jede komplexe Multiplikation erfordert vier reelle Multiplikationen und zwei reelle Additionen (eine Subtraktion wird wie eine Addition behandelt):
- (R1+j⋅I1)(R2+j⋅I2)=(R1⋅R2−I1⋅I2)+j⋅(R1⋅I2+R2⋅I1).
- Somit sind zur Berechnung aller N Koeffizienten insgesamt die folgende Anzahl M reeller Multiplikationen und die Anzahl A reeller Additionen erforderlich:
- M=4⋅N⋅(N−1),
- A=2⋅N⋅(N−1)+2⋅N⋅(N−1)=M.
- In heutigen Rechnern benötigen Multiplikationen und Additionen/Subtraktionen etwa die gleiche Rechenzeit. Es genügt, die Gesamtzahl O=M+A aller Operationen zu betrachten:
- O=8⋅N⋅(N−1)≈8⋅N2.
Fazit:
- Man benötigt bereits für eine Diskrete Fouriertransformation (DFT) mit N=1000 knapp acht Millionen Rechenoperationen. Gleiches gilt für eine IDFT.
- Mit N=16 sind immerhin noch 1920 Rechenoperationen erforderlich.
Ist der Parameter N eine Potenz zur Basis 2, so können rechenzeitgünstigere Algorithmen angewendet werden. Die Vielzahl solcher aus der Literatur bekannten Verfahren werden unter dem Sammelbegriff Fast–Fouriertransformation – abgekürzt FFT – zusammengefasst. Alle diese Methoden basieren auf dem Überlagerungssatz der DFT.
Überlagerungssatz der DFT
Die Grafik verdeutlicht den so genannten Überlagerungssatz der DFT am Beispiel N=16. Dargestellt ist hier der Übergang vom Zeit– in den Spektralbereich, also die Berechnung der Spektralbereichskoeffizienten aus den Zeitbereichskoeffizienten: ⟨D(μ)⟩∙−−(N)−−∘⟨d(ν)⟩.
Der dadurch beschriebene Algorithmus ist durch folgende Schritte gekennzeichnet:
- Die Folge ⟨d(ν)⟩ der Länge N wird in zwei Teilfolgen ⟨d1(ν)⟩ und ⟨d2(ν)⟩ jeweils halber Länge separiert (in der Garafik gelb bzw. grün hinterlegt). Mit 0≤ν<N/2 erhält man so die Folgenelemente
- d1(ν)=d(2ν),
- d2(ν)=d(2ν+1).
- Die Ausgangsfolgen ⟨D1(μ)⟩ und ⟨D2(μ)⟩ der beiden Teilblöcke ergeben sich daraus jeweils durch eine eigene DFT, aber nun nur noch mit halber Länge N/2=8:
- ⟨D1(μ)⟩∙−−(N/2)−−∘⟨d1(ν)⟩,
- ⟨D2(μ)⟩∙−−(N/2)−−∘⟨d2(ν)⟩.
- Die Ausgangswerte ⟨D2(μ)⟩ der unteren (grünen) DFT (mit 0≤μ<N/2) werden danach im rot umrandeten Block durch komplexe Drehfaktoren hinsichtlich der Phasenlage verändert:
- D2(μ)⇒D2(μ)⋅wμ,wobeiw=e−j⋅2π/N.
- Jeder einzelne Butterfly im blau umrandeten Block (in der Grafikmitte) liefert durch Addition bzw. Subtraktion zwei Elemente der gesuchten Ausgangsfolge. Mit 0≤μ<N/2 gilt dabei:
- D(μ)=1/2⋅[D1(μ)+D2(μ)⋅wμ],
- D(μ+N/2)=1/2⋅[D1(μ)−D2(μ)⋅wμ].
Durch diese erste Anwendung des Überlagerungssatzes halbiert sich somit in etwa der Rechenaufwand.
Beispiel 1: Die DFT–Koeffizienten d(ν) zur Beschreibung des Zeitverlaufs seien entsprechend der Zeile 2 der folgenden Tabelle „dreieckförmig” belegt. Beachten Sie hierbei die periodische Fortsetzung der DFT, so dass der lineare Anstieg für t<0 durch die Koeffizienten d(8), ...,d(15) ausgedrückt wird.
Durch Anwendung des DFT–Algorithmus mit N=16 erhält man die in der Zeile 3 angegebenen Spektralkoeffizienten D(μ), die bei Vernachlässigung des Aliasingfehlers gleich D(μ)=4⋅si2(π⋅μ/2) wären. Man erkennt, dass sich der Aliasingfehler nur auf die ungeradzahligen Koeffizienten auswirkt (schraffierte Felder). Beispielsweise müsste D(1)=16/π2≈1.621≠1.642 sein.
Spaltet man die Gesamtfolge ⟨d(ν)⟩ in zwei Teilfolgen ⟨d1′(ν)⟩ und ⟨d2′(ν)⟩ auf, und zwar derart, dass die erste (gelb hinterlegte) Teilfolge nur geradzahlige Koeffizienten (ν=0,2, ...,N–2) und die zweite (grün hinterlegt) nur ungeradzahlige Koeffizienten (ν=1,3, ...,N–1) beinhalten und alle anderen zu Null gesetzt sind, so erhält man die zugehörigen Folgen im Spektralbereich:
- ⟨D1′(μ)⟩∙−−(N)−−∘⟨d1′(ν)⟩,
- ⟨D2′(μ)⟩∙−−(N)−−∘⟨d2′(ν)⟩.
In den gelb bzw. grün hinterlegten Zeilen 4 ...7 erkennt man:
- Wegen d(ν)=d1′(ν)+d2′(ν) gilt auch D(μ)=D1′(μ)+D2′(μ). Dies lässt sich zum Beispiel mit dem Additionstheorem linearer Systeme begründen.
- Die Periode der Folge ⟨D1′(μ)⟩ ist aufgrund des Nullsetzens eines jeden zweiten Zeitkoeffizienten nun N/2 im Gegensatz zur Periode N der Folge ⟨D(μ)⟩:
- D1′(μ+N/2)=D1′(μ).
- ⟨D2′(μ)⟩ beinhaltet zusätzlich einen Phasenfaktor (Verschiebung um einen Abtastwert), der einen Vorzeichenwechsel zweier um N/2 auseinanderliegender Koeffizienten bewirkt:
- D2′(μ+N/2)=−D2′(μ).
- Die Berechnung von ⟨D1′(μ)⟩ und ⟨D2′(μ)⟩ ist aber jeweils ebenso aufwändig wie die Bestimmung von ⟨D(μ)⟩, da ⟨d1′(ν)⟩ und ⟨d2′(ν)⟩ ebenfalls aus N Elementen bestehen, auch wenn einige Null sind.
Beispiel 2: Zur Fortsetzung des ersten Beispiels wird nun die bisherige Tabelle um die Zeilen 8 bis 12 erweitert.
Verzichtet man auf die Koeffizienten d1′(ν)=0 mit ungeraden sowie auf d2′(ν)=0 mit geraden Indizes, so kommt man zu den Teilfolgen ⟨d1(ν)⟩ und ⟨d2(ν)⟩ entsprechend den Zeilen 9 und 11 . Man erkennt:
- Die Zeitfolgen ⟨d1(ν)⟩ und ⟨d2(ν)⟩ weisen ebenso wie die dazugehörigen Spektralfolgen ⟨D1(μ)⟩ und ⟨D2(μ)⟩ nur noch die Dimension N/2 auf.
- Ein Vergleich der Zeilen 5, 7, 10 und 12 zeigt für 0≤μ<N/2 folgenden Zusammenhang:
- D1′(μ)=1/2⋅D1(μ),
- D2′(μ)=1/2⋅D2(μ)⋅wμ.
- Entsprechend ergibt sich für N/2≤μ<N:
- D1′(μ)=1/2⋅D1(μ−N/2),
- D2′(μ)=1/2⋅D2(μ−N/2)⋅wμ=−1/2⋅D2(μ−N/2)⋅wμ−N/2.
- Zum Beispiel erhält man mit N=16 ⇒ w=e–j⋅π/8 für die Indizes μ=1 bzw. μ=9:
- D1′(1)=1.708/2=0.854,D2′(1)=1/2⋅(1.456+j0.603)⋅e−j⋅π/8=0.788
- ⇒D(1)=D1′(1)+D2′(1)=1.642.
- D9′(1)=1.708/2=0.854,D2′(9)=−1/2⋅(1.456+j0.603)⋅e−j⋅π/8=−0.788
- ⇒D(9)=D1′(9)+D2′(9)=0.066.
Fazit:
- Durch diese erste Anwendung des Überlagerungssatzes halbiert sich nahezu der Rechenaufwand.
- Statt O=1920 benötigt man nur noch O=2·448+8⋅(4+2)+16⋅2=976 reelle Operationen.
- Der erste Summand berücksichtigt die beiden DFT–Berechnungen mit N/2=8.
- Der Rest wird für die acht komplexen Multiplikationen und die 16 komplexen Additionen bzw. Subtraktionen benötigt.
Radix-2-Algorithmus nach Cooley und Tukey
Ebenso wie andere FFT–Algorithmen baut das hier vorgestellte Verfahren [CT65][1] von James W. Cooley und John W. Tukey auf dem Überlagerungssatz der DFT auf. Es funktioniert nur dann, wenn die Stützstellenzahl N eine Zweierpotenz ist.
Die Grafik verdeutlicht den Algorithmus für N=8, wobei wieder die Transformation vom Zeit– in den Frequenzbereich dargestellt ist.
- Vor dem eigentlichen FFT-Algorithmus müssen zunächst die Eingangswerte d(0),...,d(N−1) im grauen Block „Bitumkehroperation” umsortiert werden.
- Die Berechnung erfolgt in log2N=3 Stufen, wobei in jeder Stufe N/2=4 gleiche Berechnungen mit verschiedenen μ–Werten
(= Exponent des komplexen Drehfaktors) ausgeführt werden. Eine solche Basisoperation bezeichnet man auch als Butterfly. - Jeder Butterfly berechnet aus zwei (im Allgemeinen komplexen) Eingangsgrößen A und B die beiden Ausgangsgrößen A+B⋅wμ sowie A–B⋅wμ entsprechend der folgenden Skizze.
Fazit: Die komplexen Spektralkoeffizienten D(0),...,D(N−1) erhält man am Ausgang der letzten Stufe nach Division durch N.
- Wie in der Aufgabe 5.5Z gezeigt wird, ergibt sich gegenüber der DFT eine deutlich kürzere Rechenzeit, zum Beispiel für N=1024 um mehr als den Faktor 150.
- Die inverse DFT zur Berechnung der Zeit– aus den Spektralkoeffizienten geschieht mit dem gleichen Algorithmus und nur geringfügigen Modifizierungen.
Beispiel 3:
Abschließend wird ein C–Programm
- fft(N, Re, Im)
gemäß dem oben beschriebenen Radix–2–Algorithmus angegeben:
- Beim Aufruf beinhalten die beiden Float–Arrays „Re” und „Im” die N Real– und Imaginärteile der komplexen Zeitkoeffizienten d(0), ... , d(N−1).
- In den gleichen Feldern „Re” und „Im” werden am Programmende die komplexen Koeffizienten D(0), ... , D(N−1) an das aufrufende Programm zurückgegeben.
- Aufgrund der „In–Place”–Programmierung reichen somit für diesen Algorithmus N komplexe Speicherplätze aus, allerdings nur, wenn zu Beginn die Eingangswerte umsortiert werden.
- Dies geschieht durch das Programm „bitumkehr”, wobei die Inhalte von Re(ν) und Im(ν) in die Elemente Re(κ) und Im(κ) eingetragen werden. Beispiel 4 verdeutlicht die Vorgehensweise.
Beispiel 4: Bitumkehroperation
- Der neue Index κ ergibt sich, wenn man den Index ν als Dualzahl schreibt und anschließend die log2N Bit in umgekehrter Reihenfolge darstellt.
- Zum Beispiel wird aus ν=3 der neue Index κ=6.
Aufgaben zum Kapitel
Aufgabe 5.5: Fast-Fouriertransformation
Aufgabe 5.5Z: Rechenaufwand für die FFT
Quellenverzeichnis
- ↑ Cooley, J.W.; Tukey, J.W.: An Algorithm for the Machine Calculation of Complex Fourier Series. In: Mathematics of Computation, Vol. 19, No. 90. (Apr., 1965), pp. 297-301.
