% -*- mode: Noweb; noweb-code-mode: c-mode -*-

\bibliographystyle{apalike}

\subsection{The Fast Discrete Fourier Transform}

\begin{fullwidth}

The motivating problem is spectral analysis: 
Think of sampling a continuous time signal $a(t)$ at $N$ uniformly spaced times.
The discrete Fourier transform (DFT) map these samples to values in the
frequency domain, linear combinations of sine and cosine functions.

The discrete Fourier transform (DFT) transforms a sequence $\seq{a}=\sequence{a_{0},\,a_{1},\ldots,\,a_{N-1}}$
of $N=2^{n}$ complex numbers into another sequence $\seq{y}=\sequence{y_{0},\,y_{0},\ldots,\,y_{N-1}}$.
The computation for $y_{k}$ is:
\begin{equation}
  \label{fft:eq:1}
  y_{k} = \sum_{0\leq j < N}a_{j}\cdot \left(e^{-2\pi i k/N}\right)^{j}\tag{$k=0,\ldots,\,N-1$}  
\end{equation}
The discrete Fourier transform can be represented as a matrix equation
where $\omega=e^{-2\pi i/N}$ (Use $\omega_{N}$ when syntactic sugar is needed).
\[\seq{y}=W\seq{a}=
\begin{bmatrix}
 1 & 1          & 1          & 1          & \cdots & 1\\
 1 & \omega     & \omega^{2} & \omega^{3} & \cdots & \omega^{N-1}\\
 1 & \omega^{2} & \omega^{4} & \omega^{6} & \cdots & \omega ^{2(N-1)}\\
 1 & \omega^{3} & \omega^{6} & \omega^{9} & \cdots &\omega^{3(N-1)}\\
\vdots & \vdots & \vdots     & \vdots     &\ddots  & \vdots \\
 1 & \omega^{N-1} &\omega^{2(N-1)} &\omega^{3(N-1)} & \cdots & \omega^{(N-1)(N-1)}
\end{bmatrix}
\begin{bmatrix}
  a_{0} \\
  a_{1} \\
  a_{2} \\
  a_{3} \\
\vdots \\
  a_{N-1} \\
\end{bmatrix}\]
Computing $\seq{y}$ using matrix--vector multiplication 
requires not only evaluation the expressions $\omega^{k}$ for $k=0,\ldots (N-1)^{2}$, but also $N(N-1)$ additions and
$N^{2}$ multiplications.  


To see if we can do better, let's start with the two point Fourier transform, where $\omega = e^{-\pi i}=-1$.
\[ \begin{bmatrix}[rr]
   1 &  1 \\
   1 & -1 
 \end{bmatrix}\tag{The Fourier transform of $[a_{0},\,a_{1}]$ is their sum and difference.}
\]
Next, consider the four point Fourier transform matrix, where $\omega=e^{-\pi i/2}=-i$. 
Since $\omega^{2}=-\omega^{0}=-1$ any power $\omega^{m}$ can be reduced to $\pm\omega^{r}$ where
$r=m\bmod{2}$ and the sign depends on the parity mod $4$ of $m$.
Lastly, permute the columns putting the even columns first. See how this partitions the matrix into four quadrants.
\begin{align*}
 W & = \begin{bmatrix}
  \omega^{0} & \omega^{0} & \omega^{0} & \omega^{0} \\
  \omega^{0} & \omega^{1} & \omega^{2} & \omega^{3} \\
  \omega^{0} & \omega^{2} & \omega^{4} & \omega^{6} \\
  \omega^{0} & \omega^{3} & \omega^{6} & \omega^{9} 
\end{bmatrix} =
\begin{bmatrix}[rrrr]
  \omega^{0} & \omega^{0} &  \omega^{0} & \omega^{0} \\
  \omega^{0} & \omega^{1} & -\omega^{0} & -\omega^{1} \\
  \omega^{0} & -\omega^{0} &  \omega^{0} & -\omega^{0} \\
  \omega^{0} & -\omega^{1} &  -\omega^{0} & \omega^{1} 
\end{bmatrix} 
 = 
 \begin{bmatrix}[rrrr]
   1 &  1 &  1 &  1 \\
   1 & -i & -1 &  i \\
   1 & -1 &  1 & -1 \\
   1 &  i & -i & -i 
 \end{bmatrix}\\
 & \mapsto
