\documentclass[12pt]{article}
\usepackage{fullpage,amsmath,amsfonts,amsthm,html}
%\usepackage[bottomafter]{draftcopy}

\newcommand{\R}{\mathbb{R}}
\newcommand{\T}{\mathbb{T}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\C}{\mathbb{C}}

\def\e{\mathrm{e}}
\def\i{\mathrm{i}}
\def\d{\mathrm{d}}
\def\per{{\operatorname{per}}}

\newtheorem{theorem}{Theorem}
\theoremstyle{definition}
\newtheorem{definition}{Definition}
\newtheorem{remark}{Remark}
\newtheorem{problem}{Problem}
\newtheorem{example}{Example}
\newtheorem{exercise}{Exercise}

\begin{document}
\title{The discrete and fast Fourier transforms}
\author{Marcel Oliver}
\date{Revised April 27, 2012}
\maketitle

\section{Fourier series}

We begin by recalling the familiar definition of the Fourier series.  For a
periodic function $u \colon [0,2 \pi] \to \C$, we define the \emph{Fourier
transform} 
\begin{equation}
  \hat u_k 
  = \frac1{2\pi} \int_0^{2\pi} \e^{-\i kx} \, u(x) \, \d x \,.
  \label{e.ft}
\end{equation}
for every $k \in \Z$.  By direct computation, noting that the boundary
terms of the integral cancel due to periodicity, we obtain the
\emph{orthogonality relation}
\begin{equation}
  \frac1{2\pi} \int_0^{2\pi} \e^{-\i jx} \, \e^{\i kx} \, \d x
  = \delta_{ij} \,,
  \label{e.orth}
\end{equation}
where $\delta_{jk}$ denotes the \emph{Kronecker symbol}
\begin{equation}
  \delta_{jk} =
  \begin{cases}
    1 & \text{if } j=k \,, \\ 
    0 & \text{otherwise} \,.
  \end{cases}
\end{equation}
We conclude that the \emph{inverse Fourier transform} is given by
\begin{equation}
  u(x) = \sum_{k \in \Z} \hat u_k \, \e^{\i kx} \,, 
  \label{e.ift}
\end{equation}
which can be verified by inserting \eqref{e.ift} into \eqref{e.ft},
changing the order of integration and summation, and applying the
orthogonality relation \eqref{e.orth}.  In the same way, we obtain
the \emph{Parseval identity}
\begin{equation}
  \frac1{2\pi} \int_0^{2\pi} \overline{u(x)} \, v(x) \, \d x
  = \sum_{k \in \Z} \overline{\hat u_k} \, \hat v_k \,,
  \label{e.parseval} 
\end{equation}
where the overbar denotes the complex conjugate.

