%% ================================================================
\section{Fourier transforms and Fourier series}

Faris \cite{bib:faris} points out that the general context of Fourier analysis
is an abelian group and its dual group.  (The duality is in the Pontryagin
sense, the details of which are beyond the scope of this paper.  See
\cite{bib:wik} for a lucid discussion.) The elements of the group are typically
thought of as time or space variables; the elements of the dual group are
thought of as angular-frequency or wave-number variables.  The following may be
shown to be dual pairs of abelian groups:
\begin{itemize}
\item[(I)] The real line and itself.  In this case, we have Fourier transforms:
	given $f \in L^2(\R, dx)$, we compute $\hat{f} \in L^2(\R, dk/2\pi)$.
\item[(II)] The unit circle and the integers $\Z$.  In this case, we have
	Fourier series: given $f \in L^2(T, dx)$, we compute $\vecc \in \ell^2(\Z)$.
\item[(III)] The integers $\Z$ and the unit circle.  This is the case for using
	Fourier series to reconstruct the original function.
\item[(IV)] A finite cyclic group of order $N$ and itself.  Here, we have
	discrete Fourier transforms.
\end{itemize}
This paper elaborates on transforms of type IV.  One may define the discrete
Fourier transform by fiat.  Instead, though, we consider two motivating
constructions.

(Note:  This paper is similar to my \emph{Poisson summation and the discrete
Fourier transform}, of which the present work is a rough draft.  The
Poisson-summation paper omits information about the fast Fourier transform
which is included here, so I have retained this draft.)

%% ================================================================
\section{The discrete Fourier transform via Riemann sums}

First, consider an approximation to a transform of type II.  Define $T_L$ to be
the circle of circumference $L$.  Let $f(x)$ be square-integrable on $T_L$,
i.e.  let $f(x)$ be a particular representative for $L^2(T_L,dx)$.
Equivalently, think of $f(x)$ as being periodic on the real line, with period
$L$.

The basis functions are of the form $e^{ikx}$, for some values of $k$.  To see
what these values must be, look at $x=L$.  For $e^{ikx}$ to have period $L$,
$kL$ must be an integer multiple of $2\pi$.  That is, we take
$$
	k = \frac{2\pi\ell}{L}
$$
for $\ell \in \Z$.

The Fourier coefficients $\hatf(k)$ are
$$
	\hatf(k) = \langle e^{ikx}, f(x) \rangle
	= \frac{1}{L} \int_0^L e^{-ikx} f(x) dx.
$$
Form a Riemann sum with $N$ evenly spaced mesh points, deferring for a moment
the well-definedness of the sum on $L^2(T_L,dx)$.  We have
$$x_j = j \Delta x$$
where
$$\Delta x = \frac{L}{N}.$$
Then the integral is approximated by
$$
	\frac{1}{L} \int_0^L e^{-ikx} f(x) dx
	\approx
	\frac{1}{L} \sum_{j=0}^{N-1} e^{-ikx_j} f(x_j) \Delta x.
$$
Now recall that the angular frequencies $k$ were also evenly spaced, namely,
with
$$
	k = \ell \Delta k \qand \Delta k = \frac{2\pi}{L}.
$$
Thus the $\Delta x$ and $\Delta k$ are not unrelated.  In fact,
$$
	\Delta k \Delta x = \frac{L}{N} \cdot \frac{2\pi}{L} = \frac{2\pi}{N}.
$$
Making substitutions, we now have for the estimated integral
\begin{eqnarray*}
	\hatf(k) &\approx& \frac{1}{L} \sum_{j=0}^{N-1} e^{-ikx_j} f(x_j) \Delta x \\
	&=& \frac{1}{L}
		\sum_{j=0}^{N-1} e^{-i \ell \Delta k j \Delta x} f(x_j) \Delta x \\
	&=& \frac{1}{L}
		\sum_{j=0}^{N-1} e^{-i \ell \frac{2\pi}{L} j \frac{L}{N}} f(x_j)
		\frac{L}{N}\\
	&=& \frac{1}{N}
		\sum_{j=0}^{N-1} e^{-i 2\pi \ell j / N} f(x_j).
\end{eqnarray*}

Observe that now there are $N$ input points, $f(x_0)$ through $f(x_{N-1})$.
Also observe that, when restricted to the mesh $x=j\Delta x$, the basis
functions
$$
	e^{ikx} = e^{i2\pi\ell j/N}
$$
have identical values when $\ell$ is shifted by multiples of $N$.  Thus the
periodicity of the coefficients arises naturally.  There are $N$ input points
and $N$ output points.  Phrased differently, a pure sinusoidal input of
frequency 1 is indistinguishable from a pure sinusoidal input of frequency
$N+1$ when the input data are sampled only at the $N$ mesh points.  (This fact
has significant repercussions in the implementation of signal-processing
systems.  See \cite{bib:lyons} for more information.)