\begin{bmatrix}[rr|rr]
  \omega^{0} &  \omega^{0} &  \omega^{0}  & \omega^{0} \\
  \omega^{0} & -\omega^{0} &  \omega^{1}  & -\omega^{1} \\\midrule
  \omega^{0} &  \omega^{0} &  -\omega^{0} & -\omega^{0} \\
  \omega^{0} & -\omega^{0} &  -\omega^{1} & \omega^{1} 
\end{bmatrix} =
 \begin{bmatrix}[rr|rr]
   1 &  1 &  1 &  1 \\
   1 & -1 & -i &  i \\\midrule
   1 &  1 & -1 & -1 \\
   1 & -1 &  i & -i 
 \end{bmatrix}
\end{align*}
Apply this last matrix operator to even--odd samples $[a_{0},\,a_{2},\,a_{1},\,a_{3}]$.
This computes, in steps, the DFT of $[a_{0},\,a_{2}]$:
\[y'_{0}=a_{0}+a_{2},\,  y'_{2}=a_{0}-a_{2}\] 
and 
\[y'_{1}=a_{1}+a_{3},\,  y'_{3}=a_{1}-a_{3}\] 
the DFT of $[a_{1},\, a_{3}]$.
Then, combine the results.
\[ [y_{0},\,y_{1},\,y_{2},\,y_{3}] = [y'_{0}+y'_{1},\, y'_{2}+\omega y'_{3},\, y'_{0}-y'_{1},\, y'_{2}-\omega y'_{3}] \]
This is the Fourier transform, check that it is so.
\begin{align*}
y_{0} & = y'_{0}+y'_{1}        =  (a_{0}+a_{2})+(a_{1}+a_{3}) \\
y_{1} & = y'_{2}+\omega y'_{3} =  (a_{0}-a_{2})-i(a_{1}-a_{3}) \\
y_{2} & = y'_{0}-y'_{1}        =  (a_{0}+a_{2})-(a_{1}+a_{3}) \\  
y_{3} & = y'_{2}-\omega y'_{3} =  (a_{0}-a_{2})+i(a_{1}-a_{3})
\end{align*}
The eight point case further reveals the pattern.  Now, $\omega=e^{-2\pi i/8}$ and $\omega^{4}=-1$.
Therefore, any power $\omega^{m}$ can be reduced to $\pm\omega^{r}$ where 
$r=m\bmod{4}$ and the sign depends on the parity mod $8$ of $m$.
Lastly, permute the columns putting the even columns first. See how this partitions the matrix into four quadrants.
\begin{align*}
W &  = 
  %\begin{bmatrix}
  %\omega^{0} & \omega^{0} & \omega^{0} & \omega^{0} & \omega^{0} & \omega^{0} & \omega^{0} & \omega^{0}\\ 
  %\omega^{0} & \omega^{1} & \omega^{2} & \omega^{3} & \omega^{4} & \omega^{5} & \omega^{6} & \omega^{7}\\
  %\omega^{0} & \omega^{2} & \omega^{4} & \omega^{6} & \omega^{0} & \omega^{2} & \omega^{4} & \omega^{6}\\
  %\omega^{0} & \omega^{3} & \omega^{6} & \omega^{1} & \omega^{4} & \omega^{7} & \omega^{2} & \omega^{5}\\
  %\omega^{0} & \omega^{4} & \omega^{0} & \omega^{4} & \omega^{0} & \omega^{4} & \omega^{0} & \omega^{4}\\
  %\omega^{0} & \omega^{5} & \omega^{2} & \omega^{7} & \omega^{4} & \omega^{1} & \omega^{6} & \omega^{3}\\
  %\omega^{0} & \omega^{6} & \omega^{4} & \omega^{2} & \omega^{0} & \omega^{6} & \omega^{4} & \omega^{2}\\
  %\omega^{0} & \omega^{7} & \omega^{6} & \omega^{5} & \omega^{4} & \omega^{3} & \omega^{2} & \omega^{1}\\
  %\end{bmatrix}
  % =
  \begin{bmatrix}[rrrrrrrr]
  \omega^{0} & \omega^{0}  & \omega^{0}  & \omega^{0} & \omega^{0}  & \omega^{0} & \omega^{0} & \omega^{0}\\ 
  \omega^{0} & \omega^{1}  & \omega^{2}  & \omega^{3} & -\omega^{0} & -\omega^{1} & -\omega^{2} & -\omega^{3}\\
  \omega^{0} & \omega^{2}  & -\omega^{0} & -\omega^{2} & \omega^{0}  & \omega^{2} & -\omega^{0} & -\omega^{2}\\
  \omega^{0} & \omega^{3}  & -\omega^{2} & \omega^{1} & -\omega^{0}  & -\omega^{3} & \omega^{2} & -\omega^{1}\\
  \omega^{0} & -\omega^{0} & \omega^{0} & -\omega^{0} & \omega^{0}  & -\omega^{0} & \omega^{0} & -\omega^{0}\\
  \omega^{0} & -\omega^{1} & \omega^{2}  & -\omega^{3} & -\omega^{0}  & \omega^{1} & -\omega^{2} & \omega^{3}\\
  \omega^{0} & -\omega^{2} & -\omega^{0}  & \omega^{2} & \omega^{0}  & -\omega^{2} & -\omega^{0} & \omega^{2}\\
  \omega^{0} & -\omega^{3} & -\omega^{2}  & -\omega^{1} & -\omega^{0}  & \omega^{3} & \omega^{2} & \omega^{1}\\
  \end{bmatrix}\\
 &   \mapsto
  \begin{bmatrix}[rrrr|rrrr]
  \omega^{0} & \omega^{0}  & \omega^{0}  & \omega^{0}  & \omega^{0}  &  \omega^{0}  & \omega^{0}  & \omega^{0}\\ 
  \omega^{0} & \omega^{2}  & -\omega^{0} & -\omega^{2} & \omega^{1}  &  \omega^{3}  & -\omega^{1} & -\omega^{3}\\
  \omega^{0} & -\omega^{0} & \omega^{0}  & -\omega^{0} & \omega^{2}  &  -\omega^{2} & \omega^{2}  & -\omega^{2}\\
  \omega^{0} & -\omega^{2} & -\omega^{0} & \omega^{2}  & \omega^{3}  &  \omega^{1}  & -\omega^{3} & -\omega^{1}\\\midrule
  \omega^{0} & \omega^{0}  & \omega^{0}  & \omega^{0}  & -\omega^{0} &  -\omega^{0} & -\omega^{0} & -\omega^{0}\\
  \omega^{0} & \omega^{2}  & -\omega^{0} & -\omega^{2} & -\omega^{1} &  -\omega^{3} & \omega^{1}  & \omega^{3}\\
  \omega^{0} & -\omega^{0} & \omega^{0}  & -\omega^{0} & -\omega^{2} &  \omega^{2}  & -\omega^{2} & \omega^{2}\\
  \omega^{0} & -\omega^{2} & -\omega^{0} & \omega^{2}  & -\omega^{3} &  -\omega^{1} & \omega^{3}  & \omega^{1}\\
  \end{bmatrix}
  = 
  \begin{bmatrix}[rrrr|rrrr]
   1 & 1  &  1  &  1  &  1         & 1        & 1        &1        \\
   1 & -i &  -1 &  i  &  \omega    & -i\omega & -\omega  &i\omega  \\
   1 & -1 &  1  &  -1 &  -i        & i        & -i       &i        \\
   1 & i  &  -1 &  -i &  -i\omega  & \omega   & i\omega  &-\omega  \\\midrule
   1 & 1  &  1  &  1  &  -1        & -1       & -1       &-1       \\
   1 & -i &  -1 &  i  &  -\omega   & i\omega  & \omega   &-i\omega \\
   1 & -1 &  1  &  -1 &  i         & -i       & i        &-i       \\
   1 & i  &  -1 &  -i &  i\omega   & -\omega  & -i\omega &\omega 
 \end{bmatrix}\\
