Franz J. Vesely > CompPhys Tutorial > Appendix  
 
 





 
< >
Appendix A2



Appendix



A2 Discrete Fourier Transformation


A2.1 Fundamentals


We are using the convention

$ \tilde{f}(\nu)=\int \limits_{-\infty}^{\infty}f(t)\,e^{2\pi\,i \nu t}\, dt\,,\;\;\; f(t)=\int \limits_{-\infty}^{\infty}\tilde{f}(\nu)\,e^{-2\pi\,i \nu t}\, d\nu $     (8.94)

Assume that the function $f(t)$ is given only at discrete, equidistant values of its argument:

$ f_{k}\equiv f(t_{k})\equiv f(k\, \Delta t) \;\;\;k=\dots -2, -1,0,1,2, \dots $     (8.95)

The reciprocal value of the time increment $\Delta t$ is called sampling rate. The higher the sampling rate, the more details of the given function $f(t)$ will be captured by the table of discrete values $f_{k}$. This intuitively evident fact is put in quantitative terms by Nyquist's theorem: if the Fourier spectrum of $f(t)$,

$ \tilde{f}(\nu)\equiv \int \limits_{-\infty}^{\infty}f(t) e^{\textstyle 2\pi i\nu t} dt $      (8.96)

is negligible for frequencies beyond the critical (or Nyquist) frequency

$ \pm \nu_{0} \equiv \pm \frac{\textstyle 1}{\textstyle 2 \Delta t} $     (8.97)

then $f(t)$ is called a band-limited process. Such a process is completely determined by its sampled values $f_{k}$. The formula that permits the reconstruction of $f(t)$ from the sampled data reads

$ f(t)=\sum \limits_{k=-\infty}^{\infty} f_{k} \frac{\textstyle \sin[2\pi \nu_{0}(t-k\Delta t)]}{\textstyle 2\pi \nu_{0}(t-k\Delta t)} $     (8.98)

(In contrast, if $f(t)$ is not band-limited, sampling with finite time resolution results in "mirroring in" the outlying parts of the spectrum from beyond $\pm \nu_{0}$, superposing them on the correct spectrum. In signal processing this effect is known as "aliasing".)

Let us assume now that a finite set of sampled values is given:

$ f_{k}\,,\;\;k=0,1,\,\dots\,N-1 $     (8.99)

and let $N$ be an even number. Define discrete frequencies by

$ \nu_{n}\equiv\frac{\textstyle n}{\textstyle N\Delta t}\,,\;\;\; n=-\frac{\textstyle N}{\textstyle 2},\dots\,,0,\,\dots ,\frac{\textstyle N}{\textstyle 2} $     (8.100)

(The $\nu_{n}$ pertaining to $n=N/2$ is again the Nyquist frequency.) Then the Fourier transform of $f(t)$ at some frequency $\nu_{n}$ is given by

$ \tilde{f}(\nu_{n}) \approx \Delta t \sum \limits_{k=0}^{N-1}f_{k} e^{\textstyle 2\pi i \nu_{n}t_{k}} = \Delta t \sum \limits_{k=0}^{N-1}f_{k}e^{\textstyle 2\pi i kn/N} $      (8.101)

Thus it makes sense to define the discrete Fourier transform as

$ F_{n} \equiv \sum \limits_{k=0}^{N-1} f_{k}e^{\textstyle 2\pi ikn/N} \;\;{\rm with}\; N \; {\rm even,\;\; and}\;\; n=0, \pm 1, \, \dots, \, N/2 $

According to 8.101 the Fourier transform proper is just $\tilde{f}(\nu_{n}) \approx \Delta t \, F_{n}$.

From the definition of $F_{n}$ it follows that $F_{-n}=F_{N-n}$. We make use of this periodicity to renumber the $F_{n}$ such that $n$ runs from $0$ to $N-1$ (instead of $ -N/2 $ to $N/2$):

$ 0 $, $ 1 $, $ \dots $ $ \frac{\textstyle N}{\textstyle 2}-1$, $ \frac{\textstyle N}{\textstyle 2}$ $ [ - \frac{\textstyle N}{\textstyle 2}$, $-\frac{\textstyle N}{\textstyle 2}+1$, $ \dots $ $ -1, ] $
$\Longrightarrow$
$ 0 $, $ 1 $, $ \dots $ $ \frac{\textstyle N}{\textstyle 2}-1 $, $ \frac{\textstyle N}{\textstyle 2}$, $ [\equiv \frac{\textstyle N}{\textstyle 2}$, $ \frac{\textstyle N}{\textstyle 2}+1$, $ \dots$ $N-1 ]$


With this indexing convention the back transformation may be conveniently written

$ f_{k} =\frac{\textstyle 1}{\textstyle N} \sum \limits_{n=0}^{N-1}F_{n}e^{\textstyle -2\pi ikn/N} $      (8.102)