% xxx put the nice figure here ...

An advantage of this derivation of the discrete Fourier transform is that it
uses elementary and easily visualizable Riemann sums.  However, it is intuitive
only:  it is not well-defined on $L^2(T_L, dx)$.  If $f$ and $g$ differ on a
set of measure zero with respect to Lebesgue measure, then they will have the
same Fourier coefficients on $L^2(T_L, dx)$.  However, if that measure-zero set
includes one of the mesh points, then the discrete Fourier transform as
presented here will produce different results.

%% ================================================================
\section{The DFT via counting measure}

Here is a second derivation of the discrete Fourier transform.  Define a
counting measure $\chi$ on $T_L$ via
$$
	\chi(f) = \frac{1}{N} \sum_{j=0}^{N-1} f(x_j)
$$
using the same mesh $x=j\Delta x$ as before.  Then for
$$
	f \in L^2(T_L, d\chi)
$$
we have
\begin{eqnarray*}
	\hatf(k) &=& \langle e^{ikx}, f(x) \rangle \\
	&=& \int_0^L e^{-ikx} f(x) d\chi \\
	&=& \frac{1}{N} \sum_{j=0}^{N-1} e^{-ikx_j} f(x_j) \\
	&=& \frac{1}{N} \sum_{j=0}^{N-1} e^{-i2\pi\ell j/N} f(x_j).
\end{eqnarray*}
Note that due to the finiteness of the sums, the spaces $L^p(T_L, d\chi)$, $p
\in [1, +\infty]$, are all the same.

This derivation of the DFT has the advantage that it is well-defined on
$L^2(T_L, d\chi)$.  The maximal set of measure zero is precisely
$$
	T_L \setminus \{0, \Delta x, \ldots, (N-1)\Delta x\}.
$$
Thus if $f$ and $g$ differ on a set of measure zero, they must agree only at
the mesh points; they may differ everywhere else.  This is the direct opposite
of the case when Lebesgue measure is used, where it is the set of mesh points,
rather than its complement, which has measure zero.  This makes sense: in
applications, we take continuous electromagnetic or acoustic phenomena, then
acquire only a finite number of data samples for analysis.

%Another advantage of the counting-measure derivation of the DFT is that the
%$\hatf(k)$'s are equal to the above sum, rather than being mere approximations.

%% ----------------------------------------------------------------
\section{The DFT as a linear transformation on $\C^N$}

For brevity of notation, let $\vecb$ and $\vecc$ be given by
$$
	b_j = f(x_j) \qand c_\ell = \hatf(k) = \hatf(2\pi\ell/L).
$$
For the DFT as derived above, we have $N$ input data points, namely $b_0, b_1,
b_2, \ldots, b_{N-1}$.  We also have $N$ output data points, namely $c_0, c_1,
c_2, \ldots, c_{N-1}$, given by
$$
	c_\ell = \frac{1}{N}\, \sum_{j=0}^{N-1} e^{-i2\pi \ell j/N} b_j
$$
Then the input data form a vector of length $N$.  We may take $\vecb \in
\R^N$, but it might as well be in $\C^N$.  Now, define the matrix $W$, the
\emph{Fourier matrix}, by
$$
	W_{\ell j} = e^{-i 2 \pi \ell j/N}.
$$
All $N$ Fourier coefficients, for $\ell = 0, 1, 2, \ldots, N-1$, are collected
together in nothing more than the matrix product
$$
	\vecc = \frac{1}{N}\, W \vecb.
$$
This makes sense:  when we multiply a matrix times a vector, the $\ell$th slot
of the output vector is the inner product of the $\ell$th row of the matrix
with the input vector.  The inner product on functions, via $\langle f, g
\rangle = \int \overline{f}g dx$, is replaced by the standard inner product on
$\C^N$.  The first vector in the inner product is $e^{i2\pi\ell j/N}$, for
fixed $\ell$; the second is $\vecb$.  This discrete Fourier transform is then
seen to be a linear transformation from $\C^N$ to $\C^N$.

For example, take $N=4$.  Then the matrix product is
$$
	\begin{bmatrix} c_0 \\ c_1 \\ c_2 \\ c_3 \\ \end{bmatrix}
	=
	\frac{1}{4}
	\left[\begin{array}{rrrr}
		1 &  1 &  1 &  1 \\
		1 & -i & -1 &  i \\
		1 & -1 &  1 & -1 \\
		1 &  i & -1 & -i \\
	\end{array}\right]
	\begin{bmatrix} b_0 \\ b_1 \\ b_2 \\ b_3 \\ \end{bmatrix}