%& = 
%  \begin{bmatrix}
%   1 & 1         & 1  & 1        & 1  & 1        & 1  & 1        \\
%   1 & \omega    & -i & -i\omega & -1 & -\omega  & i  & i\omega  \\
%   1 & -i        & -1 & i        & 1  & -i       & -1 & i        \\
%   1 & -i\omega  & i  & \omega   & -1 & i\omega  & -i & -\omega  \\
%   1 & -1        & 1  & -1       & 1  & -1       & 1  & -1       \\
%   1 & -\omega   & -i & i\omega  & -1 & \omega   & i  & -i\omega \\
%   1 & i         & -1 & -i       & 1  & i        & -1 & -i       \\
%   1 & i\omega   & i  & -\omega  & -1 & -i\omega & -i & \omega 
% \end{bmatrix}\\
\end{align*}
The Fourier transform of $8$ samples can be computed by constructing
the Fourier transform of $[a_{0},\,a_{2},\,a_{4},\,a_{6}]$
and $[a_{1},\,a_{3},\,a_{5},\,a_{7}]$
and then combining the results
\[ 
  [y_{0},\,y_{2},\,y_{4},\,y_{6}]+[y_{1},\,\omega y_{3},\,\omega^{2} y_{5},\,\omega^{3} y_{7}],
  [y_{0},\,y_{2},\,y_{4},\,y_{6}]-[y_{1},\,\omega y_{3},\,\omega^{2} y_{5},\,\omega^{3} y_{7}]
\]
The upshot is: To compute the Fourier transform of $N=2^{n}$ samples, 
compute the Fourier transform of $N/2$ samples twice.
Then combine the results 
\[\seq{y}_{\text{even}}+\seq{\omega}\oplus\seq{y}_{\text{odd}},\, \seq{y}_{\text{even}}-\seq{\omega}\oplus\seq{y}_{\text{odd}}\]
where 
\[\seq{\omega}\oplus\seq{y}_{\text{odd}}=[y_{1}+\omega y_{3}+\omega^{2} y_{5}+\cdots+\omega^{\floor{N/2}} y_{N-1}]\]
Cooley and Tukey~\citep{cotu:65} 
are credited with discovering how to make the
discrete Fourier transform fast by exploiting facts about polynomials and
\emph{roots of unity}.
Their fast discrete Fourier transform (FDFT) computes the values
$y_{k},\,k=0,\ldots,\,N-1$ in $O(N\lg{N})$ steps.
A huge improvement. 

