Difference between revisions of "Signal Representation/Fast Fourier Transform (FFT)"

From LNTwww
 
(77 intermediate revisions by 9 users not shown)
Line 1: Line 1:
 
{{LastPage}}
 
{{LastPage}}
 
{{Header
 
{{Header
|Untermenü=Zeit- und frequenzdiskrete Signaldarstellung
+
|Untermenü=Time and Frequency-Discrete Signal Representation
|Vorherige Seite=Spektralanalyse
+
|Vorherige Seite=Spectrum Analysis
 
|Nächste Seite=
 
|Nächste Seite=
 
}}
 
}}
  
==Rechenaufwand von DFT bzw. IDFT==   
+
==Complexity of DFT and IDFT==   
+
<br>
Ein Nachteil der direkten Berechnung der (im Allgemeinen komplexen) DFT–Zahlenfolgen
+
A disadvantage of the direct calculation of the&nbsp; $($generally complex$)$&nbsp; DFT sequences
 +
 
 +
:$$\langle \hspace{0.03cm}D(\mu)\hspace{0.03cm}\rangle  \hspace{0.2cm}\bullet\!\!-\!\!\!-(N)\!-\!\!\!-\!\!\hspace{0.05cm}\circ\, \hspace{0.2cm} \langle \hspace{0.03cm}d(\nu)\hspace{0.03cm} \rangle$$
 
   
 
   
gemäß den in Kapitel 5.2 angegebenen Gleichungen ist der große Rechenaufwand. Zum Beispiel kann dieser Aufwand für die Diskrete Fouriertransformation (DFT)
+
according to the equations given in chapter&nbsp; [[Signal_Representation/Discrete_Fourier_Transform_(DFT)|&raquo;Discrete Fourier Transform&raquo;]]&nbsp;  $\rm (DFT)$&nbsp; is the large computational cost.
 +
 
 +
We consider as an example the DFT,&nbsp; i.e. the calculation of the&nbsp; $D(\mu)$&nbsp; coefficients from the&nbsp; $d(\nu)$&nbsp; coefficients:
 
   
 
   
wie folgt abgeschätzt werden:
+
:$$N \cdot D(\mu)  =  \sum_{\nu = 0 }^{N-1}
*Wir gehen davon aus, dass die Potenzen des komplexen Drehfaktors $w$ = exp(–j2 $\pi/N$) bereits in Real– und Imaginärteilform in einer Lookup–Tabelle vorliegen.
+
d(\nu) \cdot  {w}^{\hspace{0.03cm}\nu \hspace{0.03cm} \cdot \hspace{0.05cm}\mu}
*Zur Berechnung eines einzelnen Koeffizienten benötigt man dann $N$–1 komplexe Multiplikationen und ebenso viele komplexe Additionen.
+
=
*Jede komplexe Addition erfordert zwei reelle Additionen:
+
  d(0) \cdot w^{\hspace{0.03cm}0} + d(1) \cdot w^{\hspace{0.03cm}\mu}+ d(2) \cdot w^{\hspace{0.03cm}2\mu}+\hspace{0.05cm}\text{ ...} \hspace{0.05cm}+ d(N-1) \cdot w^{\hspace{0.03cm}(N-1)\cdot \mu}.$$
 +
 
 +
The computational effort required for this is to be estimated,&nbsp; assuming that the powers of the complex rotation factor&nbsp; $w = {\rm e}^{-{\rm j} \hspace{0.05cm}\cdot \hspace{0.05cm}2 \pi/N}$&nbsp; already exist in real and imaginary part form in a lookup table.&nbsp; To calculate a single coefficient,&nbsp; one then needs&nbsp; $N-1$&nbsp; complex multiplications and as many complex additions,&nbsp; observing:
 +
*Each complex addition requires two real additions:
 +
:$$(R_1 + {\rm j} \cdot I_1) + (R_2 + {\rm j} \cdot I_2) = (R_1 +
 +
R_2) + {\rm j} \cdot (I_1 + I_2)\hspace{0.05cm}.$$
 +
*Each complex multiplication requires four real multiplications and two real additions&nbsp; $($a subtraction is treated as an addition$)$:
 +
:$$(R_1 + {\rm j} \cdot I_1)  (R_2 + {\rm j} \cdot I_2) = (R_1 \cdot
 +
R_2 - I_1 \cdot I_2) + {\rm j} \cdot (R_1 \cdot I_2 + R_2 \cdot
 +
I_1)\hspace{0.05cm}.$$
 +
*Thus,&nbsp; the following number of real multiplications and the number of real additions are required to calculate all&nbsp; $N$&nbsp; coefficients in total:
 +
:$$M = 4 \cdot N \cdot (N-1),$$
 +
:$$A = 2 \cdot N \cdot
 +
(N-1)+2 \cdot N \cdot (N-1)=M \hspace{0.05cm}.$$
 +
*In today's computers,&nbsp; multiplications and additions/subtractions need about the same computing time.&nbsp; It is sufficient to consider the total number&nbsp; $\mathcal{O} = M + A$&nbsp; of all operations:
 +
:$$\mathcal{O} = 8 \cdot N \cdot (N-1) \approx 8 \cdot N^2\hspace{0.05cm}.$$
 +
 
 +
{{BlaueBox|TEXT=
 +
$\text{Conclusion:}$&nbsp;
 +
*For a&nbsp; Discrete Fourier Transform&nbsp; $\rm (DFT)$&nbsp; with&nbsp; $N = 1000$&nbsp; one already needs almost eight million arithmetic operations.&nbsp; The same applies to an IDFT.
 
   
 
   
*Jede komplexe Multiplikation erfordert vier reelle Multiplikationen und zwei reelle Additionen (eine Subtraktion wird wie eine Addition behandelt):
+
*With&nbsp; $N =16$&nbsp; still &nbsp;$1920$&nbsp; computational operations are required.
 +
 
 +
*If the parameter&nbsp; $N$&nbsp; is a power to the base&nbsp; $2$,&nbsp; more computationally efficient algorithms can be applied.&nbsp; The multitude of such methods known from the literature are summarized under the collective term&nbsp; &raquo;'''Fast Fourier Transform'''&laquo;&nbsp; &ndash; abbreviated&nbsp; $\text{FFT}$.&nbsp; All these methods are based on the&nbsp; &raquo;superposition theorem&laquo;&nbsp; of the DFT.}}
 
   
 
   
*Somit sind zur Berechnung aller N Koeffizienten insgesamt die folgende Anzahl M reeller Multiplikationen und die Anzahl A reeller Additionen erforderlich:
+
==Superposition theorem of the DFT== 
 +