$$
Note, for future reference, that although there are $N^2 = 16$ entries
in the matrix $W$, there are only $N=4$ distinct values in it.

%% ================================================================
\section{The inverse DFT}

Since the DFT is nothing more than multiplication by a particular matrix, its
inverse must be multiplication by the inverse matrix.  Recall that $W$ is given
by
$$
	(W)_{\ell j} = e^{-i2\pi \ell j/N}.
$$
One may guess from symmetry considerations that the analogous matrix
$\ovW$, given by
$$
	(\ovW)_{\ell j} = e^{i2\pi \ell j/N},
$$
is somehow involved.  The matrix product is
$$
	(W \ovW)_{\ell j} =
	\sum_{m=0}^{N-1} (W)_{km} (\ovW)_{m j} =
	\sum_{m=0}^{N-1} e^{-i2\pi km / N} e^{ i2\pi m j/ N} =
	\sum_{m=0}^{N-1} e^{i2\pi m(j-k) / N}.
$$
For $j\ne k$, we have a sum of $(j-k)$th roots of unity evenly spaced around
the circle, which by center-of-mass considerations sum to 0.  For $j=k$, on the
other hand, we have $1$ added to itself $N$ times.  Thus
$$
	W \ovW = N I
$$
where $I$ is the $N \times N$ identity matrix.  Therefore,
$$
	W^{-1} = \frac{1}{N} \ovW.
$$

%% ================================================================
\section{Scaling for the DFT}

There are two predominant conventions for scaling of the DFT, neither of which
matches $1/N$ scaling obtained by Riemann sums or counting measure.

Recall that for the Fourier transform on the real line, we measure functions in
the time domain using $dx$, whereas we measure functions in the frequency
domain using $dk/2\pi$.  That is, functions in the two domains are held up to
different rulers.

For the DFT, on the other hand, we have the Euclidean norm which is a canonical
way to measure vectors.  Take the usual Euclidean norm on $\C^N$, namely,
$$
	\|\vecb\| = \sum_{j=0}^{N-1} |b_j|^2.
$$

To see what happens to vector norms using the DFT, take the $N=4$ example as
above, and let $\vecb =  (1,0,0,0)^t$.  Then we have
$$
	\frac{1}{4}
	\left[\begin{array}{rrrr}
		1 &  1 &  1 &  1 \\
		1 & -i & -1 &  i \\
		1 & -1 &  1 & -1 \\
		1 &  i & -1 & -i \\
	\end{array}\right]
	\begin{bmatrix} 1 \\ 0 \\ 0 \\ 0 \\ \end{bmatrix}
	=
	\begin{bmatrix} 1/4 \\ 1/4 \\ 1/4 \\ 1/4 \\ \end{bmatrix}
$$
The Euclidean norm of the input vector is
$$
	\|\vecb\| = \sqrt{1^2 + 0^2 + 0^2 + 0^2} = 1,
$$
whereas the Euclidean norm of the output vector is
$$
	\|\vecc\| = \sqrt{1/16 + 1/16 + 1/16 + 1/16} = \sqrt{1/4} = 1/2.
$$

If we instead take the DFT to be scaled by $1/\sqrt N$, then we have
$$
	\frac{1}{2}
	\left[\begin{array}{rrrr}
		1 &  1 &  1 &  1 \\
		1 & -i & -1 &  i \\
		1 & -1 &  1 & -1 \\
		1 &  i & -1 & -i \\
	\end{array}\right]
	\begin{bmatrix} 1 \\ 0 \\ 0 \\ 0 \\ \end{bmatrix}.
	=
	\begin{bmatrix} 1/2 \\ 1/2 \\ 1/2 \\ 1/2 \\ \end{bmatrix}
$$
Now the Euclidean norm of the output vector is
$$
	\|\vecc\| = \sqrt{1/4 + 1/4 + 1/4 + 1/4} = \sqrt{1} = 1.
$$

Since we have $W \ovW = NI$, it is apparent that implementing a forward DFT
using $W/\sqrt N$ and implementing the reverse DFT using $\ovW/\sqrt N$
recovers the input.  That is, $W/\sqrt N$ and $\ovW/\sqrt N$ are inverses of
one another.  (In fact, both are unitary matrices.) This is a common convention
for DFTs in engineering practice.  The other common convention is to do a
forward DFT using $W$, and to do the reverse by multiplication by $\ovW/N$.

There is another pair of conventions:  one sometimes sees $\ovW$ used
for the forward DFT and $W$ used for the reverse.  Thus the lack of a single
standard in mathematical circles extends outside of academia.