To better understand the FDFT some basic facts are needed.
Consider the problem in terms of polynomials.
Let 
\[y(x)=a_{0} + a_{1}x + a_{2}x^{2} + a_{3}x^{3} + \cdots + a_{N-1}x^{N-1}\] 
be a polynomial of degree $N-1$. The polynomial can be split by the parity (even or odd) of its coefficients. 
 \[y_{\text{even}}(x)=a_{0} + a_{2}x + a_{4}x^{2} + \cdots + a_{N-2}x^{N/2-1}\]
  and 
  \[y_{\text{odd}}(x)=a_{1} + a_{3}x + a_{5}x^{2} + \cdots + a_{N-1}x^{N/2-1}\]
 and $y(x)$ can be reconstructed by
 \begin{equation}
   \label{fft:eq:2}
   y(x) = y_{\text{even}}(x^{2}) + xy_{\text{odd}}(x^{2})
 \end{equation}
  The DFT calls for evaluating $y(x)$ at $x=\omega^{k}=e^{-2\pi i k/N}$
  for $k=0,\ldots,\,N-1$. 
  This reduces to evaluating $y_{\text{even}}(x^{2})$ and $y_{\text{odd}}(x^{2})$
  at $x^{2}=\omega^{k}=e^{-4\pi i k/N}$ for $0\leq k < N/2$.

  The values $\omega^{k}$ are $N$-th complex roots of unity: They lie on the unit circle in the complex plain.
  \[\left(\omega^{k}\right)^{N} = e^{-2\pi ik} = \cos{{2\pi k}}-i\sin{2\pi k}=1\tag{$0\leq k <N-1$}\]
  The \emph{principal $N$-th root of unity} is
    \[\omega=\omega_{N} = e^{-2\pi i/N}\]
