\documentclass[12pt]{article}
\usepackage{graphics}
\usepackage{color}
\usepackage{verbatim}
\usepackage{hyperref}

\title{DFT and FFT}
\author{C.\ Kankelborg}
\date{Rev. \today}


\begin{document}
\maketitle

\section{Introduction}

The Fourier transform is a powerful tool in the solution of linear systems, including:
\begin{itemize}
   \item Inhomogeneous ODEs (e.g.\ frequency response, impulse response)
   \item Inhomogeneous PDEs (e.g.\ scattering, diffraction, diffusion)
   \item Linear integral equations  (e.g.\ deconvolution, tomography)
\end{itemize}
All of these applications can be realized numerically on coordinate grids in 1 or more dimensions using the \emph{discrete Fourier transform} (DFT), typically implemented as a \emph{fast fourier transform} (FFT). In addition to the above applications, the FFT offers a computationally efficient approach to a wide range of signal and image processing tasks:
\begin{itemize}
   \item Power spectral estimation and peak bagging
   \item Ideal interpolation of time series and images
   \item Digital filtering
   \item Compression algorithms
\end{itemize}
This paper summarizes some important practical aspects of using the FFT for newcomers or old hands who have grown rusty. For those seeking a more comprehensive treatment of continuous or discrete Fourier transforms, I recommend Ronald Bracewell's classic \textit{The Fourier Transform and Its Applications}.


\section{Definition and Properties}

\subsection{DFT and its Inverse}
\label{sec:definition}

Given a time series $f_n$ for uniformly sampled times $t_n = (n-1)\,\Delta t$, where $1 \le n \le N$, the DFT is defined as follows:
\begin{equation} \label{eq:dft}
   \tilde{f}_k = \sum_{n=1}^{N} f_n\,e^{-2\pi i(k-1)(n-1)/N}.
\end{equation}
I have reluctantly used the awkward $n-1$ and $k-1$ instead of just $n$ and $k$ because MATLAB/Octave uses unit-offset arrays (indices $1..N$). The inverse transform is:
\begin{equation} \label{eq:idft}
   f_n = \frac{1}{N}\sum_{k=1}^{N} \tilde{f}_k\,e^{2\pi i(k-1)(n-1)/N}.
\end{equation}
The frequency index $k$ takes on $N$ discrete values. Therefore arrays $f$ and $\tilde{f}$ contain an equal number of (in general, complex) elements. \emph{In principle, there is no loss of information with the transform or its inverse}. In practice, there is roundoff error, as the following example demonstrates. The example is also a good illustration of the DFT normalization (table \ref{tab:properties}).

\footnotesize
\verbatiminput{ifftExample.txt}
\normalsize

\subsection{DFT Properties}

The DFT is useful because it possesses all the useful properties of the Fourier transform. Some examples are listed in table \ref{tab:properties}. There also exist analogues to the Fourier scaling, translation, differentiation and integration theorems. Like the inverse, all of these  hold exactly in the absence of roundoff error.

\begin{table}[h]
\caption{Some properties of the DFT.}
\label{tab:properties}
\begin{minipage}{\columnwidth}
\begin{center}
\begin{tabular}{c c}
\textbf{Normalization} 
   & $\tilde{f}_0 = \sum_n f_n $  \\
\textbf{Symmetry} & real, even $ \iff $ real, even \\
    of transform  & imaginary, even $ \iff $ imaginary, even \\
        pairs     & real, odd $ \iff $ imaginary, odd \\
% Linearity & DFT$(f + g) = \tilde{f} + \tilde{g}$ \\
% Differentiation &  \\
% Integration &  \\
% Scaling & DFT$\left[ f(at) \right] = (1/a)\,\tilde{f}(\omega / a)$ \\
% Translation & DFT$\left[ f(a+t) \right] =  e^{i\omega a} \tilde{f}$\\
%             &  \\
\textbf{Convolution} & DFT$\left[ f * g \right] = \tilde{f}\,\tilde{g}$ \\
            & DFT$\left[ fg \right] = \left(\frac{1}{N}\right)
               \tilde{f} * \tilde{g}$ \\
\textbf{Correlation} & DFT$\left[ f \otimes g \right] = 
                  \tilde{f}^{^*}\tilde{g}$ \\
               & DFT$\left[ f^{^*}g \right] = \left(\frac{1}{N}\right)
                  \tilde{f} \otimes \tilde{g}$ \\
