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

\title{Power Spectral Estimation With FFT \\ (Numerical Recipes Section 13.4)}
\author{C.\ Kankelborg}
\date{Rev. \today}


\begin{document}
\maketitle




\section{The Periodogram and Windowing}
\label{sec:powerspec}

Several methods have been developed for the estimation of power spectra from data (see \textit{Numerical Recipes}, \S\S\,13.4,13.7-8). The simplest FFT estimate of the power spectrum is called the \emph{periodogram} (\textit{Numerical Recipes} eq. 13.4.5). The following example uses the periodogram technique, and illustrates the need for \emph{windowing}. The accompanying figures (\ref{fig:powerspec1}, \ref{fig:powerspec2}, \ref{fig:powerspec3}) correspond to the three plots generated by the example.

\footnotesize
\verbatiminput{programs/powerspec.m}
\normalsize

\begin{figure}[h]
	\centering
	\resizebox{0.95\columnwidth}{!}{\includegraphics{programs/powerspec1.pdf}}
	\caption{The signal and its hanning-windowed counterpart.}
	\label{fig:powerspec1}
\end{figure}

\begin{figure}[h]
	\centering
	\resizebox{0.95\columnwidth}{!}{\includegraphics{programs/powerspec2.pdf}}
	\caption{The power spectrum of the signal and its hanning-windowed counterpart.}
	\label{fig:powerspec2}
\end{figure}

\begin{figure}[h]
	\centering
	\resizebox{0.95\columnwidth}{!}{\includegraphics{programs/powerspec3.pdf}}
	\caption{Same as the previous figure, but zoomed in. The periodogram calculated from 
	   the raw data has tails that fall off like $1/k$ due to the discontinuity across 
	   the domain boundary. Vertical lines compare estimates
	   of the signal frequency. The Hanning window offers a marked improvement in frequency
	   estimation and eliminates the $1/k$ tails.}
	\label{fig:powerspec3}
\end{figure}

The example calculates the power spectrum using the FFT with and without \emph{windowing} (in optics and image processing, the customary term is \emph{apodization}). The Octave/MATLAB function \texttt{hanning(N)} produces an $N$-element array of the form $\frac{1}{2}\left[ 1-\cos \left(\frac{2\pi (n-1)}{N-1}\right) \right]$. As you can see in figure \ref{fig:powerspec1}, this eliminates the wraparound discontinuity of the data. Since the Hanning window is a smooth function with a broad peak, its only side effect is that the FFT of the windowed data will be smoothed (convolved with the Fourier transform of the window, which is itself a narrow, peaked function). Without the windowing, there is a $1/\nu^2$ tail artifact in the power spectrum; this is eliminated in the power spectrum of the windowed data, at the expense of broadening the peak (figures \ref{fig:powerspec2} \& \ref{fig:powerspec3}).


\section{Repeatability and Uncertainty}

Often, we are looking for the power spectrum of some stochastic phenomenon, like the avalanche effect in a diode or the pattern of waves on the surface of the ocean. In cases like this, the periodogram is merely an estimate of the true power spectrum. Sampled data misses what came before the sampled interval, what came after, and what came in between the samples. Ideally, we would measure the signal over an infinite interval of time, at infinitesimal sampling period. But it turns out that once we have a high enough sampling frequency and a long enough duration to catch the frequencies of interest, getting more data only adds noise to the periodogram. \textit{Numerical Recipes} asserts that the standard deviation of the periodogram at any frequency is equal to its mean (100\% error!). This does not change as N increases (the periodogram reaches for higher spectral resolution rather than tightening up the error bars). One solution to this problem is to partition the data into many segments. The optimal information per data point is obtained when the second half of each segment overlaps the first half of the next segment. 

\noindent \textbf{Question:} What assumptions underlie the 100\% error estimate given by \textit{Numerical Recipes}?


\section{Example}

Begin by reviewing the following codes.

\subsection{Power spectral estimation code}
\footnotesize
\verbatiminput{programs/spectrum.m}
\normalsize

\subsection{Example time series}
\footnotesize
\verbatiminput{programs/timeseries4.m}
\normalsize

\subsection{Script to run the example and produce figures}
\footnotesize
\verbatiminput{programs/timeseries4test.m}
\normalsize


\subsection{In-class exercise}
\begin{enumerate}
   \item Run timeseries4test.m, and examine the plots.
   \item Do you see the ``100\% error'' in the ordinary periodogram (red)?
   \item Does averaging many partitions help:
   \begin{enumerate}
      \item Estimation of the red noise spectrum?
      \item Determination of the line shape for the peak broadened by phase noise?
      \item Frequency estimation for the 0.36788 Hz peak?
   \end{enumerate}
   \item Discuss the implications for power spectral analysis of \emph{resolved} and \emph{unresolved} spectral features.
\end{enumerate}

It generally helps to increase spectral resolution when looking at unresolved peaks. But with well-resolved features like the red noise and the broadened peak in the above example, better power spectral estimation is obtained by averaging many windowed partitions.

\end{document}