Further,
\begin{align*}
  \omega^{N/2}                    & = \left(e^{{-2\pi i/N}}\right)^{N/2}=e^{-\pi i}=-1 \\
  \omega^{k+N/2}                  & = \left(e^{-2\pi i(k+N/b2}\right)=e^{-2\pi ik}e^{-\pi iN}=-\omega^{k}\\
  \left(\omega^{k+N/2}\right)^{2} & = \left(e^{-2\pi i(2k+N}\right)=e^{-2\pi i(2k)}=(\omega^{k})^{2}
\end{align*}

\end{fullwidth}

\subsection{Developing the Cooley--Tukey FFT in Haskell}

Here's Haskell code~\thanks{I am indebted to a former student,
  Cristhian Cruz Espinoza, for the Haskell code structure.} for the
fast discrete Fourier transform. 
The outline for the {\tt{}fft} module is simple: Import other needed modules, define other functions as needed, and include
a {\tt{}main} function for IO and testing.

\begin{fullwidth}
\nwfilename{FFT.nw}\nwbegincode{1}\sublabel{NWFFT6-FFT3-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWFFT6-FFT3-1}}}\moddef{FFT~{\nwtagstyle{}\subpageref{NWFFT6-FFT3-1}}}\endmoddef
module FFT ( fft ) where
\LA{}Import external modules~{\nwtagstyle{}\subpageref{NWFFT6-ImpN-1}}\RA{}
\LA{}Auxiliary functions~{\nwtagstyle{}\subpageref{NWFFT6-AuxJ-1}}\RA{}
\LA{}Main function~{\nwtagstyle{}\subpageref{NWFFT6-MaiD-1}}\RA{}
\nwnotused{FFT}\nwendcode{}\nwbegindocs{2}\nwdocspar
The {\tt{}Complex} data type is needed for constructing complex numbers and performing arithmetic on values of this type. 
\nwenddocs{}\nwbegincode{3}\sublabel{NWFFT6-ImpN-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWFFT6-ImpN-1}}}\moddef{Import external modules~{\nwtagstyle{}\subpageref{NWFFT6-ImpN-1}}}\endmoddef
import Data.Complex 
\nwalsodefined{\\{NWFFT6-ImpN-2}}\nwused{\\{NWFFT6-FFT3-1}}\nwendcode{}\nwbegindocs{4}\nwdocspar
Looking forward, here are descriptions of the functions that will be needed.
The value {\tt{}n} is what has been called $N$ previously, a quirk of Haskell here.
\nwenddocs{}\nwbegincode{5}\sublabel{NWFFT6-AuxJ-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWFFT6-AuxJ-1}}}\moddef{Auxiliary functions~{\nwtagstyle{}\subpageref{NWFFT6-AuxJ-1}}}\endmoddef
\LA{}Compute a root of unity $e^{(-2\pi k/n)j}$~{\nwtagstyle{}\subpageref{NWFFT6-n)jg-1}}\RA{}
\LA{}Split the list into even and odd lists~{\nwtagstyle{}\subpageref{NWFFT6-Splc-1}}\RA{}
\LA{}Combine the lists~{\nwtagstyle{}\subpageref{NWFFT6-ComH-1}}\RA{}
\LA{}Cooley-Tukey FDFT~{\nwtagstyle{}\subpageref{NWFFT6-CooH-1}}\RA{} 
\nwalsodefined{\\{NWFFT6-AuxJ-2}}\nwused{\\{NWFFT6-FFT3-1}}\nwendcode{}\nwbegindocs{6}\nwdocspar
The {\tt{}mkPolar} function maps an amplitude $r$ and phase $\theta$ to complex number $re^{\theta i}=r(\cos\theta + i\sin\theta)$.
For roots of unity, the amplitude is $r=1$, and the phase is $\theta=-2\pi i jk/n$.
\nwenddocs{}\nwbegincode{7}\sublabel{NWFFT6-n)jg-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWFFT6-n)jg-1}}}\moddef{Compute a root of unity $e^{(-2\pi k/n)j}$~{\nwtagstyle{}\subpageref{NWFFT6-n)jg-1}}}\endmoddef
cmpExp n k j = mkPolar 1 (-2 * pi * fromIntegral k / fromIntegral n) * j 
\nwused{\\{NWFFT6-AuxJ-1}}\nwendcode{}\nwbegindocs{8}\nwdocspar
To apply the recursion, the list of samples must be split into two pairs: The first list contains 
sample with even indexes; the second contains odd indexed samples.
\nwenddocs{}\nwbegincode{9}\sublabel{NWFFT6-Splc-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWFFT6-Splc-1}}}\moddef{Split the list into even and odd lists~{\nwtagstyle{}\subpageref{NWFFT6-Splc-1}}}\endmoddef
splitEvenOdd [] = ([], []) 
splitEvenOdd [x] = ([x],[]) 
splitEvenOdd (x:y:xs) = (x:xe, y:yo) 
    where (xe, yo) = splitEvenOdd xs
\nwused{\\{NWFFT6-AuxJ-1}}\nwendcode{}\nwbegindocs{10}\nwdocspar
Combine the two lists into one.
The recursion requires concatenating the sum and difference of the two lists:

