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

 < >