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 (\nu = 0, 2, \hspace{0.03cm}\text{ ...} \hspace{0.1cm}, N–2) und die zweite (grün hinterlegt) nur ungeradzahlige Koeffizienten (\nu = 1, 3, \hspace{0.03cm}\text{ ...} \hspace{0.1cm} , N–1) beinhalten und alle anderen zu Null gesetzt sind, so erhält man die zugehörigen Folgen im Spektralbereich:
- \langle \hspace{0.1cm}{D_1}'(\mu)\hspace{0.1cm} \rangle \hspace{0.2cm}\bullet\!\!-\!\!\!-(N)\!-\!\!\!-\!\!\hspace{0.05cm}\circ\, \hspace{0.2cm} \langle \hspace{0.1cm} {d_1}'(\nu) \hspace{0.1cm}\rangle ,
- \langle \hspace{0.1cm}{D_2}'(\mu) \hspace{0.1cm}\rangle \hspace{0.2cm}\bullet\!\!-\!\!\!-(N)\!-\!\!\!-\!\!\hspace{0.05cm}\circ\, \hspace{0.2cm} \langle\hspace{0.1cm} {d_2}'(\nu) \rangle \hspace{0.1cm}\hspace{0.05cm}.
In den gelb bzw. grün hinterlegten Zeilen 4\hspace{0.05cm}\text{ ...} \hspace{0.05cm}7 erkennt man:
- Wegen d(\nu) = {d_1}'(\nu) + {d_2}'(\nu) gilt auch D(\mu ) = {D_1}'(\mu ) + {D_2}'(\mu ). Dies lässt sich zum Beispiel mit dem Additionstheorem linearer Systeme begründen.
- Die Periode der Folge \langle \hspace{0.1cm}{D_1}'(\mu )\hspace{0.1cm}\rangle ist aufgrund des Nullsetzens eines jeden zweiten Zeitkoeffizienten nun N/2 im Gegensatz zur Periode N der Folge \langle \hspace{0.1cm} D(\mu )\hspace{0.1cm}\rangle:
- {D_1}'(\mu + {N}/{2}) ={D_1}'(\mu)\hspace{0.05cm}.
- \langle \hspace{0.1cm} {D_2}'(\mu )\hspace{0.1cm}\rangle beinhaltet zusätzlich einen Phasenfaktor (Verschiebung um einen Abtastwert), der einen Vorzeichenwechsel zweier um N/2 auseinanderliegender Koeffizienten bewirkt:
- {D_2}'(\mu + {N}/{2}) = - {D_2}'(\mu)\hspace{0.05cm}.
- Die Berechnung von \langle \hspace{0.1cm}{D_1}'(\mu )\hspace{0.1cm}\rangle und \langle \hspace{0.1cm} {D_2}'(\mu )\hspace{0.1cm}\rangle ist aber jeweils ebenso aufwändig wie die Bestimmung von \langle \hspace{0.1cm}D(\mu )\hspace{0.1cm}\rangle, da \langle \hspace{0.1cm}{d_1}'(\nu)\hspace{0.1cm}\rangle und \langle \hspace{0.1cm}{d_2}'(\nu)\hspace{0.1cm}\rangle ebenfalls aus N Elementen bestehen, auch wenn einige Null sind.
\text{Beispiel 2:} Zur Fortsetzung des ersten Beispiels wird nun die bisherige Tabelle um die Zeilen 8 bis 12 erweitert.
Verzichtet man auf die Koeffizienten {d_1}'(\nu) = 0 mit ungeraden sowie auf {d_2}'(\nu) = 0 mit geraden Indizes, so kommt man zu den Teilfolgen \langle \hspace{0.1cm}d_1(\nu)\hspace{0.1cm}\rangle und \langle \hspace{0.1cm}d_2(\nu)\hspace{0.1cm}\rangle entsprechend den Zeilen 9 und 11 . Man erkennt:
- Die Zeitfolgen \langle \hspace{0.1cm}{d_1}(\nu )\hspace{0.1cm}\rangle und \langle \hspace{0.1cm}{d_2}(\nu )\hspace{0.1cm}\rangle weisen ebenso wie die dazugehörigen Spektralfolgen \langle \hspace{0.1cm}{D_1}(\mu )\hspace{0.1cm}\rangle und \langle \hspace{0.1cm}{D_2}(\mu )\hspace{0.1cm}\rangle nur noch die Dimension N/2 auf.
- Ein Vergleich der Zeilen 5, 7, 10 und 12 zeigt für 0 \le \mu \lt N/2 folgenden Zusammenhang:
- {D_1}'(\mu) = {1}/{2}\cdot {D_1}(\mu)\hspace{0.05cm},
- {D_2}'(\mu) = {1}/{2}\cdot {D_2}(\mu)\cdot w^{\hspace{0.04cm}\mu}\hspace{0.05cm}.
- Entsprechend ergibt sich für N/2 \le \mu \lt N:
- {D_1}'(\mu) = {1}/{2}\cdot {D_1}(\mu - {N}/{2})\hspace{0.05cm},
- {D_2}'(\mu) = {1}/{2}\cdot {D_2}(\mu {-} {N}/{2})\cdot w^{\hspace{0.04cm}\mu} = { - } {1}/{2}\cdot {D_2}(\mu-N/2)\cdot w^{\hspace{0.04cm}\mu {-} N/2}\hspace{0.05cm}.
- Zum Beispiel erhält man mit N = 16 ⇒ w = {\rm e}^{ – {\rm j}\hspace{0.04cm} \cdot \hspace{0.04cm}\pi/8} für die Indizes \mu = 1 bzw. \mu = 9:
- {D_1}'(1) = {1.708}/{2} = 0.854,\hspace{0.8cm} {D_2}'(1) ={1}/{2}\cdot (1.456 + {\rm j} 0.603) \cdot {\rm e}^{ - {\rm j} \hspace{0.05cm}\cdot \hspace{0.05cm} \pi/8} = 0.788
- \Rightarrow D(1) = {D_1}'(1)+ {D_2}'(1)= 1.642 \hspace{0.05cm}.
- {D_9}'(1) = {1.708}/{2} = 0.854,\hspace{0.8cm} {D_2}'(9) = - {1}/{2}\cdot (1.456 + {\rm j} 0.603) \cdot {\rm e}^{ - {\rm j} \hspace{0.05cm}\cdot \hspace{0.05cm} \pi/8} = - 0.788
- \Rightarrow D(9) = {D_1}'(9)+ {D_2}'(9)= 0.066 \hspace{0.05cm}.
\text{Fazit:}
- Durch diese erste Anwendung des Überlagerungssatzes halbiert sich nahezu der Rechenaufwand.
- Statt \mathcal{O}= 1920 benötigt man nur noch \mathcal{O} = 2 · 448 + 8 \cdot (4+2) + 16 \cdot 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), \hspace{0.05cm}\text{...} \hspace{0.1cm}, d( N - 1) im grauen Block „Bitumkehroperation” umsortiert werden.
- Die Berechnung erfolgt in \text{log}_2 N = 3 Stufen, wobei in jeder Stufe N/2 = 4 gleiche Berechnungen mit verschiedenen \mu–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 \cdot w^{\mu} sowie A – B \cdot w^{\mu} entsprechend der folgenden Skizze.
\text{Fazit:} Die komplexen Spektralkoeffizienten D(0), \hspace{0.05cm}\text{...} \hspace{0.1cm}, 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.
\text{Beispiel 3:}
Abschließend wird ein C–Programm
- \text{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 {\rm Re}( \nu) und {\rm Im}( \nu) in die Elemente {\rm Re}( \kappa) und {\rm Im}( \kappa) eingetragen werden. \text{Beispiel 4} verdeutlicht die Vorgehensweise.
\text{Beispiel 4: Bitumkehroperation}
- Der neue Index \kappa ergibt sich, wenn man den Index \nu als Dualzahl schreibt und anschließend die \text{log}_2 \hspace{0.05cm} N Bit in umgekehrter Reihenfolge darstellt.
- Zum Beispiel wird aus \nu = 3 der neue Index \kappa = 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.