In summary, desired properties for the DFT scaling are as follows:
\begin{itemize}
\item Forward DFT followed by reverse must recover the original.
\item The Parseval identity:  $\|\vecc\| = \|\vecb\|$.  This requires $\sqrt N$
scaling for forward and reverse transformations, as long as we insist on
holding input and output vectors up to the same Euclidean measure.
\end{itemize}

%% ================================================================
\section{Input and output folding}

It was seen above that the DFT output frequencies $k$ are taken mod $2\pi N/L$.
For example, frequency $2\pi/L$ is indistinguishable from $2\pi(N+1)/L$.
Commonly one takes the frequencies to be from $0$ to $N-1$.  However, one may
wish to have the frequences be zero-centered, i.e. from $-N/2$ to $(N/2)-1$.
Clearly one may accomplish this by performing the DFT, then interchanging the
first $N/2$ elements of the output vector with the last $N/2$ samples.  This
bulk data movement is called \emph{output folding}.

Likewise, when performing an inverse DFT, if the frequencies have been folded,
then the processing should account for that.  That is, suppose a time-domain
vector has been forward-transformed and then output-folded.  One may undo these
operations by interchanging the first and second halves of the frequency
vector, then performing the inverse DFT.  From the point of view of the inverse
DFT, then input needs to be folded.

It is perhaps less apparent that \emph{output folding} may be accomplished by
negating every other sample of the \emph{input}, and vice versa.  (In fact, for
certain implementations, these negations are more efficient than the data
movement required for bulk folding.)  This claim needs some justification.

First, see what happens when every other sample of the input is negated.
\begin{eqnarray*}
	c_\ell &=& \sum_{j=0}^{N-1} e^{-i2\pi \ell j/N} b_j\\
	\sum_{j=0}^{N-1} e^{-i2\pi \ell j/N} (-1)^j b_j
	&=& \sum_{j=0}^{N-1} e^{-i2\pi \ell j/N} e^{-i\pi j} b_j \\
	&=& \sum_{j=0}^{N-1} e^{-i2\pi \ell j/N} e^{-i2\pi j(N/2)/N} b_j \\
	&=& \sum_{j=0}^{N-1} e^{-i2\pi (\ell j + j(N/2))/N} b_j \\
	&=& \sum_{j=0}^{N-1} e^{-i2\pi (\ell + (N/2))j/N} b_j \\
	&=& c_{\ell+N/2}
\end{eqnarray*}
where the $\ell$ index is taken mod $N$.

Second, suppose every other sample of the output is negated.
\begin{eqnarray*}
	(-1)^\ell c_\ell &=& (-1)^\ell \sum_{j=0}^{N-1} e^{-i2\pi \ell j/N} b_j\\
	&=& e^{-i\pi\ell}\sum_{j=0}^{N-1} e^{-i2\pi \ell j/N} b_j\\
	&=& e^{-i\pi\ell(N/2)/N} \sum_{j=0}^{N-1} e^{-i2\pi \ell j/N} b_j\\
	&=& \sum_{j=0}^{N-1} e^{-i2\pi \ell (j+N/2)/N} b_j\\
\end{eqnarray*}
where the $j$ index is taken mod $N$.

Thus the folding claims are proved.

%% ================================================================
\section{The fast Fourier transform}

This section is purely expository.  Here I motivate, for the pure-math reader,
the derivation of the fast Fourier transform.  There is little mathematical
content in this section, and in fact that is the point:  the FFT is nothing
more than a clever trick for obtaining precisely the same end result as the
DFT, using less computation time.  When doing theoretical analysis, one uses
the DFT rather than the FFT, since the DFT has an elegant expression as a
matrix product.

Note that multiplying an $N \times N$ matrix by a vector of length $N$ requires
on the order of $N^2$ arithmetic operations.  To see this, observe that the
$j$th element of the output vector is nothing more than the dot product of the
$j$th row of the matrix with the input vector.  There are $N$ of these dot
products, and each of them contains $N$ products and $N$ sums.  Thus one often
hears that the DFT is $O(N^2)$.

The fast Fourier transform is a trick that was originally discovered by Gauss,
and has been rediscovered at least half a dozen times over the intervening
centuries.  The rediscovery which stuck in the minds of the technical community
was that of Cooley and Tukey in 1965, at which point computing hardware was
available to make practical use of such trickery.  More historical information
may be found in \cite{bib:nr}.  The trick is illustrated by example.