Fast Fourier Transform (FFT) If we were to use the definition 8.102 "as is" to calculate the discrete Fourier transform, we would have to perform some $N^{2}$ operations. Cooley and Tukey (and before them Danielson and Lanczos; see PRESS 86]) have demonstrated how, by smart handling of data, the number of operations may be pushed down to $\approx N\, \log_{2}N$. Note that for $N=1000$ this is an acceleration of $100:1$. Indeed, many algorithms of modern computational physics hinge on this possibility of rapidly transforming back and forth long tables of function values.

In the following it is always assumed that $N=2^{m}$. If $N$ is not a power of $2$, simply "pad" the table, putting $f_{k}=0$ up to the next useful table length. Defining

$ W_{N}\equiv e^{\textstyle 2\pi i/N} $     (8.103)

we realize that $W_{N}^{2}=W_{N/2}$ etc. The discrete Fourier transform is therefore

$ \begin{eqnarray} F_{N}&=&\sum \limits_{k=0}^{N-1}W_{N}^{nk}f_{k} \\ &=&\sum \limits_{l=0}^{N/2-1}W_{N/2}^{nl}f_{2l} +W_{N}^{n}\sum \limits_{l=0}^{N/2-1}W_{N/2}^{nl}f_{2l+1} \\ &\equiv&\quad F_{n}^{e}\quad+\quad W_{N}^{n} F_{n}^{o} \end{eqnarray} $      (8.104-8.106)


where the indices e and o stand for "even" and "odd". Next we treat each of the two terms to the right of 8.106 by the same pattern, finding

$ \begin{eqnarray} F_{n}^{e}&=&F_{n}^{ee}+W_{N/2}^{n}F_{n}^{eo} \\ F_{n}^{o}&=&F_{n}^{oe}+W_{N/2}^{n}F_{n}^{oo} \end{eqnarray} $     (8.107-8.108)


By iterating this procedure $m=\log_{2} N$ times we finally arrive at terms $F_{n}^{(...)}$ that are identical to the given table values $f_{k}$.


EXAMPLE: Putting $N=4$ we have $W_{4} \equiv exp[2\pi i/4]$ and

$ \begin{eqnarray} F_{n}&=&\sum \limits_{k=0}^{3}W_{4}^{nk}f_{k} \;\;n=0,\dots 3 \\ &=&\sum \limits_{l=0}^{1}W_{2}^{nl}f_{2l}+W_{4}^{n}\sum \limits_{l=0}^{1}W_{2}^{nl}f_{2l+1} \\ &\equiv& \;\; F_{n}^{e}\;\;\;+\;\;\;W_{4}^{n}F_{n}^{o} \\ &=&F_{n}^{ee}+W_{2}^{n}F_{n}^{eo} +W_{4}^{n}\left[ F_{n}^{oe}+W_{2}^{n}F_{n}^{oo}\right] \\ &=&f_{0}+W_{2}^{n}f_{2} +W_{4}^{n} \left[f_{1}+W_{2}^{n}f_{3} \right] \end{eqnarray} $      (8.109-8.113)


Thus the correspondence between the table values $f_{k}$ and the terms $F_{n}^{ee}$ etc. is as follows:

$ ee $ $ eo $ $ oe $ $ oo $
$ 0 $ $ 2 $ $ 1 $ $ 3 $


EXERCISE: Demonstrate that a similar analysis as above leads for $N=8$ to the correspondences

$ eee $ $ eeo $ $ eoe $ $ eoo $ $ oee $ $ oeo $ $ ooe $ $ ooo $
$ 0 $ $ 4 $ $ 2 $ $ 6 $ $ 1 $ $ 5 $ $ 3 $ $ 7 $


It is easy to see that this correspondence is reproduced by the following rule:
  1. put $e \leftrightarrow 0$ and $o \leftrightarrow 1$, such that $eeo \leftrightarrow 001$ etc.;

  2. reverse the bit pattern thus obtained and interpret the result as an integer number: $eeo \leftrightarrow 001 \leftrightarrow 100 \,=\,4$.


In other words, arrange the table values $f_{k}$ in bit-reversed order. (For example, $k=4$ is at position $1$ since $4=100 \longrightarrow 001=1$.)

The correctly arranged $f_{k}$ are now combined in pairs according to 8.114. The rule to follow in performing this "decimation" step is sketched in Fig. 8.4.

  $ 0 $ $ 4 $ $ 2 $ $ 6 $ $ 1 $ $ 5 $ $ 3 $ $ 7 $
$ m=1 $ $ a $ $ b $ $ a $ $ b $ $ a $ $ b $ $ a $ $ b $
$ m=2 $ $ a $ $ b $ $ a $ $ b $
$ m=3 $ $ a $ $ b $

Figure 8.4: Decimation for $N=8$


On each level ($m$) the terms $a,b$ are combined according to

$ a+W_{2^{m}}^{n}b $     (8.114)

It is evident that the number of operations is of order $N \, \log_{2}N$.

Further details of the method, plus sample programs in Pascal, Fortran, or C are given in PRESS 86.


EXERCISE: Sketch the pattern of Figure 8.4 for $N=4$ and perform the "decimation". Compare your result to equ. 8.113.



vesely 2006

< >