\nwenddocs{}\nwbegincode{11}\sublabel{NWFFT6-ComH-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWFFT6-ComH-1}}}\moddef{Combine the lists~{\nwtagstyle{}\subpageref{NWFFT6-ComH-1}}}\endmoddef
combine [] [] = [] 
combine xs [] = xs 
combine [] ys = ys 
combine xs ys = add ++ sub 
    where add = zipWith (+) xs ys 
          sub = zipWith (-) xs ys
\nwused{\\{NWFFT6-AuxJ-1}}\nwendcode{}\nwbegindocs{12}\nwdocspar

\nwenddocs{}\nwbegincode{13}\sublabel{NWFFT6-CooH-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWFFT6-CooH-1}}}\moddef{Cooley-Tukey FDFT~{\nwtagstyle{}\subpageref{NWFFT6-CooH-1}}}\endmoddef
fft :: [Complex Double] -> [Complex Double] 
fft [] = [] 
fft [x] = [x] 
fft xs = combine even cmpOdd 
    where (e, o) = splitEvenOdd xs 
          even = fft e 
          odd = fft o 
          size = length xs 
          cmpOdd = zipWith (cmpExp size) [0..] odd 
\nwused{\\{NWFFT6-AuxJ-1}}\nwendcode{}\nwbegindocs{14}\nwdocspar

\subsection{Testing the code}

Okay, we've got some code, let's test that it works as expected.
Interactive tests are good for unit testing functions.
%Of course, I've been testing it throughout its development.

\subsection{Unit Tests with {\tt{}ghci}}