Again take $N=4$.  Writing out the matrix product in full detail, we have
$$
	\begin{array}{rrrrrrrrr}
		c_0 &=&   b_0 &+&   b_1 &+&   b_2 &+&   b_3 \\
		c_1 &=&   b_0 &-& i b_1 &-&   b_2 &+& i b_3 \\
		c_2 &=&   b_0 &-&   b_1 &+&   b_2 &-&   b_3 \\
		c_3 &=&   b_0 &+& i b_1 &-&   b_2 &-& i b_3 \\
	\end{array}
$$
Not all of the 16 sums and 16 products appear to the eye; for this small
value of $N$, often we are multiplying by 1.

Note that the sum $b_0+b_2$ is computed twice, as is the difference $b_0-b_2$.
Re-ordering the sums, we see
$$
	\begin{array}{rrrrrrrrr}
		c_0 &=&   (b_0 &+&   b_2) &+&   (b_1 &+&   b_3) \\
		c_1 &=&   (b_0 &-&   b_2) &-& i (b_1 &-&   b_3) \\
		c_2 &=&   (b_0 &+&   b_2) &-&   (b_1 &+&   b_3) \\
		c_3 &=&   (b_0 &-&   b_2) &+& i (b_1 &-&   b_3) \\
	\end{array}
$$
Here we have the four sums $(b_0 \pm b_2)$ and $(b_1 \pm b_3)$.  These are
followed by four products, namely multiplying by $\pm 1$ and $\pm i$, and four
more sums to finish up.  Thus there are now 8 sums and 4 products, down from 16
of each.  This in the essence of the fast Fourier transform.

To generalize the concept beyond $N=4$, define the \emph{two-point butterfly
operation} (this is a standard term) as follows.  This is formally a map from
$\C^3 \to \C^2$, given by
$$
	(a, b, w) \mapsto (a+wb, a-wb).
$$
One writes

\setlength{\unitlength}{1.0in}
\begin{picture}(6,0.6)
%\put(0,0){\framebox(6,0.6)}

\put(2.00, 0.40){$a$}
\put(2.00, 0.15){$b$}
\put(2.25, 0.15){$w$}

\put(4.00, 0.40){$a+wb$}
\put(4.00, 0.15){$a-wb$}

\put(2.95, 0.32){\line(-2, 1){.25}}
\put(2.95, 0.28){\line(-2,-1){.25}}
\put(3.05, 0.32){\line( 2, 1){.25}}
\put(3.05, 0.28){\line( 2,-1){.25}}
\put(3.00, 0.30){\circle{.10}}
\end{picture}

Using this notation, one may diagram the 4-point FFT as follows:

\setlength{\unitlength}{1.0in}
\begin{picture}(6,1.2)
%\put(0,0){\framebox(6,1.2)}

\put(0.10, 0.90){$b_0$}
\put(0.10, 0.65){$b_2$}
\put(0.35, 0.65){$1$}
%
\put(0.10, 0.40){$b_1$}
\put(0.10, 0.15){$b_3$}
\put(0.35, 0.15){$1$}
%
\put(0.95, 0.82){\line(-2, 1){.25}}
\put(0.95, 0.78){\line(-2,-1){.25}}
\put(1.05, 0.82){\line( 2, 1){.25}}
\put(1.05, 0.78){\line( 2,-1){.25}}
\put(1.00, 0.80){\circle{.10}}
%
\put(0.95, 0.32){\line(-2, 1){.25}}
\put(0.95, 0.28){\line(-2,-1){.25}}
\put(1.05, 0.32){\line( 2, 1){.25}}
\put(1.05, 0.28){\line( 2,-1){.25}}
\put(1.00, 0.30){\circle{.10}}
%
\put(1.60, 0.90){$b_0+b_2$}
\put(1.60, 0.65){$b_0-b_2$}
\put(2.30, 0.40){$1$}
%
\put(1.60, 0.40){$b_1+b_3$}
\put(1.60, 0.15){$b_1-b_3$}
\put(2.25, 0.15){$-i$}


\put(2.75, 0.72){\line(-1, 1){.25}}
\put(2.75, 0.68){\line(-1,-1){.25}}
\put(2.85, 0.72){\line( 1, 1){.25}}
\put(2.85, 0.68){\line( 1,-1){.25}}
\put(2.80, 0.70){\circle{.10}}
%
\put(2.75, 0.42){\line(-1, 1){.25}}
\put(2.75, 0.38){\line(-1,-1){.25}}
\put(2.85, 0.42){\line( 1, 1){.25}}
\put(2.85, 0.38){\line( 1,-1){.25}}
\put(2.80, 0.40){\circle{.10}}

\put(3.50, 0.90){$(b_0 + b_2) +  (b_1 + b_3)$}
\put(3.50, 0.65){$(b_0 - b_2) - i(b_1 - b_3)$}
%
\put(3.50, 0.40){$(b_0 + b_2) -  (b_1 + b_3)$}
\put(3.50, 0.15){$(b_0 - b_2) + i(b_1 - b_3)$}