\textbf{Parseval's theorem} & $\sum_n {f_n}^{^*}f_n = 
   \left(\frac{1}{N}\right)
   \sum_k \tilde{f_k}^{^*}\tilde{f_k}$ \\
\end{tabular}
\end{center}
\end{minipage}
\end{table}

\section{Being Discrete}

Though the DFT works much like the Fourier transform, there are some subtleties intrinsic to working with discrete data over a finite interval.


\subsection{Finite Domain}

Since the time domain is of finite extent, and the basis functions all satisfy periodic boundary conditions over the domain, the DFT is really more like a Fourier series than a Fourier transform. The time series $f_n$ is therefore treated as if it corresponds to a periodic function. Consequently, a discontinuity (in $f$ or its derivatives) from $f_N$ to $f_1$ can produce artifacts in $\tilde{f}$. For this reason, \emph{detrending} and \emph{windowing} are commonly used. See comments in \S\S\,\ref{sec:implementations}, \ref{sec:powerspec}, and \ref{sec:interpolation}.



\subsection{Nyquist Frequency}

The shortest period that can be meaningfully represented in our time series is $2\,\Delta t$. This corresponds to $k=\frac{N}{2}+1$. The corresponding highest frequency is called the Nyquist frequency,
\begin{equation} \label{eq:nyquist}
   \nu_c = \frac{1}{2\,\Delta t}.
\end{equation}


\subsection{Arrangement of Frequencies}
\label{sec:frequencies}

First-time users of the DFT may be surprised by the order in which the frequencies are stored (figure \ref{fig:fftsketch}). Where do the so-called ``negative frequencies'' come from? We will see that this is a natural consequence of the definition of the DFT.

The left half of the frequency diagram, from DC (zero frequency) to the Nyquist frequency, looks reasonable. In physical units, the frequencies are evidently
\begin{equation}\label{eq:plusnu}
   \nu = \frac{k-1}{N\,\Delta t}.
\end{equation}
The Nyquist frequency (as we saw in the last section) occurs only about halfway though the progression of frequencies. How can we physically interpret a frequency ``higher'' than the Nyquist frequency? The range of interest is
\[
   \frac{N}{2}+1 < k \le N.
\]
On the above range, let us use a new variable $k'$:
\[
   (k-1) = N - (k'-1), \quad   2 \le k' < \frac{N}{2}+1 
\]
In terms of the new variable, the DFT is
\[
   \tilde{f}_k = \sum_{n=1}^N f_n e^{-2\pi i [N-(k'-1)](n-1)/N}
      = \sum_{n=1}^N f_n e^{-2\pi i [-(k'-1)](n-1)/N}.
\]
It is now evident that $k'$ represents a negative frequency, in direct analogy to equation \ref{eq:plusnu}:
\[
   \nu = \frac{-(k'-1)}{N\,\Delta t}, \quad -\nu_c < \nu < 0.
\]
Note that the function $e^{-i\omega t}$ is orthogonal to $e^{i\omega t}$. Therefore the negative frequencies cannot be ignored in general. However, if the signal $f$ is known to be real-valued, then the information in the negative frequencies is redundant. 

\noindent Questions for study:
\begin{enumerate}
   \item For a positive frequency of index $k$, what is the index of the
      corresponding negative frequency?
   \item If $f$ is real, and the Fourier coefficient $\tilde{f}_k$ is 
      known, what is the Fourier coefficient for the 
      corresponding negative frequency?
   \item Sketch the waveforms of the sine and cosine at the Nyquist
      frequency. When these are sampled, can they both be measured?
   \item Why is there only one coefficient in the DFT at the Nyquist
      frequency? Does it correspond to a positive frequency, or a
      negative one?
\end{enumerate}

\begin{figure}[h]
	\centering
	\input{fftsketch.pdf_t}
	\caption{Arrangement of frequencies in the DFT of a time series
	  with an even number of elements. If there are an odd number of
	  elements, the Nyquist frequency is omitted.}
	\label{fig:fftsketch}
\end{figure}

Following is another simple example using Octave. In the first case, {\tt v} is a signal at the Nyquist frequency, but with a DC offset (mean value) of $\frac{1}{2}$. The transform therefore contains two Kronecker delta functions, one in the first element and one in element 5. In the second part of the example, the frequency is half the Nyquist frequency. Since the original signal is real, the values in the positive frequency bins are complex conjugates of the corresponding negative frequency bins. The same symmetry is evident in the example in \S\,\ref{sec:definition}.

\footnotesize
\verbatiminput{simpleFFTexample.txt}
\normalsize




\subsection{Aliasing}
\label{sec:aliasing}

A given sampling rate has limited bandwidth --- a definite range of frequencies that can be represented by the sampled data: $|\nu| < \nu_c$. Unfortunately, this does not mean that a signal with a frequency outside this range will remain undetected. In \S\,\ref{sec:frequencies}, we saw that there is a correspondence between negative frequencies and frequencies in the range $\nu_c < \nu < 2 \nu_c$. Figure \ref{fig:aliasing} shows explicitly how a frequency above the Nyquist cutoff will be \emph{aliased} back into the sampling passband. 

\begin{figure}[h]
	\centering
	\input{aliasing.pdf_t}
	\caption{Two signals sampled at interval $\Delta t$. The sampled data,
	  which are identical, are circled by yellow-shaded circles. A frequency of 
	  $\nu=\frac{3}{4\Delta t}$ has the same observational signature as 
	  $\nu' = \frac{-1}{4\Delta t}$.}
	\label{fig:aliasing}
\end{figure}

It turns out that every frequency outside the passband gets aliased back into the passband. The program given in Appendix \ref{sec:aliastest} is a numerical experiment used to explore this mapping (figure \ref{fig:aliastest}). If you really understand the example, then you will be able to predict how it would look if the complex exponential in {\tt aliastest.m} were replaced by a sine or a cosine. 

\begin{figure}[h]
	\centering
	\resizebox{0.95\columnwidth}{!}{\includegraphics{aliastest_figure.pdf}}
	\caption{Aliasing of any frequency into the range $|\nu | < \nu_c$,
	  as marked by the dotted lines. The input signal is $e^{2\pi i\nu t}$.
	  The source code for this example is in Appendix \ref{sec:source}.}
	\label{fig:aliastest}
\end{figure}



\subsection{Shannon's Sampling Theorem}

Any periodic function $f(t)$ (or equivalently, a function defined only on the interval $0 < t < T$) can be represented as a Fourier series,
\[
   f(t) = \sum_{k=-\infty}^{\infty} a_k e^{2\pi i(k-1) t/T}.
\]
Let us suppose that all we know of $f(t)$ is its value at $N$ evenly spaced points:
\[
   f(t_n) = f_n, \quad t_n = \frac{n-1}{N} T, \quad 1 \le n \le N.
\]
Let us further suppose that
\[
   a_k = 0 \quad \mbox{for} \quad \left| k-1 \right| > \frac{N}{2},
      \quad \mbox{and} \quad a_{1+N/2} = -a_{1-N/2}.
\]
In other words, there are no frequencies in $f$ above the Nyquist frequency. This is the definition of a \emph{band-limited} signal. Under this condition, there are only as many \emph{distinct, nonzero} Fourier coefficients as there are data points. Indeed, we note by comparison with equation \ref{eq:idft} that $a_k = f_k/N$. It follows that we can use the DFT (equation \ref{eq:dft}) to construct the exact Fourier series for $f(t)$. This demonstrates the validity of Shannon's theorem: 

\begin{quotation}
\noindent If a function $f(t)$ contains no frequencies higher than $\nu$, it is completely determined by giving its values $f_n$ at a series of points spaced $\Delta t = 1/(2\nu)$ apart.
\end{quotation}

The application of Shannon's sampling theorem to interpolation is investigated in \S\,\ref{sec:interpolation}.





\section{Fast Fourier Transform (FFT)}

The DFT as defined in equation \ref{eq:dft} uses a sum of $N$ terms to find $f_k$ for each value of $k$. The whole transform ($1 \le k \le N$) therefore requires $O(N^2)$ terms, each involving the evaluation of either a complex exponential or trig functions. This is computationally daunting for large data sets.

\subsection{The ``Fast'' Part}

Fortunately, there exists a fast Fourier transform (FFT) algorithm that can be completed in $O(N \log N)$ operations. Since the FFT is widely available in every programming language, there is little point in dwelling on the details (see \textit{Numerical Recipes}, Ch.\ 12). The key points are
\begin{enumerate}
   \item Factoring of $N$ to split $f_n$ into smaller sub-arrays.
   \item Combining the DFTs of the sub-arrays to form the DFT of
      the whole array.
   \item Clever use of multiple-angle formulae
      to minimize calls to trig functions.
\end{enumerate}
This all works most efficiently if $N$ is a power of 2. 


\subsection{FFT implementations}
\label{sec:implementations}

The only important qualities in an FFT implementation are speed, speed, and speed. Well, OK, speed, versatility, and minimizing \emph{Roundoff error}. 
Appendix \ref{sec:roundoff} briefly considers roundoff error because it can afflict practically any sort of numerical computation.


There are many implementations of the FFT. Here are several that you might encounter:
\begin{itemize}
   \item \textbf{Numerical Recipes} Ch.\ 12 contains a simple
      FFT and several variations. The \textit{NR} FFT is mainly
      a pedagogical tool, uses single precision arithmetic, and 
      is not especially versatile 
      ($N$ must be a power of 2). Nevertheless, it would presumably
      work well enough for many applications.
   \item \textbf{FFTW} (the ``Fastest Fourier Transform in the West,''
      see \url{http://www.fftw.org}) was developed, somewhat ironically,
      at MIT.\footnote{The
      FFTW website points out the MIT is west of Italy, the home of
      the spaghetti western.} FFTW is a highly optimized, open source 
      package that is widely used.
   \item \textbf{GNU Octave} uses FFTW. Relevant functions include
      {\tt fft}, {\tt ifft}, and {\tt fftw}. Detrending and windowing
      are provided through functions such as {\tt detrend}, 
      {\tt hanning}, and {\tt hamming}.
   \item \textbf{MATLAB} has its own proprietary FFT code that is
      well-optimized.
\end{itemize}

Many FFT implementations such as FFTW and \textit{Numerical Recipes} leave out the $1/N$ in the inverse, so that it is not quite the inverse. Although the FFT requires $O(N \log N)$ operations, multiplication is more computationally expensive than addition; therefore elimintating $N$ multiplication operations results in non-negligible time savings for some high performance tasks. Since many FFT applications involve multiplying the transform by a filter function, it is most efficient to include the normalization as part of the filter function. Although Octave uses FFTW, they have included the $1/N$ factor (presumably for compatibility with MATLAB).



\section{Application: Fourier Interpolation}
\label{sec:interpolation}

Shannon's theorem implies that a band-limited signal $f(t)$ for which we have $N_1$ samples, $f_n$, can be reconstructed exactly. All we need to do is supply the higher frequency Fourier coefficients associated with measuring $f(t)$ at some higher resolution, $N_2 > N_1$. Since the signal is band-limited, this is a trivial exercise: all the new Fourier coefficients will be zero! The algorithm is simply:
\begin{enumerate}
   \item Transform the time series $f_n$ to obtain $\tilde{f}_k$.
   \item Form an expanded transform, $\tilde{f}'_k$, by snipping 
      the array $\tilde{f}_k$ at the Nyquist frequency,
      and inserting $N_2 - N_1$ zeros in the middle.
   \item Inverse transform $\tilde{f}'_k$ to form the interpolated time 
      series $f'_n$, which will now have $N_2$ elements.
\end{enumerate}

As a first example, look at figure \ref{fig:fstretch_undetrended}. Recall that the DFT treats the time series as periodic. The  large change from $f_N$ to $f_1$ 
is impossible to accommodate in a band-limited signal without ringing. What we must do is \emph{detrend} the data. The simplest way is to remove a linear trend with the right slope so that the end data ponts are brought to the same level. After calculating the interpolants, we will then add the linear trend that was removed. For this purpose, the \texttt{detrend} function in Octave/MATLAB is unfortunately useless. I have instead built a detrend algorithm into \texttt{fstretch}, my general purpose FFT interpolation routine (\S\,\ref{sec:fstretch}). A better result using \texttt{fstretch} with detrending is shown in figure \ref{fig:fstretch}. The example script \texttt{fstretchExample.m} is in Appendix \ref{sec:fstretchExample}.

\begin{figure}[h]
	\centering
	\resizebox{0.95\columnwidth}{!}{\includegraphics{fstretch_undetrended.pdf}}
	\caption{Fourier interpolation of data. The discontinuity between 
	$f_N$ and $f_1$ causes ringing.}
	\label{fig:fstretch_undetrended}
\end{figure}

\begin{figure}[h]
	\centering
	\resizebox{0.95\columnwidth}{!}{\includegraphics{fstretch_detrended.pdf}}
	\caption{Fourier interpolation of detrended data.}
	\label{fig:fstretch}
\end{figure}


\appendix

\section{Roundoff Error}
\label{sec:roundoff}

Roundoff error results from approximating real numbers with a finite number of digits. As operations are performed, the roundoff error gradually accumulates. Using our discrete Fourier transform as an example, suppose that we have $N=10^8$ samples, the signal is merely a constant, $f_n = 1$, and the internal representation of the machine is good to about 7 significant figures (this is roughly true of single precision floating point numbers; numerical codes nearly always employ double precision---except \textit{Numerical Recipes}!). When we evaluate the first (zero frequency, $k=1$) Fourier coefficient, it is simply
\[
   \tilde{f}_1 = \sum_{n=1}^{N} 1.
\]
After the first 10 million terms, we will be adding 1 to approximately 10 million. Because of the limited precision however, $10^7+1$ rounds to $10^7$, and all the remaining terms are similarly ignored. We end up adding $10^8$ ones to obtain only $10^7$! A more realistic example would show a modest loss of precision even for a smaller array at double precision. Fortunately, the FFT tends to minimize roundoff error because it subdivides the data into small segments and combines them hierarchically.


\section{Source Code}
\label{sec:source}

\subsection{\texttt{aliastest.m}}
\label{sec:aliastest}
\footnotesize
\verbatiminput{programs/aliastest.m}
\normalsize

\subsection{\texttt{fstretch.m}}
\label{sec:fstretch}
\footnotesize
\verbatiminput{programs/fstretch.m}
\normalsize

\subsection{\texttt{fstretchExample.m}}
\label{sec:fstretchExample}
\footnotesize
\verbatiminput{programs/fstretchExample.m}
\normalsize

\end{document}