% Except otherwise noted,
% this work is © 2009 Dmitri Fedorov, under CC-BY-SA license, version 3+.

\chapter{Fast Fourier transform}
Fast Fourier transform (FFT) is an efficient algorithm to compute the
discrete Fourier transform (DFT).

Computing DFT of $N$ points in the naive way, using the definition,
takes $O(N^2)$ arithmetic operations, while an FFT can compute the same
result in only $O(N\log N)$ operations. The difference in speed can be
substantial, especially for large data sets. This improvement made many
DFT-based algorithms practical.

Since the inverse of a DFT is also a DFT any FFT algorithm can be
used in for the inverse DFT as well.

The most well known FFT algorithms, like the Cooley-Tukey algorithm,
depend upon the factorization of $N$. However, there are FFTs with
$O(N\log N)$ complexity for all $N$, even for prime $N$.

\section{Discrete Fourier Transform}
For a set of complex numbers $x_n$, $n=0,\dots,N-1$, 
the DFT is defined as
a set of complex numbers $c_k$,
	\begin{equation}
c_k=
\sum_{n=0}^{N-1} x_n e^{-2\pi i\frac{nk}{N}}\;,
k=0,\dots,N-1\;.
	\end{equation}
The inverse DFT is given by
	\begin{equation}
x_n=\frac{1}{N}\sum_{k=0}^{N-1} c_k e^{+2\pi i\frac{nk}{N}}\;.
	\end{equation}

These transformations can be viewed as expansion of the vector $x_n$
in terms of the orthogonal basis of vectors $e^{2\pi i{kn\over N}}$,
\begin{equation}
\sum_{n=0}^{N-1}
\left(e^{2\pi i{kn\over N}}\right)
\left(e^{-2\pi i{k^\prime n\over N}}\right)
=N\delta_{kk^\prime} \;.
\end{equation}

The DFT represent the amplitude and phase of the different sinusoidal
components in the input data $x_n$.

The DFT is widely used in different fields, like spectral analysis,
data compression, solution of partial differential equations and others.

\subsection{Applications}
\subsubsection{Data compression}
Several lossy (that is, with certain loss of data) image and sound
compression methods employ DFT as an approximation for the Fourier
series. The signal is discretized and transformed, and then the Fourier
coefficients of high/low frequencies, which are assumed to be unnoticeable,
are discarded. The decompressor computes the inverse transform based on
this reduced number of Fourier coefficients.

\subsubsection{Partial differential equations}
Discrete Fourier transforms are often used to solve partial differential
equations, where the DFT is used as an approximation for the Fourier
series (which is recovered in the limit of infinite $N$). The
advantage of this approach is that it expands the signal in complex
exponentials $e^{inx}$, which are eigenfunctions of differentiation:
$\frac{d}{dx} e^{inx} = in e^{inx}$. Thus, in the Fourier representation,
differentiation is simply multiplication by $in$.

A linear differential equation with constant coefficients is
transformed into an easily solvable algebraic equation. One then uses
the inverse DFT to transform the result back into the ordinary spatial
representation. Such an approach is called a {\em spectral method}.

\subsubsection{Convolution and Deconvolution}
FFT can be used to efficiently compute convolutions of two sequences. A
convolution is the pairwise product of elements from two different
sequences, such as in multiplying two polynomials or multiplying two
long integers.

Another example comes from data acquisition processes where the detector
introduces certain (typically Gaussian) blurring to the sampled signal.
A reconstruction of the original signal can be obtained by deconvoluting
the acquired signal with the detector's blurring function.

\section{Cooley-Tukey algorithm}
In its simplest incarnation this algorithm re-expresses the DFT of
size $N=2M$ in terms of two DFTs of size $M$,
	\begin{eqnarray}
c_k&=& \sum_{n=0}^{N-1} x_n e^{-2\pi i\frac{nk}{N}} \nonumber\\
&=&\sum_{m=0}^{M-1}x_{2m}e^{-2\pi i{mk\over M}}
+e^{-2\pi i{k\over N}}
\sum_{m=0}^{M-1}x_{2m+1}e^{-2\pi i{mk\over M}} \nonumber\\
&=&\left\{\begin{array}{ll}
c^{(\mathrm{even})}_k+
e^{-2\pi i{k\over N}}c^{(\mathrm{odd})}_k &, k< M\\
c^{(\mathrm{even})}_{k-M}-
e^{-2\pi i{k-M\over N}}c^{(\mathrm{odd})}_{k-M} &, k\ge M\\
\end{array}\right. \;,
	\end{eqnarray}
where $c^{(\mathrm{even})}$ and $c^{(\mathrm{odd})}$ are the DFTs of
the even- and odd-numbered sub-sets of $x$.

This re-expression of a size-$N$ DFT as two size-$\frac{N}{2}$ DFTs is
sometimes called the Danielson-Lanczos lemma\index{Danielson-Lanczos
lemma}. The exponents $e^{-2\pi i{k\over N}}$ are called {\em twiddle
factors}.

The operation count by application of the lemma is reduced from the
original $N^2$ down to $2(N/2)^2+N/2=N^2/2+N/2<N^2$.

For $N=2^p$ Danielson-Lanczos lemma can be applied recursively until the
data sets are reduced to one datum each. The number of operations is then
reduced to $O(N\ln N)$ compared to the original $O(N^2)$. The established
library FFT routines, like FFTW and GSL, further reduce the operation
count (by a constant factor) using advanced programming techniques like
precomputing the twiddle factors, effective memory management and others.

\section{Multidimensional DFT}
For example, a two-dimensional set of data $x_{n_1n_2}$, $n_1=1\dots N_1$,
$n_2=1\dots N_2$ has the discrete Fourier transform
\begin{equation}
c_{k_1k_2} =
\sum_{n_1=0}^{N_1-1}
\sum_{n_2=0}^{N_2-1}
x_{n_1n_2}
e^{-2\pi i{n_1k_1\over N_1}}
e^{-2\pi i{n_2k_2\over N_2}}\;.
\end{equation}

\section{C implementation}
%\lstinputlisting
%{/usr/users/fedorov/public_html/numeric/c/fft/fft.c}