\begin{remark}
The choice of the interval of definition of the periodic function $u$
is arbitrary, but choosing $[0,2\pi[$ has notational advantages.
Other intervals can be obtained by translation and scaling.  The
prefactor on the right of \eqref{e.ft} is sometimes fully or partially
moved into the inverse Fourier transform \eqref{e.ift}.
\end{remark}

\begin{remark}
More elegantly, we may speak of functions on the torus $\T =
\R/2\pi\Z$.  The right hand notation means that two real numbers are
considered equivalent if they differ by an integer multiple of $2\pi$.
\end{remark}

\begin{remark}
The precise functional setting is not important for our purposes, and
we will not dwell on this issue.  Let us remark that this theory is
developed by noting that (a) the inversion formula \eqref{e.ift} is a
continuous linear map from the sequence space $\ell^1$ into the space
of essentially bounded functions $L^\infty(\T)$ and that (b) the
Parseval identity \eqref{e.parseval} shows that the Fourier transform
is an isometry with respect to the norms of the larger respective
spaces $\ell^2$ and $L^2(\T)$, and can thus be extended to those
larger spaces using a density argument.  The details depend on
advanced Analysis and may be found, for example, in \cite{Folland}.
\end{remark}

\begin{remark}
The orthogonality relation \eqref{e.orth} expresses that $\{ \e^{\i
kx} \colon k \in \Z \}$ is an orthonormal set with respect to the
inner product
\begin{equation}
  \langle u, v \rangle
  = \frac1{2\pi} \int_0^{2\pi} \overline{u(x)} \, v(x) \, \d x \,.
\end{equation}
It is a nontrivial result that this set is a \emph{Hilbert basis} of
the function space $L^2(\T)$, which means that every $u \in L^2(\T)$
can be expressed as a countable linear combination of basis functions
\cite{Folland}.  From this point of view, the inverse Fourier
transform \eqref{e.ift} is nothing but the representation of a
function in terms of a Hilbert basis.
\end{remark}


\section{The discrete Fourier transform}

Let $v_0, \dots, v_{N-1}$ be an $N$-tuple of real or complex numbers.
We define the \emph{discrete Fourier transform} (DFT) to be the linear
map
\begin{equation}
  \tilde v_k = \frac1N \sum_{j=0}^{N-1} \e^{-\i k j h} \, v_j
  \label{e.dft}
\end{equation}
where $k=0, \dots, N-1$ and
\begin{equation}
  h = \frac{2\pi}N \,.
\end{equation}

\begin{remark}
We can think of $v_0, \dots, v_{N-1}$ as the samples of a
$2\pi$-periodic function $v$ on the equidistantly spaced nodes
\begin{equation}
  x_j = jh \,.
\end{equation}
In this case, we can write the formula for the DFT as
\begin{equation}
  \tilde v_k = \frac1N \sum_{j=0}^{N-1} \e^{-\i k x_j} \, v(x_j) \,,
\end{equation}
which looks like a Riemann sum approximation to the continuous Fourier
transform \eqref{e.ft}.  This interpretation, however, has some
subtleties which will be discussed further in Section~\ref{s.approx}.
\end{remark}

We now investigate the inverse of the discrete Fourier transform
\eqref{e.dft}.  To this end, we first derive the \emph{discrete
orthogonality relation}
\begin{align}
  \frac1N \sum_{j=0}^{N-1} \e^{-\i k jh} \, \e^{\i l jh}
  & = \frac1N \sum_{j=0}^{N-1} \e^{\i (l-k) jh} 
    = \frac1N \sum_{j=0}^{N-1}
      \biggl(
        \e^{\i 2 \pi \tfrac{l-k}N}_{\vphantom 1}
      \biggr)^j
      \notag \\
  & = \begin{cases}
        1 & \text{if } l=k +m N \text{ for some } m \in \Z  \\ 
        \dfrac{1-q^N}{1-q} = 0 & \text{otherwise} 
      \end{cases}
      \notag \\
  & = \delta_{kl}^\per \,,
  \label{e.d-orth}
\end{align}
where
\begin{equation}
  q = \e^{\i 2 \pi \tfrac{l-k}N}_{\vphantom 1} 
\end{equation}
so that $q^N = \e^{2\pi\i (l-k)}=1$, and where $\delta_{kl}^\per$
denotes the \emph{periodic Kronecker symbol}
\begin{equation}
  \delta_{kl}^\per =
  \begin{cases}
    1 & \text{if } l=k +m N \text{ for some } m \in \Z \,, \\ 
    0 & \text{otherwise} \,.
  \end{cases}
\end{equation}

\begin{remark}
Along the same lines, it is easy to show that for any $m \in \Z$,
\begin{equation}
  \frac1N \sum_{j=m+1}^{m+N} \e^{-\i k jh} \, \e^{\i l jh} 
  = \delta_{kl}^\per \,.
  \label{e.d-orth2}
\end{equation}
\end{remark}

\begin{remark}
In the language of linear algebra, \eqref{e.d-orth} states that the
family of vectors $\{ \boldsymbol e_k = N^{-1/2} \, (\e^{\i k x_0},
\dots, \e^{\i k x_{N-1}}) \colon k = 0, \dots, N-1 \}$ form an
orthonormal basis of $\C^N$ with respect to the standard inner
product.
\end{remark}

\begin{exercise}
Prove the \emph{discrete Parseval identity}
\begin{equation}
  \sum_{k=0}^{N-1} \lvert \tilde v_k \rvert^2
  = \frac1N \sum_{j=0}^{N-1} \lvert v_j \rvert^2 \,.
\end{equation}
\end{exercise}

\begin{exercise} \label{ex.2} 
Let $\tilde v_k$ be defined by \eqref{e.dft} for any $k\in\Z$. Show that
$\tilde v_{k+N} = \tilde v_k$.  Vice versa, set $v_{j+N}=v_j$, then show
that for every $m \in \Z$,
\begin{equation}
  \tilde v_k = \frac1N \sum_{j=m+1}^{m+N} \e^{-\i k x_j} \, v_j \,.
\end{equation}
\end{exercise}

\begin{exercise} 
Let
\begin{equation}
  w_j = \frac1N \sum_{l=0}^{N-1} u_l \, v_{j-l}
\end{equation}
with the understanding that the three $N$-tuples are periodically
extended beyond their basic range of definition on which the index
varies from $0$ to $N-1$; see Exercise~\ref{ex.2}.  Show that
\begin{equation}
  \tilde w_k = \tilde u_k \, \tilde v_k \,.
\end{equation}
\end{exercise}

\begin{remark}
Exercise~\ref{ex.2} shows, in particular, that our choice of range for
$k$, the ``wavenumber'' index, is arbitrary.  Any range of $N$
consecutive wavenumbers will do.
\end{remark}

We conclude that the \emph{inverse discrete Fourier transform} is
given by
\begin{equation}
  v_j = \sum_{k=0}^{N-1} \tilde v_k \, \e^{\i k x_j} \,,
  \label{e.idft}
\end{equation}
as can be directly verified by substituting the discrete Fourier
transform \eqref{e.dft} into the right hand sum and applying the
orthogonality relation \eqref{e.d-orth}.

\begin{remark}
The normalization convention is arbitrary.  Many software
implementations of the discrete Fourier transform, including those in
\textsf{Matlab} and \textsf{Octave}, have the factor $1/N$ in the
inverse transform \eqref{e.idft} rather than in the forward transform
\eqref{e.dft}.  The convention above has the advantage that the DFT
can be seen as the Riemann sum approximation of the continuous Fourier
transform.  If both the forward and the inverse transform get a factor
of $1/\sqrt N$, the transform is unitary, which has certain
theoretical advantages, but is usually avoided in actual code.
\end{remark}

\section{Sampling and approximation}
\label{s.approx}

From the linear algebra perspective, the previous section already
gives a complete description of the discrete Fourier transform.  In
practice, however, the DFT is often used as an approximation to the
continuous Fourier transform.  Thus, we need to define the
reconstruction of a function on $[0,2\pi]$ and must then discuss the
quality of the approximation.

To begin, assume that $u\colon [0, 2\pi] \to \C$ is continuous and
periodic, and write $u_j=u(x_j)$ for $j=0, \dots, N-1$ to denote the
values of $u$ sampled on the grid.  For simplicity, we assume that $N$
is even.  We now define a function $v$, the \emph{trigonometric
interpolant} of $u$ on the grid, by letting $\tilde v_k$ be the DFT of
the samples $u_j$,
\begin{equation}
  \tilde u_k = \frac1N \sum_{j=0}^{N-1} \e^{-\i k x_j} \, u_j 
  \label{e.hatvk}
\end{equation}
for $k=-N/2, \dots, N/2-1$, then set
\begin{equation}
  v(x) = \sum_{k=-N/2}^{N/2-1} \tilde u_k \, \e^{\i k x} \,.
  \label{e.reconstruction}
\end{equation}
Note that $k$ still ranges over a set of $N$ consecutive wave numbers.
The particular choice here seems odd, but will turn out to be crucial
for keeping the approximation error small.

The name ``trigonometric interpolant'' is justified by the fact that
$v$ and $u$ coincide on the grid, to be verified in the following
straightforward exercise.

\begin{exercise}
Show that $v(x_j)=u_j$ for $j=0, \dots, N-1$ and that $\hat v_k =
\tilde u_k$ for $k = -N/2, \dots, N/2-1$.
\end{exercise}

However, the discrete Fourier coefficients $\hat v_k = \tilde u_k$ are
generally different from the Fourier coefficients $\hat u_k$ defined
by \eqref{e.ft}!  To quantify how they differ, we compute, using
\eqref{e.hatvk}, \eqref{e.ift}, and \eqref{e.d-orth}, that
\begin{equation}
  \tilde u_k 
  = \frac1N \sum_{j=0}^{N-1} \e^{-\i kx_j} \, 
    \sum_{l \in \Z} \hat u_l \, \e^{\i l x_j}
  = \sum_{l \in \Z} \hat u_l \, \delta_{kl}^\per \,.
\end{equation}
Moving the term where $k=l$ onto the left side of the equality and
re-indexing the sum, we obtain the \emph{Poisson summation formula} or
\emph{aliasing formula}
\begin{equation}
  \tilde u_k - \hat u_k
  = \sum_{\substack{m \in \Z \\ m \neq 0}} \hat u_{k+mN} \,.
\end{equation}

This result shows that the error in the Fourier coefficients with
wavenumbers within the range $k=-N/2, \dots, N/2-1$ depends on the
magnitude of $\hat u_l$ outside of this range.  Since smooth functions
have Fourier coefficients whose magnitude decreases with increasing
$\lvert k \rvert$ (a more concrete version of this statement is given
below), and real-valued functions have $\lvert \hat u_k \rvert =
\lvert \hat u_{-k} \rvert$, the maximally symmetric choice of
wavenumber range where $k=-N/2, \dots, N/2-1$ can be expected to
minimize the approximation error for ``typical'' functions.

\begin{remark}
The aliasing formula implies a version of the \emph{Shannon sampling
theorem}, often stated as ``Exact reconstruction of a continuous-time
baseband signal from its samples is possible if the signal is
bandlimited and the sampling frequency is greater than twice the
signal bandwidth'' \cite{Shannon-wiki}.
\end{remark}

It is often useful to quantify the approximation error $u-v$ directly.
It is easiest to do under the assumption that
\begin{equation}
  \lVert u \rVert_{H^m}^2
  = \sum_{k \in \Z} (1 + \lvert k \rvert^2)^m \, 
    \lvert \hat u_k \rvert^2 < \infty \,.
  \label{e.hm}
\end{equation}

\begin{remark}
The condition above is a statement about the smoothness of $u$.  More
specifically, it says that derivatives of $u$, while they may not
exist everywhere, do not diverge too rapidly.  For $m=1$, for example,
\eqref{e.hm} is equivalent to saying that
\begin{equation}
  \int_0^{2 \pi} \bigl(
                  \lvert u(x) \rvert^2 + \lvert u'(x) \rvert^2
                \bigl) \, \d x < \infty \,.
\end{equation}
\end{remark}

\begin{remark}
The set of $L^2$ functions, i.e.\ square integrable functions, for
which \eqref{e.hm} holds, is called the Sobolev space $H^m$.
\end{remark}

We compute
\begin{align}
  \lVert u-v \rVert^2_{L^2}
  & = \int_0^{2\pi} \lvert u-v \rvert^2 \, \d x 
    = 2 \pi \sum_{k \in \Z} \lvert \hat{u}_k - \hat{v}_k \rvert^2 
      \notag \\
  & = 2 \pi \sum_{k = -N/2}^{N/2-1} 
        \lvert \hat{u}_k - \hat{v}_k \rvert^2 
      + 2 \pi \left( \sum_{k < -N/2} + \sum_{k \geq N/2} \right)
        \lvert \hat u_k \rvert^2 
      \notag \\
  & = 2 \pi \sum_{k = -N/2}^{N/2-1} 
        \left\lvert 
          \sum_{j \neq 0} \hat u_{k+jN} 
        \right\rvert^2 
      + 2 \pi \left( \sum_{k < -N/2} + \sum_{k \geq N/2} \right)
        \lvert \hat u_k \rvert^2 
      \notag \\
  & \leq 2 \pi 
      \left(
         \sum_{\lvert k \rvert \geq N/2} \lvert \hat u_k \rvert
      \right)^2
      + 2 \pi \sum_{\lvert k \rvert \geq N/2} 
        \lvert \hat u_k \rvert^2 \,,
  \label{e.error}
\end{align}
where the last inequality is based on $\sum a_i^2 \leq \bigl( \sum a_i
\bigr)^2$ for nonnegative $a_i$, and where we have symmetrized the
range of the summation index taking the more conservative bound.  The
second term on the right of \eqref{e.error} is the \emph{truncation
error}, which can be estimated as
\begin{align}
  \sum_{\lvert k \rvert \geq N/2} \lvert \hat u_k \rvert^2
  & = \sum_{\lvert k \rvert \geq N/2} 
      \lvert k \rvert^{-2m} \,  \lvert k \rvert^{2m} \,
      \lvert \hat u_k \rvert^2 
      \notag \\
  & \leq \biggl( \frac{N}2 \biggr)^{-2m}
      \sum_{\lvert k \rvert \geq N/2} \lvert k \rvert^{2m} \,
      \lvert \hat u_k \rvert^2
      \notag \\
  & = c_1(m) \, N^{-2m} \, \lVert u \rVert_{H^m}^2 \,.
  \label{e.error1}
\end{align}
The first term in \eqref{e.error} is the \emph{aliasing error}, which
is estimated as
\begin{align}
   \left(
     \sum_{\lvert k \rvert \geq N/2} \lvert \hat u_k \rvert
   \right)^2
   & = \left(
         \sum_{\lvert k \rvert \geq N/2} 
         \lvert k \rvert^{-m} \, \lvert k \rvert^m \, 
         \lvert \hat u_k \rvert
       \right)^2
       \notag \\
   & \leq \sum_{\lvert k \rvert \geq N/2} \lvert k \rvert^{-2m}
        \sum_{\lvert k \rvert \geq N/2} 
        \lvert k \rvert^{2m} \, \lvert \hat u_k \rvert^2
        \notag \\
   & \leq c_2(m) \, N^{-2m+1} \, \lVert u \rVert_{H^m}^2 \,,
  \label{e.error2}
\end{align}
where the first inequality is due to the Cauchy--Schwarz inequality,
and the second inequality is a result of estimating the sum in terms
of an integral,
\begin{equation}
  \sum_{\lvert k \rvert \geq \kappa} \lvert k \rvert^{-2m}
  \leq 2 \int_{\kappa-1}^\infty k^{-2m} \, \d k
  = \frac2{2m-1} \, (\kappa-1)^{1-2m} \,.
\end{equation}
Inserting \eqref{e.error1} and \eqref{e.error2} into \eqref{e.error},
we have proved the following.

\begin{theorem}
For every $m>\tfrac12$ there exists a constant $c(m)$ such that for
every $u \in H^m$ 
\begin{equation}
  \lVert u - v \rVert_{L^2} 
  \leq c(m) \, N^{1/2-m} \, \lVert u \rVert_{H^m} \,,
\end{equation}
where $v$ is the trigonometric interpolant of $u$ as given by
\eqref{e.reconstruction}.
\end{theorem}

\section{The fast Fourier transform}

The fast Fourier transform is mathematically equivalent to the
discrete Fourier transform, but organizes the operations in a
tree-like structure which reduces to the computational complexity from
the $O(N^2)$ when viewed as a generic vector-matrix multiplication to
$O(N\ln N)$ provided that $N$ is completely factorizable into small
prime factors.

The construction can be understood as follows.  Assume that $N$ can be
written as a product of integers $N=pq$.  Now re-index the discrete
Fourier transform, writing
\begin{gather}
  j = q \, j_1 + j_2
  \quad \text{with } j_1=0,\dots,p-1 \text{ and } j_2=0,\dots,q-1 \,,
  \\ 
  k = p \, k_2 + k_1
  \quad \text{with } k_1=0,\dots,p-1 \text{ and } k_2=0,\dots,q-1 \,.
\end{gather}
The discrete Fourier transform \eqref{e.dft} then reads
\begin{align}
  \tilde v_k 
  & = \frac1N \sum_{j_2=0}^{q-1} \sum_{j_1=0}^{p-1}
        \e^{-\i (q j_1 + j_2)(p k_2 + k_1)h} \, v_j
      \notag \\
  & = \frac1N \sum_{j_2=0}^{q-1} \e^{-\i j_2 (p k_2 + k_1)h}
      \sum_{j_1=0}^{p-1} \e^{-\i j_1 k_2 pqh} \, 
        \e^{-\i k_1 j_1 qh} \, v_j 
      \notag \\
  & = \frac1N \sum_{j_2=0}^{q-1} \e^{-\i k_2 j_2 2 \pi/q} \, 
        \e^{-\i k_1 j_2 h} 
      \sum_{j_1=0}^{p-1} \e^{-\i k_1 j_1 2 \pi/p} \, v_j \,.
  \label{e.factorization}
\end{align}
The crucial observation is that the inner sum is now itself, for $j_2$
fixed, a discrete Fourier transform of length $p$.  Thus, if $T_p$
denotes the number of terms needed to compute a length $p$ Fourier
transform, the computation of the inner sum for all values of $j_2$
requires $q \, T_p$ terms. 

Similarly, the outer sum is, for $k_1$ fixed, a discrete Fourier
transform of length $q$ acting on the values of the inner sum
pre-multiplied by $\exp(-\i k_1 j_2 h)$.  Thus, its computation
involves $p \, T_q$ terms.  

We conclude that, in total,
\begin{equation}
  T_N = q \, T_p + p \, T_q
  \label{e.recurrence}
\end{equation}
which, in general, is less than $N^2=p^2 q^2$.  Moreover, the
factorization of the DFT is recursive so long as $p$ or $q$ can
themselves be factored.

In the important special case where $N=2^k$, it is easy to solve this
recurrence.  Write $t_k = T_N$, $q=2$ and $p=2^{k-1}$, whence
\eqref{e.recurrence} takes the form
\begin{equation}
  t_k = 2 \, t_{k-1} + 2^{k-1} \, t_1 \,.
\end{equation}
Noting that self-consistency of this expression requires that $t_0=0$
(computing a length-one discrete Fourier transform is trivial and does
not require any floating point operations), we can easily find the
solution to this linear difference equation \cite{MO},
\begin{equation}
  t_k = k \, 2^{k-1} \, t_1 \,.
\end{equation}
Thus,
\begin{equation}
  T_N = \log_2 N \, 2^{\log_2 N - 1} \, t_1
      = N \log_2 N \, t_1/2
      = O(N \ln N) \,.
\end{equation}
More generally, if $N=p_1 \dots p_m$, then the discrete Fourier
transform can be computed in $N(p_1 + \dots p_m)$ operations
\cite{CT}. 



\begin{thebibliography}{99}

\bibitem{CT}
  J.W. Cooley and J.W. Tukey, 
  \textit{An algorithm for the machine calculation of complex Fourier
  series}, 
  Math. Comput. \textbf{19} (1965), 297--301.\\
  {\scriptsize\url{http://www.jstor.org/view/00255718/di970514/97p0011i/0}}

\bibitem{Folland}
  G.B. Folland, 
  ``Real Analysis.'' Wiley, 1999.

\bibitem{MaslenRockmore}
  D.K. Maslen and D.N. Rockmore,
  \textit{The Cooley-Tukey FFT and group theory},
  Notices of the AMS \textbf{48} (2001), 1151--1160. \\
  {\scriptsize\url{http://www.ams.org/notices/200110/fea-maslen.pdf}}

\bibitem{MO}
  M. Oliver, 
  \textit{Difference Equations}, 
  General Mathematics and Computational Science I Lecture Notes. \\
  {\scriptsize\url{http://math.jacobs-university.de/oliver/teaching/iub/fall2006/cps101/handouts/difference-equations.pdf}}
  
\bibitem{Trefethen}
  L.N. Trefethen,
  ``Spectral Methods in Matlab.''
  SIAM, 2001.

\bibitem{Shannon-wiki}
  ``Nyquist--Shannon sampling theorem.'' 
  Wikipedia, The Free Encyclopedia. June 13, 2007, 05:39 UTC. \\
  {\scriptsize\url{http://en.wikipedia.org/w/index.php?title=Nyquist%E2%80%93Shannon_sampling_theorem&oldid=137837501}}

\end{thebibliography}



\end{document}