\end{picture}

The general implementation of a fast Fourier transform for $N$ a power of 2
follows this model.  The following computational steps are required:

\begin{itemize}

\item Re-ordering the input.  For $N=4$, we changed $b_0, b_1, b_2, b_3$ to
$b_0, b_2, b_1, b_3$.  In general, one writes the input in \emph{bit-reversed
order}.  (See \cite{bib:nr} for full details.  A discussion of bits and bytes
is tempting, but outside the scope of this paper.)  There are no sums or
products here, only data movement and/or fancy indexing of arrays.

\item Each column of butterflies is called a \emph{stage}.

\item Each stage requires on the order of $N$ sums and products.

\item There are $\log_2(N)$ stages.

\item At each stage, one accesses $N$th roots of unity.  One does not need a
$W$ matrix with $N^2$ elements; there are only $N$ distinct complex numbers in
the matrix.  Furthermore, it may be shown that only the first $N/2$ of them are
actually used during the butterfly operations.  These values may be precomputed
and stored in a table; the technical term is \emph{twiddle table}.

\end{itemize}

As a result of the above facts, the FFT, for $N$ a power of 2, has
computational complexity $O(N \log_2 N)$.  This is a vast improvement on the
$O(N^2)$ of the DFT.  For example, take $N=1024$.  Then $N^2 = 1,048,576$
whereas $N \log_2 N = 10240$, i.e. a 1024-point FFT may be executed in a
hundredth the time required for a 1024-point DFT.

For $N$ a power of 2, scaling by $1/\sqrt N$ may be accomplished by multiplying
intermediate resultss by $1/\sqrt 2$ at each stage, say, after the butterfly
operations.  Likewise, scaling by $1/N$ may be done by scaling by $1/2$ at each
stage.

Last, suppose $N$ is not a power of two.  One may perform a derivation as given
above, but using, say, 3-point butterflies when $3$ divides $N$.  In fact, some
simplification may be done using \emph{any} non-trivial factorization of $N$.
When $N$ is prime, only a DFT is possible; the more small factors $N$ has, the
better.  Some software packages will factor $N$ and do the best-possible DFT;
other packages permit $N$ to only be a power of 2, or perhaps 3 times a power
of two, etc.

Full information about implementations of FFTs may be found in
\cite{bib:nr}.  See also \cite{bib:kerl} for a particular power-of-two
implementation.

%% ================================================================
\section{The FFT via factoring the Fourier matrix}

Implementations of the FFT use bit-reverse, twiddle tables, and butterfly
operations as described in the previous section; an implementor needs nothing
more.  Yet one may wish to provide some mathematical justification, beyond a
worked example for a specific value of $N$.

To motivate the theory, examine the 4-point FFT from the previous section.
Writing the bit-reverse, stage-1 and stage-2 operations as matrices, we have
$$
	\left[\begin{array}{r}
		c_0 \\ c_1 \\ c_2 \\ c_3 \\
	\end{array}\right]
	=
	\left[\begin{array}{rrrr}
		1 & 0 & 1 & 0 \\
		0 & 1 & 0 &-i \\
		1 & 0 &-1 & 0 \\
		0 & 1 & 0 & i \\
	\end{array}\right]
	\left[\begin{array}{rrrr}
		1 & 1 & 0 & 0 \\
		1 &-1 & 0 & 0 \\
		0 & 0 & 1 & 1 \\
		0 & 0 & 1 &-1 \\
	\end{array}\right]
	\left[\begin{array}{rrrr}
		1 & 0 & 0 & 0 \\
		0 & 0 & 1 & 0 \\
		0 & 1 & 0 & 0 \\
		0 & 0 & 0 & 1 \\
	\end{array}\right]
	\left[\begin{array}{r}
		b_0 \\ b_1 \\ b_2 \\ b_3 \\
	\end{array}\right]
$$
That is, the Fourier matrix has been factored into sparser components.

Now let $N$ be an arbitrary even positive integer.  Separating even- and
odd-numbered inputs, we may write the DFT as
\begin{eqnarray*}
	c_\ell &=& \sum_{j=0}^{N-1} e^{-i2\pi\ell j/N} b_j \\
	&=& \sum_{j=0}^{N/2-1} e^{-i2\pi\ell 2j    /N} b_{2j}
	   +\sum_{j=0}^{N/2-1} e^{-i2\pi\ell (2j+1)/N} b_{2j+1} \\
	&=& \sum_{j=0}^{N/2-1} e^{-i2\pi\ell 2j/N} b_{2j}
	   + e^{-i2\pi \ell/N}
	    \sum_{j=0}^{N/2-1} e^{-i2\pi\ell 2j/N} b_{2j+1} \\
