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:
-
put $e \leftrightarrow 0$ and $o \leftrightarrow 1$, such that
$eeo \leftrightarrow 001$ etc.;
-
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