<br>
 +
The graph illustrates the so-called&nbsp; &raquo;superposition theorem&laquo;&nbsp; of the DFT using the example of&raquo; $N = 16$.&nbsp; Shown here is the transition from the time domain to the spectral domain,&nbsp; i.e. the calculation of the spectral domain coefficients from the time domain coefficients: &nbsp;    $\langle \hspace{0.03cm}D(\mu)\hspace{0.03cm}\rangle  \hspace{0.2cm}\bullet\!\!-\!\!\!-(N)\!-\!\!\!-\!\!\hspace{0.05cm}\circ\, \hspace{0.2cm} \langle \hspace{0.03cm} d(\nu) \hspace{0.03cm}\rangle.$
 +
 
 +
[[File:EN_Sig_T_5_5_S2.png|right|frame|Superposition theorem of the DFT]]
 +
 
 +
The algorithm described thereby is characterized by the following steps:
 +
*The sequence&nbsp; $\langle \hspace{0.1cm}d(\nu)\hspace{0.1cm}\rangle$&nbsp; of length&nbsp; $N$&nbsp; is divided into two subsequences&nbsp; $\langle \hspace{0.03cm} d_1(\nu)\hspace{0.03cm}\rangle$&nbsp;  and&nbsp; $\langle \hspace{0.03cm} d_2(\nu)\hspace{0.03cm}\rangle$&nbsp; each of half length&nbsp; $($highlighted in yellow and green respectively in the graphic$)$.&nbsp; With&nbsp; $0 \le \nu \lt N/2$&nbsp; one thus obtains the sequence elements
 +
:$$d_1(\nu) = d(2\nu), $$
 +
:$$d_2(\nu) = d(2\nu+1)
 +
\hspace{0.05cm}.$$
 +
*The initial sequences&nbsp; $\langle \hspace{0.03cm}D_1(\mu )\hspace{0.03cm}\rangle$&nbsp; and&nbsp; $\langle \hspace{0.03cm}D_2(\mu )\hspace{0.03cm}\rangle$&nbsp; of the two sub-blocks result from this each by its own DFT,&nbsp; but now only with half length&nbsp; $N/2 = 8$:
 +
:$$\langle \hspace{0.03cm}D_1(\mu) \hspace{0.03cm}\rangle  \hspace{0.2cm}\bullet\!\!-\!\!\!-(N/2)\!-\!\!\!-\!\!\hspace{0.05cm}\circ\, \hspace{0.2cm} \langle \hspace{0.03cm}d_1(\nu) \hspace{0.03cm}\rangle , $$
 +
:$$ \langle \hspace{0.03cm}D_2(\mu)\hspace{0.03cm} \rangle \hspace{0.2cm}\bullet\!\!-\!\!\!-(N/2)\!-\!\!\!-\!\!\hspace{0.05cm}\circ\, \hspace{0.2cm} \langle \hspace{0.03cm}d_2(\nu) \hspace{0.03cm}\rangle \hspace{0.05cm}.$$
 +
*The initial values&nbsp; $\langle \hspace{0.03cm} D_2(\mu )\hspace{0.03cm}\rangle$&nbsp; of the lower (green) DFT $($with&nbsp; $0 \le \mu \lt N/2)$&nbsp; are then changed in the block outlined in red by complex rotation factors with respect to phase:
 +
:$$D_2(\mu) \hspace{0.3cm} \Rightarrow \hspace{0.3cm}D_2(\mu) \cdot w^{\hspace{0.04cm}\mu}, \hspace{0.2cm}{\rm with}\hspace{0.1cm}w =
 +
{\rm e}^{-{\rm j} \hspace{0.05cm}\cdot \hspace{0.05cm}2 \pi/N} \hspace{0.05cm}.$$
 +
*Each single&nbsp; &raquo;'''butterfly'''&laquo;&nbsp; in the blue bordered block&nbsp; $($in the middle of the graph$)$&nbsp; yields two elements of the searched sequence by addition or subtraction.&nbsp; With&nbsp; $0 \le \mu \lt N/2$&nbsp; applies:
 +
:$$D(\mu) =  {1}/{2}\cdot \big[D_1(\mu) + D_2(\mu) \cdot w^{\hspace{0.04cm}\mu}\big],$$
 +
:$$D(\mu +{N}/{2})  =  {1}/{2}\cdot \big[D_1(\mu) - D_2(\mu) \cdot w^{\hspace{0.04cm}\mu}\big]\hspace{0.05cm}.$$
 +
 
 +