\end{eqnarray*}
Observe that we now have, for each $\ell$, a pair of sums, each corresponding
to a DFT of half the size as the original.  This process may be applied
recursively for each power of 2 dividing $N$; the recursion ends with a DFT as
described in previous sections.  Optimally, when $N$ is a power of two, the
final DFT is a 1-point DFT.  The construction just outlined is known as the
\emph{Danielson-Lanczos lemma}.

%% ================================================================
\section{Higher-dimensional DFTs}

Let $T^d_\vecL$ be a $d$-dimensional torus, where $\vecL = (L_1, \ldots, L_d)$
for positive real $L_i$'s.  Write $\vecx$ and $\veck$ for $(x_1, \ldots, x_d)$
and $(k_1, \ldots, k_d)$, respectively.  To find the possible values of
$\veck$, we again require
$$
	k_r = \frac{2\pi\ell_r}{L_r}
$$
for $\ell_r \in \Z$.
Then for $f(\vecx) \in L^2(T^N_\vecL,d\vecx)$ we have the $d$-dimensional
type-II transform (i.e. Fourier series)
\begin{eqnarray*}
	\hatf(\veck)
	= \int_{T^d_\vecL} e^{-i \veck \cdot \vecx} f(\vecx) d\vecx
\end{eqnarray*}
If in addition $f(\vecx) \in L^1(T^N_\vecL,d\vecx)$, by Fubini's theorem we may
iterate the integrals to obtain
\begin{eqnarray*}
	\hatf(k_1, \ldots, k_d)
	&=& \int_0^{L_1} \cdots \int_0^{L_d} e^{-i(k_1 x_1 + \ldots + k_d x_d)}
		f(x_1, \ldots, x_d) dx_d \cdots dx_1 \\
	&=& \int_0^{L_1} e^{-ik_1 x_1}
		\left[ \cdots \left[
		\int_0^{L_d} e^{-ik_d x_d}
		f(x_1, \ldots, x_d) dx_d
		\right] \cdots \right]
		dx_1.
\end{eqnarray*}

Now form a discrete mesh
$$
	\vecx_\vecj = (j_1\Delta x_1, \ldots, j_d\Delta x_d)
$$
with, for $r=1, 2, \ldots, d$,
$$
	L_r = N_r \Delta x_r.
$$
Define a counting measure $\chi$ via
$$
	\chi(f(\vecx)) = \frac{1}{\prod_{r=1}^d N_d}
		\sum_{j_1=0}^{N_1-1} \cdots \sum_{j_d=0}^{N_d-1}
		f(x_{j_1}, \ldots, x_{j_d}).
$$
As in the 1-dimensional case, due to the finiteness of the sums the spaces
$L^p(T^N_\vecL,d\chi)$, $p \in [1, +\infty]$, are all the same.

One may again form a Riemann sum, or take the integral using the counting
measure.  One obtains
\begin{eqnarray*}
	\hatf(\veck)
	&=& \int_{T^d_\vecL} e^{-i \veck \cdot \vecx} f(\vecx) d\chi \\
	&=& \sum_{j_1=0}^{N_1}
		\frac{e^{-i2\pi \ell_1 j_1/N_1}}{N_1}
		\cdots
		\sum_{j_d=0}^{N_d}
		\frac{e^{-i2\pi \ell_d j_d/N_d}}{N_d}
		f(x_{j_1}, \ldots, x_{j_d})
\end{eqnarray*}

%% ================================================================
\section{Decomposition of higher-dimensional DFTs}

Take the DFT as in the previous section to be a linear transformation on
$\C^{N_1} \times\cdots\times \C^{N_d}$.
For $r = 1, 2, \ldots, d$ and $j_r = 0, 1, 2, \ldots, N_r-1$ write
$$
	b_{j_1, \ldots, j_d} = f(x_{j_1}, \ldots, x_{j_d}) \qqand
	c_{\ell_1, \ldots, \ell_d} = \hatf(k_{\ell_1}, \ldots, k_{\ell_d}).
$$
Then
\begin{eqnarray*}
	c_{\ell_1, \ldots, \ell_d}
	&=& \sum_{j_1=0}^{N_1}
		\frac{e^{-i2\pi \ell_1 j_1/N_1}}{N_1}
		\cdots
		\sum_{j_d=0}^{N_d}
		\frac{e^{-i2\pi \ell_d j_d/N_d}}{N_d}
		b_{j_1, \ldots, j_d}.