The Haskell source \hlink{(should be) here}{http://cs.fit.edu/~wds/classes/aa/Code/FFT.hs}.
The small sample script below will get you started. The script shows the results of tests on computing roots of unity, 
splitting and combining lists, and the computing the Fourier transform.  
\end{fullwidth}
\todo[noline]{Verify the results in the scripts and test other computations yourself.}
Recall the values {\tt{}cmpExp\ n\ k\ j}$=(\omega^{k})^{j}=e^{-2\pi kj/n}$ are used in evaluating
values $y_{k}$ in the Fourier transform.
\[y_{k} = \sum_{0\leq j < n} a_{j}\omega^{kj}\]
\begin{fullwidth}
Therefore,
{\tt{}cmpExp\ 2\ 0\ 1}$=e^{-2\pi 0\cdot1/2}=1$ is computed correctly,
but there is round-off error in the imaginary part of {\tt{}cmpExp\ 2\ 1\ 1}.

\begin{verbatim}
% ghci
GHCi, version (your version \#): http://www.haskell.org/ghc/  :? for help
Prelude> :l FFT.hs 
[1 of 1] Compiling Main             ( FFT.hs, interpreted )
Ok, modules loaded: Main.

*FFT> cmpExp 2 0 1
1.0 :+ 0.0
*FFT> cmpExp 2 1 1
(-1.0) :+ (-1.2246467991473532e-16)
*FFT> splitEvenOdd [0,1,2,3]
([0,2],[1,3])
combine [1,2] [3,4]
[4,6,-2,-2]
*FFT> fft [1,1]
[2.0 :+ 0.0,0.0 :+ 0.0]
\end{verbatim}

\subsection{Integration Tests with {\tt{}ghc}}

\end{fullwidth}
\todo[noline]{Analyze the time complexities of {\tt{}splitEvenOdd}, {\tt{}combine}, and {\tt{}fft}.}
The topic of the class is analysis of algorithms.  
Theory says the time complexity is modeled by some function in $\Theta(N\lg{N})=\Theta(2^{n}n)$ where $N=2^{n}$ is the
number of samples. Does running the {\tt{}fft} code produce running times that fit this model?

\hypertarget{subsec:testfft}{The \mynoweb source file \texttt{FFT.nw} is
\hlink{here}{http://cs.fit.edu/~wds/classes/aa/Code/FFT.nw}.}
Running {\tt{}notangle\ FFT.nw\ >\ FFT.hs} generates the Haskell source
\hlink{FFT.hs}{http://cs.fit.edu/~wds/classes/aa/Code/FFT.hs}
Running{\tt{}noweave\ -index\ -delay\ FFT.nw\ >\ FFT.tex} generates 
the \LaTeX~source \hlink{FFT.tex}{http://cs.fit.edu/~wds/classes/aa/Code/FFT.tex}
\todo[noline]{Compile the code with profiling options, collect run time data and compare your experiment with theory.}
The overall idea of the {\tt{}main} function is simple: Read an integer $n$; create a list of
$N=2^{n}$ (random) complex numbers; compute their Fourier transform; and show the results.

\begin{fullwidth}
\nwenddocs{}\nwbegincode{15}\sublabel{NWFFT6-MaiD-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWFFT6-MaiD-1}}}\moddef{Main function~{\nwtagstyle{}\subpageref{NWFFT6-MaiD-1}}}\endmoddef
main :: IO() 
main = interact (show . fft . createList . read) 
\nwused{\\{NWFFT6-FFT3-1}}\nwendcode{}\nwbegindocs{16}\nwdocspar

A list of $N=2^{n}$ random complex numbers will be generated.
To do this, the {\tt{}System.Random} module will be used.
\nwenddocs{}\nwbegincode{17}\sublabel{NWFFT6-ImpN-2}\nwmargintag{{\nwtagstyle{}\subpageref{NWFFT6-ImpN-2}}}\moddef{Import external modules~{\nwtagstyle{}\subpageref{NWFFT6-ImpN-1}}}\plusendmoddef
import System.Random 
\nwendcode{}\nwbegindocs{18}\nwdocspar

To generate a list of complex numbers, create lists of 
\emph{real} and \emph{imaginary} parts, and zip them together with 
the {\tt{}Complex} number constructor {\tt{}:+}.
\nwenddocs{}\nwbegincode{19}\sublabel{NWFFT6-AuxJ-2}\nwmargintag{{\nwtagstyle{}\subpageref{NWFFT6-AuxJ-2}}}\moddef{Auxiliary functions~{\nwtagstyle{}\subpageref{NWFFT6-AuxJ-1}}}\plusendmoddef
\LA{}Create a list of random complex numbers~{\nwtagstyle{}\subpageref{NWFFT6-Cred-1}}\RA{}

\nwendcode{}\nwbegincode{20}\sublabel{NWFFT6-Cred-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWFFT6-Cred-1}}}\moddef{Create a list of random complex numbers~{\nwtagstyle{}\subpageref{NWFFT6-Cred-1}}}\endmoddef
createList :: Int -> [Complex Double] 
createList n = take (2^n) $ zipWith (:+) (randoms (mkStdGen 10)) (randoms (mkStdGen 10)) 
\nwused{\\{NWFFT6-AuxJ-2}}\nwendcode{}

\nwixlogsorted{c}{{Auxiliary functions}{NWFFT6-AuxJ-1}{\nwixu{NWFFT6-FFT3-1}\nwixd{NWFFT6-AuxJ-1}\nwixd{NWFFT6-AuxJ-2}}}%
\nwixlogsorted{c}{{Combine the lists}{NWFFT6-ComH-1}{\nwixu{NWFFT6-AuxJ-1}\nwixd{NWFFT6-ComH-1}}}%
\nwixlogsorted{c}{{Compute a root of unity $e^{(-2\pi k/n)j}$}{NWFFT6-n)jg-1}{\nwixu{NWFFT6-AuxJ-1}\nwixd{NWFFT6-n)jg-1}}}%
\nwixlogsorted{c}{{Cooley-Tukey FDFT}{NWFFT6-CooH-1}{\nwixu{NWFFT6-AuxJ-1}\nwixd{NWFFT6-CooH-1}}}%
\nwixlogsorted{c}{{Create a list of random complex numbers}{NWFFT6-Cred-1}{\nwixu{NWFFT6-AuxJ-2}\nwixd{NWFFT6-Cred-1}}}%
\nwixlogsorted{c}{{FFT}{NWFFT6-FFT3-1}{\nwixd{NWFFT6-FFT3-1}}}%
\nwixlogsorted{c}{{Import external modules}{NWFFT6-ImpN-1}{\nwixu{NWFFT6-FFT3-1}\nwixd{NWFFT6-ImpN-1}\nwixd{NWFFT6-ImpN-2}}}%
\nwixlogsorted{c}{{Main function}{NWFFT6-MaiD-1}{\nwixu{NWFFT6-FFT3-1}\nwixd{NWFFT6-MaiD-1}}}%
\nwixlogsorted{c}{{Split the list into even and odd lists}{NWFFT6-Splc-1}{\nwixu{NWFFT6-AuxJ-1}\nwixd{NWFFT6-Splc-1}}}%
\nwbegindocs{21}\nwdocspar

\bibliography{cs}
  
\end{fullwidth}
\nwenddocs{}