'''This first application of the superposition theorem thus roughly halves the computational effort.'''
 +
 
 +
{{GraueBox|TEXT= 
 +
$\text{Example 1:}$&nbsp;
 +
Let the DFT coefficients&nbsp; $d(\nu)$&nbsp; for the description of the time course be&nbsp; "triangular"&nbsp; according to&nbsp; '''line 2'''&nbsp; of the following table.&nbsp; Note here the periodic continuation of the DFT,&nbsp; so that the linear increase for&nbsp; $t \lt 0$&nbsp; is given by the coefficients&nbsp; $d(8), \hspace{0.05cm}\text{ ...} \hspace{0.05cm}, d(15)$.
 +
 
 +
Applying the DFT algorithm with&nbsp; $N = 16$&nbsp; one obtains the spectral coefficients&nbsp; $D(\mu )$&nbsp; given in&nbsp; '''line 3'''&nbsp; which would be equal&nbsp; $D(\mu ) = 4 \cdot \text{sinc}^2(\mu/2)$&nbsp; if the aliasing error were neglected.&nbsp; We can see that the aliasing error only affects the odd coefficients&nbsp; $($shaded boxes$)$.&nbsp; For example,&nbsp; $D(1) = 16/ \pi^2 \approx 1.621\neq 1.642$&nbsp; should be.
 +
 
 +
[[File:Sig_T_5_5_S2b_Version2.png|right|frame|Result table for &nbsp;$\text{Example 1}$&nbsp; for the superposition theorem of the DFT&nbsp; $($Note: "Zeile" &nbsp; &rArr;&nbsp; "row"$)$]]
 +
 
 +
If we split the total sequence&nbsp; $\langle \hspace{0.03cm}d(\nu)\hspace{0.03cm}\rangle$&nbsp; into two subsequences  such that the first subsequence&nbsp; $\langle \hspace{0.03cm}{d_1}'(\nu)\hspace{0.03cm}\rangle$ &nbsp; &rArr; &nbsp; yellow marked has only even coefficients&nbsp; $(\nu = 0, 2, \hspace{0.03cm}\text{ ...} \hspace{0.1cm}, N–2)$&nbsp;  and the second subsequence&nbsp; $\langle \hspace{0.03cm}{d_2}'(\nu)\hspace{0.03cm}\rangle$&nbsp; &rArr; &nbsp; green marked contains only odd coefficients&nbsp; $(\nu = 1, 3, \hspace{0.03cm}\text{ ...} \hspace{0.1cm} , N-1)$&nbsp; and all others are set to zero.&nbsp; The corresponding sequences in the spectral domain are obtained:
 
   
 
   
*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:
+
:$$ \langle \hspace{0.03cm}{D_1}'(\mu)\hspace{0.03cm} \rangle  \hspace{0.2cm}\bullet\!\!-\!\!\!-(N)\!-\!\!\!-\!\!\hspace{0.05cm}\circ\, \hspace{0.2cm} \langle \hspace{0.03cm} {d_1}'(\nu) \hspace{0.03cm}\rangle , $$
+
:$$ \langle \hspace{0.03cm}{D_2}'(\mu) \hspace{0.03cm}\rangle \hspace{0.2cm}\bullet\!\!-\!\!\!-(N)\!-\!\!\!-\!\!\hspace{0.05cm}\circ\, \hspace{0.2cm} \langle\hspace{0.03cm} {d_2}'(\nu) \rangle \hspace{0.03cm}\hspace{0.05cm}.$$
*Daraus folgt: Man benötigt bereits für eine DFT (oder eine IDFT) mit $N$ = 1000 knapp acht Millionen Rechenoperationen. Für $N$ = 16 sind noch 1920 Rechenoperationen erforderlich.
+
 
 +
In the yellow or green lines&nbsp; $4\hspace{0.05cm}\text{ ...} \hspace{0.05cm}7$&nbsp; you can see:
 +
*Because of&nbsp; $d(\nu) = {d_1}'(\nu) + {d_2}'(\nu)$&nbsp; also holds&nbsp;
 +
:$$D(\mu ) = {D_1}'(\mu ) + {D_2}'(\mu ).$$
 +
:This can be justified,&nbsp; for example,&nbsp; with the&nbsp; [[Signal_Representation/Fourier_Transform_Theorems#Multiplication_with_a_factor_-_Addition_Theorem|&raquo;Addition Theorem of Linear Systems&laquo;]].
 +
*The period of the sequence&nbsp; $\langle \hspace{0.03cm}{D_1}'(\mu )\hspace{0.03cm}\rangle$&nbsp; due to the zeroing of every second time coefficient is now&nbsp; $N/2$&nbsp; unlike the period&nbsp; $N$&nbsp; of the sequence&nbsp; $\langle \hspace{0.03cm} D(\mu )\hspace{0.03cm}\rangle$:
 +
:$${D_1}'(\mu + {N}/{2}) ={D_1}'(\mu)\hspace{0.05cm}.$$
 +
* The sequence&nbsp; $\langle \hspace{0.03cm} {D_2}'(\mu )\hspace{0.03cm}\rangle$&nbsp; additionally contains a phase factor&nbsp; $($shift by one sample$)$&nbsp; which causes a sign change of two coefficients separated by&nbsp; $N/2$:
 +
:$${D_2}'(\mu + {N}/{2}) = - {D_2}'(\mu)\hspace{0.05cm}.$$
 +
*The calculation of&nbsp; $\langle \hspace{0.03cm}{D_1}'(\mu )\hspace{0.03cm}\rangle$&nbsp; and&nbsp; $\langle \hspace{0.03cm} {D_2}'(\mu )\hspace{0.1cm}\rangle$&nbsp; is,&nbsp; however,&nbsp; in each case as time-consuming as the determination of&nbsp; $\langle \hspace{0.03cm}D(\mu )\hspace{0.03cm}\rangle$,&nbsp; since&nbsp; $\langle \hspace{0.03cm}{d_1}'(\nu)\hspace{0.03cm}\rangle$&nbsp; and&nbsp; $\langle \hspace{0.03cm}{d_2}'(\nu)\hspace{0.03cm}\rangle$&nbsp;  also consist of&nbsp; $N$&nbsp; elements,&nbsp; even if half of them are zeros.}}
 +
 
 +
 
 +
{{GraueBox|TEXT=
 +
$\text{Example 2:}$&nbsp;
 +
To continue the first example, the previous table is now extended by the rows &nbsp;$8$&nbsp; to &nbsp;$12$&nbsp;.
 +
                     
 +
[[File:EN_Sig_T_5_5_S2c.png|right|frame|Result table for &nbsp;$\text{Example 2}$&nbsp; for the superposition theorem of the DFT]]
  
 +
Omitting the coefficients&nbsp; ${d_1}'(\nu) = 0$&nbsp; with odd indices and&nbsp; ${d_2}'(\nu) = 0$&nbsp; with even indices,&nbsp; we arrive at the subsequences&nbsp; $\langle \hspace{0.03cm}d_1(\nu)\hspace{0.03cm}\rangle$&nbsp; and&nbsp; $\langle \hspace{0.03cm}d_2(\nu)\hspace{0.03cm}\rangle$&nbsp; corresponding to lines &nbsp;$9$&nbsp; and &nbsp;$11$.&nbsp;
  
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.
+
You can see:
+
*The time sequences&nbsp; $\langle \hspace{0.03cm}{d_1}(\nu )\hspace{0.03cm}\rangle$&nbsp; and&nbsp; $\langle \hspace{0.03cm}{d_2}(\nu )\hspace{0.03cm}\rangle$&nbsp; exhibit as well as the corresponding spectral sequences&nbsp; $\langle \hspace{0.03cm}{D_1}(\mu )\hspace{0.03cm}\rangle$&nbsp; and&nbsp; $\langle \hspace{0.03cm}{D_2}(\mu )\hspace{0.03cm}\rangle$&nbsp; only have the dimension&nbsp; $(N/2)$.
==Überlagerungssatz der DFT== 
 
  
Die Grafik verdeutlicht den so genannten Überlagerungssatz der DFT am Beispiel N = 16.
+
*A comparison of the lines&nbsp; $5$,&nbsp; $7$,&nbsp; $10$&nbsp; and&nbsp; $12$&nbsp; shows the following relationship for&nbsp; $0 \le \mu \lt N/2$&nbsp;:
 +
:$${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}.$$
 +
*Correspondingly,&nbsp; for&nbsp; $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}$$
 +
:$$ \Rightarrow \hspace{0.3cm}{D_2}'(\mu) = { - } {1}/{2}\cdot {D_2}(\mu-N/2)\cdot w^{\hspace{0.04cm}\mu {-} N/2}\hspace{0.05cm}.$$
  
Der dadurch beschriebene Algorithmus ist durch folgende Schritte gekennzeichnet:
+
*For example,&nbsp; with&nbsp; $N = 16$ &nbsp; ⇒ &nbsp; $w = {\rm e}^{ - {\rm j}\hspace{0.04cm} \cdot \hspace{0.04cm}\pi/8}$&nbsp; for the indices&nbsp; $\mu = 1$&nbsp; respectively&nbsp; $\mu = 9$:&nbsp;
*Die Folge $\langle d(ν)\rangle$ der Länge $N$ wird in zwei Teilfolgen $\langle d_1(ν)\rangle$ und $\langle d_2(ν)\rangle$ jeweils halber Länge separiert (gelb bzw. grün hinterlegt). Mit 0 $ν$ < $N$/2 erhält man die Folgenelemente
+
:$${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}.$$}}
 
   
 
   
*Die Ausgangsfolgen $\langle D_1(\mu )\rangle$ und $\langle D_2(\mu )\rangle$ der beiden Teilblöcke ergeben sich daraus jeweils durch eine eigene DFT, aber nun nur noch mit halber Länge $N$/2 = 8:
+
 
 +
{{BlaueBox|TEXT=
 +
$\text{Conclusion:}$&nbsp; 
 +
*This first application of the superposition theorem almost halves the computational effort.
 
   
 
   
*Die Ausgangswerte $\langle D_2(\mu )\rangle$ der unteren (grünen) DFT (mit 0 ≤ $\mu$ < $N$/2) werden danach im rot umrandeten Block durch komplexe Drehfaktoren hinsichtlich Phasenlage verändert:
+
*Instead of&nbsp; $\mathcal{O}= 1920$&nbsp; one only needs &nbsp;$\mathcal{O} = 2 - 448 + 8 \cdot (4+2) + 16 \cdot 2 = 976$&nbsp; real operations.
 
   
 
   
*Jeder einzelne '''Butterfly''' im blau umrandeten Block liefert durch Addition bzw. Subtraktion zwei Elemente der gesuchten Ausgangsfolge. Mit 0 ≤ $\mu$ < $N$/2 gilt dabei:
+
*The first summand accounts for the two DFT calculations with&nbsp; $N/2 = 8$.
 
   
 
   
Durch diese erste Anwendung des Überlagerungssatzes halbiert sich in etwa der Rechenaufwand.
+
*The remainder is needed for the eight complex multiplications and the&nbsp; $16$&nbsp; complex additions and subtractions, respectively.}}
  
{{Beispiel}}
+
==Radix-2 algorithm according to Cooley and Tukey== 
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$(9) ... $d$(15) ausgedrückt wird.
+
<br>
Durch Anwendung des DFT–Algorithmus mit $N$ = 16 erhält man die in der ebenfalls blau hinterlegten Zeile 3 angegebenen Spektralkoeffizienten $D(\mu )$, die bei Vernachlässigung des Aliasingfehlers gleich $D(\mu ) = 4 \cdot \text{si}^2(\pi \cdot \mu/2)$ wären. Man erkennt, dass sich der Aliasingfehler nur auf die ungeradzahligen Koeffizienten auswirkt (schraffierte Felder). Beispielsweise müsste $D$(1) = 16/ $\pi^2 \approx$ 1.621 $\neq$ 1.642 sein.
+
Like other FFT algorithms,&nbsp; the method &nbsp; [CT65]<ref name ='CT65'>Cooley, J.W.; Tukey, J.W.:&nbsp; An Algorithm for the Machine Calculation of Complex Fourier Series.&nbsp; <br>In:&nbsp; Mathematics of Computation, Vol. 19, No. 90. (Apr., 1965), pp. 297-301.</ref>&nbsp;  from &nbsp; [https://en.wikipedia.org/wiki/James_Cooley &raquo;$\text{James W. Cooley}$&laquo;]&nbsp; and&nbsp; [https://en.wikipedia.org/wiki/John_Tukey &raquo;$\text{John W. Tukey}$&laquo;]&nbsp; is based on the superposition theorem of the DFT.&nbsp; It only works if the number of interpolation points is a power of two.  
  
Spaltet man die Gesamtfolge $\langle d(ν)\rangle$ in zwei Teilfolgen $\langle d_1'(ν)\rangle$ und $\langle d_2'(ν)\rangle$ 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 0 gesetzt sind, so erhält man die zugehörigen Folgen im Spektralbereich:
+
The diagram illustrates the algorithm for&nbsp; $N = 8$,&nbsp; again showing the transformation from the time to the frequency domain.
 +
[[File:EN_Sig_T_5_5_S3a.png|right|frame|Radix-2 algorithm&nbsp; $($flow diagram$)$;&nbsp; note:&nbsp; "level"&nbsp; and&nbsp; stage" are synonymous terms]]
 +
 
 +
*Before the FFT algorithm,&nbsp; the input values&nbsp; $d(0), \hspace{0.05cm}\text{...} \hspace{0.1cm}, d( N - 1)$&nbsp; have to be reordered in the grey block&nbsp; &raquo;Bit Reversal Operation&laquo;.
 
   
 
   
In den gelb bzw. grün hinterlegten Zeilen 4 ... 7 erkennt man:
+
*The computation is done in&nbsp; $\text{log}_2 N = 3$&nbsp; stages,&nbsp; where in each stage&nbsp; $N/2 = 4$&nbsp; equal computations are performed with different&nbsp; $\mu$&ndash;values as exponent of the complex rotation factor.&nbsp; Such a basic operation is called&nbsp; &raquo;'''butterfly'''&laquo;.
*Wegen $d(ν) = d_1'(ν) + d_2'(ν)$ 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 D_1'(\mu )\rangle$ beträgt aufgrund des Nullsetzens eines jeden zweiten Zeitkoeffizienten nun $N$/2 im Gegensatz zur Periode $N$ der Ursprungsfolge $\langle D(\mu )\rangle$:
+
*Each butterfly calculates from two&nbsp; $($generally complex$)$&nbsp; input variables&nbsp; $A$&nbsp; and&nbsp; $B$&nbsp; the two output variables&nbsp; $A + B \cdot w^{\mu}$&nbsp; and&nbsp; $A - B \cdot w^{\mu}$&nbsp; according to the following sketch.
+
 
* $\langle D_2'(\mu )\rangle$ beinhaltet zusätzlich einen Phasenfaktor (Verschiebung um einen Abtastwert), der einen Vorzeichenwechsel zweier um $N$/2 auseinanderliegender Koeffizienten bewirkt:
 
 
*Die Berechnung von $\langle D_1'(\mu )\rangle$ und $\langle D_2'(\mu )\rangle$ ist aber jeweils ebenso aufwändig wie die Bestimmung von $\langle D(\mu )\rangle$ , da $\langle d_'(ν)\rangle$ und $\langle d_2'(ν)\rangle$ ebenfalls aus $N$ Elementen bestehen, auch wenn einige 0 sind.                        
 
  
 +
[[File:P_ID1174__Sig_T_5_5_S3b_neu.png|center|frame|Butterfly of the DFT algorithm]]
 +
<br clear=all>
 +
{{BlaueBox|TEXT=
 +
$\text{Conclusion:}$&nbsp; 
 +
The complex spectral coefficients&nbsp; $D(0), \hspace{0.05cm}\text{...} \hspace{0.1cm}, D( N - 1)$&nbsp; are obtained at the output of the last stage after division by&nbsp; $N$.
 +
*As shown in&nbsp; [[Aufgaben:Exercise_5.5Z:_Complexity_of_the_FFT|$\text{Exercise 5.5Z}$]],&nbsp;  this results in a much shorter computation time compared to the DFT,&nbsp; for example for&nbsp; $N = 1024$&nbsp; by more than a factor $150$.
  
Verzichtet man auf die Koeffizienten $d_1'(ν)$ = 0 mit ungeraden sowie auf $d_2'(ν)$ = 0 mit geraden Indizes, so kommt man zu den Teilfolgen $\langle d_1(ν)\rangle$ und $\langle d_2(ν)\rangle$  ⇒  Zeilen 9 und 11. Man erkennt:
+
*The inverse DFT for calculating the time coefficients from the spectral coefficients is done with the same algorithm and only slight modifications.}}
*Die beiden Zeitfolgen 〈d1(ν)〉 und 〈d2(ν)〉 weisen damit 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 ≤ $\mu$ < $N$/2 folgenden Zusammenhang:
 
 
*Entsprechend erhält man für $N$/2 ≤ $\mu$ < $N$:
 
 
*Zum Beispiel erhält man mit $N$ = 16  ⇒  $w$ = exp(–j $\cdot \pi$ /8) für die Indizes $\mu$ = 1 bzw. $\mu$ = 9:
 
 
{{end}}
 
 
  
Durch diese erste Anwendung des Überlagerungssatzes halbiert sich nahezu der Rechenaufwand. Statt $\mathcal{O}$ = 1920 benötigt man nur $\mathcal{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 die acht komplexen Multiplikationen und die 16 komplexen Additionen bzw. Subtraktionen.
 
  
==Radix-2-Algorithmus nach Cooley und Tukey== 
+
{{GraueBox|TEXT=
+
$\text{Example 3:}$&nbsp; 
Ebenso wie andere FFT–Algorithmen baut das hier vorgestellte Verfahren von Cooley und Tukey auf dem Überlagerungssatz der DFT auf. Es funktioniert nur dann, wenn die Stützstellenzahl $N$ eine Zweierpotenz ist. Das folgende Bild verdeutlicht den Algorithmus – siehe auch [CT65] – für das Beispiel $N$ = 8, wobei die Transformation vom Zeit– in den Frequenzbereich dargestellt ist.
+
Finally,&nbsp; the C program&nbsp; $\text{fft(N, Re, Im)}$&nbsp; according to the Radix-2 algorithm described above is given:
 +
[[File:EN_Sig_Programm.png|right|frame|Radix-2 algorithm (C program)]]
  
Vor dem eigentlichen FFT-Algorithmus müssen die Eingangswerte $d$(0), ... , $d$( $N$–1) im grauen Block Bitumkehroperation umsortiert werden. Weiter erkennt man aus obiger Darstellung:
+
*The two float arrays&nbsp; &raquo;Re&laquo;&nbsp; and&nbsp; &raquo;Im&raquo;&nbsp; contain the&nbsp; $N$&nbsp; real and imaginary parts of the complex time coefficients&nbsp; $d(0)$, ... , $d( N - 1)$.
*Die Berechnung erfolgt in $\text{log}_2 N$ = 3 Stufen, wobei in jeder Stufe genau $N$/2 = 4 prinzipiell gleiche Berechnungen mit unterschiedlichem $\mu$ (= 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 folgender Skizze:
 
  
*Die komplexen Spektralkoeffizienten $D$(0), ... , $D$( $N$–1) erhält man am Ausgang der letzten Stufe nach Division durch $N$. Wie in Aufgabe Z5.5 gezeigt wird, ergibt sich gegenüber der DFT eine deutlich kürzere Rechenzeit – z.B. für $N$ = 1024 um mehr als den Faktor 150.
+
*In the same fields&nbsp; &raquo;Re&laquo;&nbsp; and&nbsp; &raquo;Im&raquo;&nbsp; the complex coefficients&nbsp; $D(0)$, ... , $D( N - 1)$&nbsp; are returned to the main program.
Die inverse DFT zur Berechnung der Zeitkoeffizienten aus den Spektralkoeffizienten lässt sich mit dem gleichen Algorithmus und nur geringfügigen Modifizierungen bewerkstelligen.
 
  
 +
*Due to the&nbsp; &raquo;in-place&nbsp; programming&laquo;, &nbsp; $N$&nbsp; complex storage elements are thus sufficient for this algorithm but only if the input values are reordered at the beginning.
  
Nachfolgend sehen Sie das C–Programm '''fft(N, Re, Im)''' gemäß dem oben beschriebenen Radix–2–Algorithmus:
+
*This is done by the program&nbsp; &raquo;bit-reversal&laquo;,&nbsp; where the contents of&nbsp; ${\rm Re}( \nu)$&nbsp; and&nbsp; ${\rm Im}( \nu)$&nbsp; are entered into the elements&nbsp; ${\rm Re}( \kappa)$&nbsp; and&nbsp; ${\rm Im}( \kappa)$.&nbsp; $\text{Example 4}$&nbsp; illustrates this procedure.}}
 +
<br clear=all>
  
*Beim Aufruf beinhalten die beiden Float–Arrays „Re” und „Im” die jeweils $N$ Real– und Imaginärteile der komplexen Zeitkoeffizienten $d$(0), ... , $d$( $N$ – 1).
+
{{GraueBox|TEXT=
*In den gleichen Feldern „Re” und „Im” werden die N komplexen Spektralkoeffizienten $D$(0), ... , $D$( $N$–1) am Programmende an das aufrufende Programm zurückgegeben.
+
[[File:EN_Sig_T_5_5_S3d_v2.png|100px|right|frame|Radix-2 algorithm $($Bit-reversal operation for&nbsp; $N = 8)$]]
*Aufgrund dieser „In–Place”–Programmierung reichen für diesen Algorithmus $N$ komplexe Speicherplätze aus, allerdings nur, wenn zu Beginn die Eingangswerte umsortiert werden.
+
$\text{Example 4: Bit-reversal operation}$&nbsp;
*Dies geschieht durch das Programm „bitumkehr”, wobei die Inhalte von Re( $\nu$) und Im( $\nu$) in die Elemente Re( $\kappa$) und Im( $\kappa$) eingetragen werden. Die folgende Tabelle gilt für $N$ = 8.
 
  
Schreibt man den Index $\nu$ als Dualzahl und stellt die $\text{log}_2 N$ Bits in umgekehrter Reihenfolge dar, so ergibt sich der neue Index $\kappa$. Beispielsweise wird so aus $\nu$ = 3 der neue Index $\kappa$ = 6.
+
#The new index&nbsp; $\kappa$&nbsp; is obtained by writing the index&nbsp; $\nu$&nbsp; as a dual number and then representing the&nbsp; $\text{log}_2 \hspace{0.05cm} N$&nbsp; bits in reverse order.<br><br>
 +
#For example,&nbsp; $\nu = 3$&nbsp; becomes the new index&nbsp; $\kappa = 6$.}}
  
  
==Aufgaben zu Kapitel 5.5== 
+
==Exercises for the chapter== 
 +
<br>
 +
[[Aufgaben:Exercise 5.5: Fast Fourier Transform|Exercise 5.5: Fast Fourier Transform]]
  
 +
[[Aufgaben:Exercise 5.5Z: Complexity of The FFT|Exercise 5.5Z: Complexity of the FFT]]
  
 +
==References==
  
  

Latest revision as of 15:38, 28 June 2023

Complexity of DFT and IDFT


A disadvantage of the direct calculation of the  $($generally complex$)$  DFT sequences

$$\langle \hspace{0.03cm}D(\mu)\hspace{0.03cm}\rangle \hspace{0.2cm}\bullet\!\!-\!\!\!-(N)\!-\!\!\!-\!\!\hspace{0.05cm}\circ\, \hspace{0.2cm} \langle \hspace{0.03cm}d(\nu)\hspace{0.03cm} \rangle$$

according to the equations given in chapter  »Discrete Fourier Transform»  $\rm (DFT)$  is the large computational cost.

We consider as an example the DFT,  i.e. the calculation of the  $D(\mu)$  coefficients from the  $d(\nu)$  coefficients:

$$N \cdot D(\mu) = \sum_{\nu = 0 }^{N-1} d(\nu) \cdot {w}^{\hspace{0.03cm}\nu \hspace{0.03cm} \cdot \hspace{0.05cm}\mu} = d(0) \cdot w^{\hspace{0.03cm}0} + d(1) \cdot w^{\hspace{0.03cm}\mu}+ d(2) \cdot w^{\hspace{0.03cm}2\mu}+\hspace{0.05cm}\text{ ...} \hspace{0.05cm}+ d(N-1) \cdot w^{\hspace{0.03cm}(N-1)\cdot \mu}.$$

The computational effort required for this is to be estimated,  assuming that the powers of the complex rotation factor  $w = {\rm e}^{-{\rm j} \hspace{0.05cm}\cdot \hspace{0.05cm}2 \pi/N}$  already exist in real and imaginary part form in a lookup table.  To calculate a single coefficient,  one then needs  $N-1$  complex multiplications and as many complex additions,  observing:

  • Each complex addition requires two real additions:
$$(R_1 + {\rm j} \cdot I_1) + (R_2 + {\rm j} \cdot I_2) = (R_1 + R_2) + {\rm j} \cdot (I_1 + I_2)\hspace{0.05cm}.$$
  • Each complex multiplication requires four real multiplications and two real additions  $($a subtraction is treated as an addition$)$:
$$(R_1 + {\rm j} \cdot I_1) (R_2 + {\rm j} \cdot I_2) = (R_1 \cdot R_2 - I_1 \cdot I_2) + {\rm j} \cdot (R_1 \cdot I_2 + R_2 \cdot I_1)\hspace{0.05cm}.$$
  • Thus,  the following number of real multiplications and the number of real additions are required to calculate all  $N$  coefficients in total:
$$M = 4 \cdot N \cdot (N-1),$$
$$A = 2 \cdot N \cdot (N-1)+2 \cdot N \cdot (N-1)=M \hspace{0.05cm}.$$
  • In today's computers,  multiplications and additions/subtractions need about the same computing time.  It is sufficient to consider the total number  $\mathcal{O} = M + A$  of all operations:
$$\mathcal{O} = 8 \cdot N \cdot (N-1) \approx 8 \cdot N^2\hspace{0.05cm}.$$

$\text{Conclusion:}$ 

  • For a  Discrete Fourier Transform  $\rm (DFT)$  with  $N = 1000$  one already needs almost eight million arithmetic operations.  The same applies to an IDFT.
  • With  $N =16$  still  $1920$  computational operations are required.
  • If the parameter  $N$  is a power to the base  $2$,  more computationally efficient algorithms can be applied.  The multitude of such methods known from the literature are summarized under the collective term  »Fast Fourier Transform«  – abbreviated  $\text{FFT}$.  All these methods are based on the  »superposition theorem«  of the DFT.

Superposition theorem of the DFT


The graph illustrates the so-called  »superposition theorem«  of the DFT using the example of» $N = 16$.  Shown here is the transition from the time domain to the spectral domain,  i.e. the calculation of the spectral domain coefficients from the time domain coefficients:   $\langle \hspace{0.03cm}D(\mu)\hspace{0.03cm}\rangle \hspace{0.2cm}\bullet\!\!-\!\!\!-(N)\!-\!\!\!-\!\!\hspace{0.05cm}\circ\, \hspace{0.2cm} \langle \hspace{0.03cm} d(\nu) \hspace{0.03cm}\rangle.$

Superposition theorem of the DFT

The algorithm described thereby is characterized by the following steps:

  • The sequence  $\langle \hspace{0.1cm}d(\nu)\hspace{0.1cm}\rangle$  of length  $N$  is divided into two subsequences  $\langle \hspace{0.03cm} d_1(\nu)\hspace{0.03cm}\rangle$  and  $\langle \hspace{0.03cm} d_2(\nu)\hspace{0.03cm}\rangle$  each of half length  $($highlighted in yellow and green respectively in the graphic$)$.  With  $0 \le \nu \lt N/2$  one thus obtains the sequence elements
$$d_1(\nu) = d(2\nu), $$
$$d_2(\nu) = d(2\nu+1) \hspace{0.05cm}.$$
  • The initial sequences  $\langle \hspace{0.03cm}D_1(\mu )\hspace{0.03cm}\rangle$  and  $\langle \hspace{0.03cm}D_2(\mu )\hspace{0.03cm}\rangle$  of the two sub-blocks result from this each by its own DFT,  but now only with half length  $N/2 = 8$:
$$\langle \hspace{0.03cm}D_1(\mu) \hspace{0.03cm}\rangle \hspace{0.2cm}\bullet\!\!-\!\!\!-(N/2)\!-\!\!\!-\!\!\hspace{0.05cm}\circ\, \hspace{0.2cm} \langle \hspace{0.03cm}d_1(\nu) \hspace{0.03cm}\rangle , $$
$$ \langle \hspace{0.03cm}D_2(\mu)\hspace{0.03cm} \rangle \hspace{0.2cm}\bullet\!\!-\!\!\!-(N/2)\!-\!\!\!-\!\!\hspace{0.05cm}\circ\, \hspace{0.2cm} \langle \hspace{0.03cm}d_2(\nu) \hspace{0.03cm}\rangle \hspace{0.05cm}.$$
  • The initial values  $\langle \hspace{0.03cm} D_2(\mu )\hspace{0.03cm}\rangle$  of the lower (green) DFT $($with  $0 \le \mu \lt N/2)$  are then changed in the block outlined in red by complex rotation factors with respect to phase:
$$D_2(\mu) \hspace{0.3cm} \Rightarrow \hspace{0.3cm}D_2(\mu) \cdot w^{\hspace{0.04cm}\mu}, \hspace{0.2cm}{\rm with}\hspace{0.1cm}w = {\rm e}^{-{\rm j} \hspace{0.05cm}\cdot \hspace{0.05cm}2 \pi/N} \hspace{0.05cm}.$$
  • Each single  »butterfly«  in the blue bordered block  $($in the middle of the graph$)$  yields two elements of the searched sequence by addition or subtraction.  With  $0 \le \mu \lt N/2$  applies:
$$D(\mu) = {1}/{2}\cdot \big[D_1(\mu) + D_2(\mu) \cdot w^{\hspace{0.04cm}\mu}\big],$$
$$D(\mu +{N}/{2}) = {1}/{2}\cdot \big[D_1(\mu) - D_2(\mu) \cdot w^{\hspace{0.04cm}\mu}\big]\hspace{0.05cm}.$$

This first application of the superposition theorem thus roughly halves the computational effort.

$\text{Example 1:}$  Let the DFT coefficients  $d(\nu)$  for the description of the time course be  "triangular"  according to  line 2  of the following table.  Note here the periodic continuation of the DFT,  so that the linear increase for  $t \lt 0$  is given by the coefficients  $d(8), \hspace{0.05cm}\text{ ...} \hspace{0.05cm}, d(15)$.

Applying the DFT algorithm with  $N = 16$  one obtains the spectral coefficients  $D(\mu )$  given in  line 3  which would be equal  $D(\mu ) = 4 \cdot \text{sinc}^2(\mu/2)$  if the aliasing error were neglected.  We can see that the aliasing error only affects the odd coefficients  $($shaded boxes$)$.  For example,  $D(1) = 16/ \pi^2 \approx 1.621\neq 1.642$  should be.

Result table for  $\text{Example 1}$  for the superposition theorem of the DFT  $($Note: "Zeile"   ⇒  "row"$)$

If we split the total sequence  $\langle \hspace{0.03cm}d(\nu)\hspace{0.03cm}\rangle$  into two subsequences such that the first subsequence  $\langle \hspace{0.03cm}{d_1}'(\nu)\hspace{0.03cm}\rangle$   ⇒   yellow marked has only even coefficients  $(\nu = 0, 2, \hspace{0.03cm}\text{ ...} \hspace{0.1cm}, N–2)$  and the second subsequence  $\langle \hspace{0.03cm}{d_2}'(\nu)\hspace{0.03cm}\rangle$  ⇒   green marked contains only odd coefficients  $(\nu = 1, 3, \hspace{0.03cm}\text{ ...} \hspace{0.1cm} , N-1)$  and all others are set to zero.  The corresponding sequences in the spectral domain are obtained:

$$ \langle \hspace{0.03cm}{D_1}'(\mu)\hspace{0.03cm} \rangle \hspace{0.2cm}\bullet\!\!-\!\!\!-(N)\!-\!\!\!-\!\!\hspace{0.05cm}\circ\, \hspace{0.2cm} \langle \hspace{0.03cm} {d_1}'(\nu) \hspace{0.03cm}\rangle , $$
$$ \langle \hspace{0.03cm}{D_2}'(\mu) \hspace{0.03cm}\rangle \hspace{0.2cm}\bullet\!\!-\!\!\!-(N)\!-\!\!\!-\!\!\hspace{0.05cm}\circ\, \hspace{0.2cm} \langle\hspace{0.03cm} {d_2}'(\nu) \rangle \hspace{0.03cm}\hspace{0.05cm}.$$

In the yellow or green lines  $4\hspace{0.05cm}\text{ ...} \hspace{0.05cm}7$  you can see:

  • Because of  $d(\nu) = {d_1}'(\nu) + {d_2}'(\nu)$  also holds 
$$D(\mu ) = {D_1}'(\mu ) + {D_2}'(\mu ).$$
This can be justified,  for example,  with the  »Addition Theorem of Linear Systems«.
  • The period of the sequence  $\langle \hspace{0.03cm}{D_1}'(\mu )\hspace{0.03cm}\rangle$  due to the zeroing of every second time coefficient is now  $N/2$  unlike the period  $N$  of the sequence  $\langle \hspace{0.03cm} D(\mu )\hspace{0.03cm}\rangle$:
$${D_1}'(\mu + {N}/{2}) ={D_1}'(\mu)\hspace{0.05cm}.$$
  • The sequence  $\langle \hspace{0.03cm} {D_2}'(\mu )\hspace{0.03cm}\rangle$  additionally contains a phase factor  $($shift by one sample$)$  which causes a sign change of two coefficients separated by  $N/2$:
$${D_2}'(\mu + {N}/{2}) = - {D_2}'(\mu)\hspace{0.05cm}.$$
  • The calculation of  $\langle \hspace{0.03cm}{D_1}'(\mu )\hspace{0.03cm}\rangle$  and  $\langle \hspace{0.03cm} {D_2}'(\mu )\hspace{0.1cm}\rangle$  is,  however,  in each case as time-consuming as the determination of  $\langle \hspace{0.03cm}D(\mu )\hspace{0.03cm}\rangle$,  since  $\langle \hspace{0.03cm}{d_1}'(\nu)\hspace{0.03cm}\rangle$  and  $\langle \hspace{0.03cm}{d_2}'(\nu)\hspace{0.03cm}\rangle$  also consist of  $N$  elements,  even if half of them are zeros.


$\text{Example 2:}$  To continue the first example, the previous table is now extended by the rows  $8$  to  $12$ .

Result table for  $\text{Example 2}$  for the superposition theorem of the DFT

Omitting the coefficients  ${d_1}'(\nu) = 0$  with odd indices and  ${d_2}'(\nu) = 0$  with even indices,  we arrive at the subsequences  $\langle \hspace{0.03cm}d_1(\nu)\hspace{0.03cm}\rangle$  and  $\langle \hspace{0.03cm}d_2(\nu)\hspace{0.03cm}\rangle$  corresponding to lines  $9$  and  $11$. 

You can see:

  • The time sequences  $\langle \hspace{0.03cm}{d_1}(\nu )\hspace{0.03cm}\rangle$  and  $\langle \hspace{0.03cm}{d_2}(\nu )\hspace{0.03cm}\rangle$  exhibit as well as the corresponding spectral sequences  $\langle \hspace{0.03cm}{D_1}(\mu )\hspace{0.03cm}\rangle$  and  $\langle \hspace{0.03cm}{D_2}(\mu )\hspace{0.03cm}\rangle$  only have the dimension  $(N/2)$.
  • A comparison of the lines  $5$,  $7$,  $10$  and  $12$  shows the following relationship for  $0 \le \mu \lt N/2$ :
$${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}.$$
  • Correspondingly,  for  $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}$$
$$ \Rightarrow \hspace{0.3cm}{D_2}'(\mu) = { - } {1}/{2}\cdot {D_2}(\mu-N/2)\cdot w^{\hspace{0.04cm}\mu {-} N/2}\hspace{0.05cm}.$$
  • For example,  with  $N = 16$   ⇒   $w = {\rm e}^{ - {\rm j}\hspace{0.04cm} \cdot \hspace{0.04cm}\pi/8}$  for the indices  $\mu = 1$  respectively  $\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{Conclusion:}$ 

  • This first application of the superposition theorem almost halves the computational effort.
  • Instead of  $\mathcal{O}= 1920$  one only needs  $\mathcal{O} = 2 - 448 + 8 \cdot (4+2) + 16 \cdot 2 = 976$  real operations.
  • The first summand accounts for the two DFT calculations with  $N/2 = 8$.
  • The remainder is needed for the eight complex multiplications and the  $16$  complex additions and subtractions, respectively.

Radix-2 algorithm according to Cooley and Tukey


Like other FFT algorithms,  the method   [CT65][1]  from   »$\text{James W. Cooley}$«  and  »$\text{John W. Tukey}$«  is based on the superposition theorem of the DFT.  It only works if the number of interpolation points is a power of two.

The diagram illustrates the algorithm for  $N = 8$,  again showing the transformation from the time to the frequency domain.

Radix-2 algorithm  $($flow diagram$)$;  note:  "level"  and  stage" are synonymous terms
  • Before the FFT algorithm,  the input values  $d(0), \hspace{0.05cm}\text{...} \hspace{0.1cm}, d( N - 1)$  have to be reordered in the grey block  »Bit Reversal Operation«.
  • The computation is done in  $\text{log}_2 N = 3$  stages,  where in each stage  $N/2 = 4$  equal computations are performed with different  $\mu$–values as exponent of the complex rotation factor.  Such a basic operation is called  »butterfly«.
  • Each butterfly calculates from two  $($generally complex$)$  input variables  $A$  and  $B$  the two output variables  $A + B \cdot w^{\mu}$  and  $A - B \cdot w^{\mu}$  according to the following sketch.


Butterfly of the DFT algorithm


$\text{Conclusion:}$  The complex spectral coefficients  $D(0), \hspace{0.05cm}\text{...} \hspace{0.1cm}, D( N - 1)$  are obtained at the output of the last stage after division by  $N$.

  • As shown in  $\text{Exercise 5.5Z}$,  this results in a much shorter computation time compared to the DFT,  for example for  $N = 1024$  by more than a factor $150$.
  • The inverse DFT for calculating the time coefficients from the spectral coefficients is done with the same algorithm and only slight modifications.


$\text{Example 3:}$  Finally,  the C program  $\text{fft(N, Re, Im)}$  according to the Radix-2 algorithm described above is given:

Radix-2 algorithm (C program)
  • The two float arrays  »Re«  and  »Im»  contain the  $N$  real and imaginary parts of the complex time coefficients  $d(0)$, ... , $d( N - 1)$.
  • In the same fields  »Re«  and  »Im»  the complex coefficients  $D(0)$, ... , $D( N - 1)$  are returned to the main program.
  • Due to the  »in-place  programming«,   $N$  complex storage elements are thus sufficient for this algorithm but only if the input values are reordered at the beginning.
  • This is done by the program  »bit-reversal«,  where the contents of  ${\rm Re}( \nu)$  and  ${\rm Im}( \nu)$  are entered into the elements  ${\rm Re}( \kappa)$  and  ${\rm Im}( \kappa)$.  $\text{Example 4}$  illustrates this procedure.


Radix-2 algorithm $($Bit-reversal operation for  $N = 8)$

$\text{Example 4: Bit-reversal operation}$ 

  1. The new index  $\kappa$  is obtained by writing the index  $\nu$  as a dual number and then representing the  $\text{log}_2 \hspace{0.05cm} N$  bits in reverse order.

  2. For example,  $\nu = 3$  becomes the new index  $\kappa = 6$.


Exercises for the chapter


Exercise 5.5: Fast Fourier Transform

Exercise 5.5Z: Complexity of the FFT

References

  1. 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.