\end{eqnarray*}
As before, each of the $j_r$'s and $\ell_r$'s are taken mod $N_r$.  There are
$\prod_{r=1}^d N_r$ output values $c_{\ell_1, \ldots, \ell_d}$, and for each of
them there are on the order of $\prod_{r=1}^d N_r$ arithmetic operations to
evaluate the nested sum.  Thus, it would appear that the $d$-dimensional DFT
has complexity $O(N_1^2 \cdots N_d^2)$.  In the particular case when all the
$N_r$'s are the same, say $N_r=N$, this is $O(N^{2d})$.

However, one may decompose the $d$-dimensional DFT as a sequence of
1-dimensional DFTs, avoiding some redundant computation.  The key is that the
inner sum is invariant under the outer sums, and thus does not need to be
recomputed, and so on for each sum.  Let
$$
	\omega_m = e^{-i2\pi/m}.
$$
Form the sequence $I^{(0)}, \ldots, I^{(d)}$ of $d$-dimensional arrays via
\begin{eqnarray*}
	I^{(0)}_{j_1, \ldots, j_d} &=& b_{j_1, \ldots, j_d}\\
	I^{(1)}_{\ell_1, j_2, \ldots, j_d} &=&
		\sum_{j_1=0}^{N_1}\frac{\omega_{N_1}^{\ell_1, j_1}}{N_1}
			I^{(0)}_{j_1, \ldots, j_d}\\
	I^{(2)}_{\ell_1, \ell_2, j_3, \ldots, j_d} &=&
		\sum_{j_2=0}^{N_2}\frac{\omega_{N_2}^{\ell_2, j_2}}{N_2}
			I^{(1)}_{\ell_1, j_2, \ldots, j_d}\\
	&\vdots& \\
	I^{(d)}_{\ell_1, \ldots, \ell_d}
		&=& \sum_{j_d=0}^{N_d}\frac{\omega_{N_d}^{\ell_d, j_d}}{N_d}
			I^{d-1}_{\ell_1, \ldots, \ell_{d-1}, j_d}\\
		&=& b_{\ell_1, \ldots, \ell_d}\\
\end{eqnarray*}
Each sum visibly consists of $O(N_d)$ operations, and there are $\prod_{r=1}^d
N_r$ elements in each of the $I^{(r)}$ arrays.  Thus the total
complexity is
$$
	O(N_1^2 N_2 \cdots N_d + N_1 N_2^2 \cdots N_d + \ldots + N_1
		\cdots N_{d-1} N_d^2).
$$
In the particular case when all $N_r=N$, this is $O(dN^{d+1})$.

Another special case is when $d=2$.  Here, the input data form a $N_1 \times
N_2$ matrix, and the decomposition reduces to doing 1-dimensional DFTs along
each row of the matrix, then doing 1-dimensional DFTs along each column.  (For
this reason, decomposition of the multi-dimensional DFT into 1-dimensional DFTs
is often referred to as the \emph{row-column algorithm}.) That is, there are
$N_1$ 1-dimensional DFTs with complexity $O(N_2^2)$, followed by $N_2$
1-dimensional DFTs with complexity $O(N_1^2)$.

The 1-dimensional DFTs may be replaced with FFTs, with optimal performance when
each of the $N_r$'s is a power of two.  In this case, the complexity becomes
\begin{eqnarray*}
	&&O(N_1 \log_2(N_1) N_2 \cdots N_d + N_1 N_2 \log_2(N_2) \cdots N_d + \ldots
		+ N_1 \cdots N_{d-1} N_d \log_2(N_d)) \\
	&=& O\left( [N_1 \cdots N_d)\right] [\log_2(N_1) +\ldots+ \log_2(N_d)]
	= O\left(\left[\prod_{r=1}^d N_d\right]
		\left[\sum_{r=1}^d \log_2(N_d)\right]\right) \\
	&=& O\left(\left[\prod_{r=1}^d N_d\right]
		\left[\log_2\left(\prod_{r=1}^d N_d \right)\right]\right)
	= O(M \log_2(M))
\end{eqnarray*}
if we define $M = \prod_{r=1}^d N_r$.  Thus, the performance benefit of the
one-dimensional fast Fourier transform extends to higher dimensions.

% d  bad good
% 1  N^2  N^2 =  1  (N^2)
% 2  N^4 2N^3 = 2N  (N^2)
% 3  N^6 3N^4 = 3N^2(N^2)
% 4  N^8 4N^5 = 4N^3(N^2)

%As one would imagine, FFTs of various dimensions arise in practice: for
%example, there are 1-dimensional problems in audio processing, 2-dimensional
%problems in image processing, and 3-dimensional problems in geophysics.

% xxx curse of dimensionality still applies.
