%% ================================================================
\section{Fourier transforms and Fourier series}
\label{sec:intro}

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.) 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.

%% ================================================================
\section{The discrete Fourier transform via Riemann sums}
\label{sec:riemann_sums}

Consider an approximation to a transform of type II.  Define $T_L$ to be the
circle of circumference $L$.  Let $f(x)$ be continuous on $T_L$.  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 $\fhat(k)$ are
$$
	\fhat(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 \twopiL = \frac{2\pi}{N}.
$$
Making substitutions, we now have for the estimated integral
\begin{equation}
\label{eqn:dft}
	\begin{array}{ccc}
	\fhat(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 \twopiL 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{array}
\end{equation}

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.)

An advantage of this derivation of the discrete Fourier transform is that it
uses elementary and easily visualizable Riemann sums.  One asks, however, about
well-definedness.  Since $f \in C(T_L)$, and since $T_L$ is compact, $f \in
L^2(T_L, dx)$.  Likewise, $f \in L^1(T_L, dx)$.  Due to the continuity, we may
without loss of generality take $f$ to be a particular representative for an
equivalence class with respect to Lebesgue measure on $T_L$.

An alternative derivation of the discrete Fourier transform may be obtained by
defining 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.  This approach results in the same
DFT formula, namely, equation \ref{eqn:dft}.

%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.

%%% ----------------------------------------------------------------
%\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 = \fhat(k) = \fhat(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
%\begin{equation}
%\label{eqn:fourmxdft}
%	\vecc = \frac{1}{N}\, W \vecb.
%\end{equation}
%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.
%
%There are, to a degree depending on the prime factorization of $N$, redundant
%computations in equation \ref{eqn:fourmxdft}.  These are at their most
%pronounced when $N$ is a power of 2.  In practice, the \emph{fast Fourier
%transform} algorithm is used to reduce the computational complexity of the
%discrete Fourier transform.  See \cite{bib:nr} and \cite{bib:lyons} for
%information on this important topic.

%%% ================================================================
%\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{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*}
%	\fhat(\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*}
%	\fhat(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*}
%	\fhat(\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*}
%
%As with the 1-dimensional DFT, the computational complexity of the
%multi-dimensional DFT reduces drastically when the $N_r$'s consist of small
%prime factors.  See \cite{bib:nr} for more information on the multi-dimensional
%fast Fourier transform.

%%% ================================================================
%\section{Error analysis using the mean-value theorem}
%
%Following any elementary reference on numerical integration, e.g. \cite{bib:bf}
%or \cite{bib:wik}, suppose $h(x)$ is differentiable on $T_L$ with bounded
%derivative.  Let $M$ be that bound, i.e.
%$$
%	M = \sup_{x \in T_L} |h'(x)|.
%$$
%Use the $N$-point uniform mesh as before.  Then
%$$
%	\int_0^L h(x) dx \approx \sum_{j=0}^{N-1} h(x_j) \Delta x
%	= \sum_{j=0}^{N-1} h\left(\frac{jL}{N}\right) \frac{L}{N}.
%$$
%
%By the mean value theorem, for each $j$ and for each $x \in [x_j, x_{j+1})$,
%there exists a $\xi_x \in [x_j, x]$ such that
%$$
%	h'(\xi_x) = \frac{h(x) - h(x_j)}{x - x_j}
%$$
%i.e.
%$$
%	h(x) - h(x_j) = (x - x_j) h'(\xi_x).
%$$
%Integrating both sides and taking absolute values gives
%\begin{eqnarray*}
%	\left|
%		\int_{x_j}^{x_{j+1}} h(x) dx - (x_{j+1} - x_j) h(x_j)
%	\right|
%	&=& \left|
%			\int_{x_j}^{x_{j+1}} (x-x_j) h'(\xi_x) dx
%		\right| \\
%	&\le& \int_{x_j}^{x_{j+1}} \left| (x-x_j) h'(\xi_x) \right| dx \\
%	&\le& \int_{x_j}^{x_{j+1}} (x-x_j) M dx \\
%	&=& M\frac{(x_{j+1} - x_j)^2}{2} \\
%	&=& \frac{ML^2}{2N^2}.
%\end{eqnarray*}
%
%Now take the above estimate for the error on the $j$th interval.  Summing these
%up and using the triangle inequality yields an upper bound for the error on the
%Riemann sum, namely,
%
%\begin{eqnarray*}
%	\left|
%		\int_0^L h(x) dx - \sum_{j=0}^{N-1} (x_{j+1} - x_j) h(x_j)
%	\right|
%	&=&
%	\left|
%		\sum_{j=0}^{N-1} \left[
%			\int_{x_j}^{x_{j+1}} h(x) dx - (x_{j+1} - x_j) h(x_j)
%		\right]
%	\right| \\
%	&\le& \sum_{j=0}^{N-1}
%		\left|
%			\int_{x_j}^{x_{j+1}} h(x) dx - (x_{j+1} - x_j) h(x_j)
%		\right| \\
%	&\le& \sum_{j=0}^{N-1} \frac{ML^2}{2N^2}
%	= \frac{ML^2}{2N}.
%\end{eqnarray*}
%
%This makes sense intuitively:  a constant function is exactly approximated by
%any Riemann sum; steeper slopes may contribute errors.  Note that one may
%further refine the error estimate, if higher degrees of differentiability are
%assumed.
%
%Now recall that we are taking $h(x) = e^{-ikx}f(x)$ for $k \in \ZtwopiL$.  Then
%\begin{eqnarray*}
%	h'(x) &=& e^{-ikx}(f'(x) - ikf(x)) \\
%	|h'(x)| &=& |f'(x) - ikf(x)| \ge |f'(x)| - |k||f(x)|
%\end{eqnarray*}
%such that in general (barring any preconceptions about $f(x)$), error worsens
%for higher $|k|$.

%% ================================================================
\section{Poisson summation on the real line}

In section \ref{sec:poierr}, Poisson summation will be used to analyze
approximation error in the DFT.  First, however, the theory will be examined on
the real line, to fix concepts.  The presentation here elaborates on section
24.6 of \cite{bib:faris}.

\begin{prop}
Let $L$ be a positive real number.  Let $f(x) \in L^1(\R, dx) \cap L^2(\R, dx)$
with Fourier transform $\fhat \in L^1(\R, dk/2\pi)$.  Further suppose that
$$
	\sum_{k \in \ZtwopiL} \fhat(k) < +\infty.
$$
Then
$$
	\sum_{x \in \Z L}f(x) = \frac{1}{L} \sum_{k \in \ZtwopiL} \fhat(k).
$$
\begin{proof}
Recall that for $k \in \R$ we have
$$
	\fhat(k) = \int_{-\infty}^{+\infty} e^{-ikx}f(x) dx.
$$
Now define (this utility function being the central trick of the proof)
$$
	S(t) = \sum_{x \in \Z L} f(x+t).
$$
This sum is finite since $f$ is in $L^1(\R, dx)$.  Since $S(t)$ is well-defined
on $\R/\Z L = T_L$, it has Fourier coefficients
$$
	a_k = \frac{1}{L} \int_0^L e^{-ikt} S(t) dt
$$
for $k \in \ZtwopiL$.
%Since [xxx need to justify this --- I think we need $S \in L^2(T_L, dt)$?],
We may reconstruct $S(t)$ via
\begin{equation}
	S(t) = \sum_{k \in \ZtwopiL} a_k e^{ikt}.
\end{equation}
But
\begin{eqnarray*}
	a_k &=& \frac{1}{L} \int_0^L e^{-ikt} S(t) dt
	= \frac{1}{L} \int_0^L e^{-ikt} \left[\sum_{x \in \Z L} f(x+t)\right] dt.
\end{eqnarray*}
Since $f \in L^1(\R, dx)$, by Fubini's theorem we may interchange the integral
and sum to obtain
\begin{eqnarray*}
	a_k &=& \frac{1}{L} \sum_{x \in \Z L} \int_0^L e^{-ikt} f(x+t) dt.
\end{eqnarray*}
Since $x \in \Z L$, let $x = j L$ for $j \in \Z$.  Using change of variable $u
= jL+t$ in the integrand, and recalling that we are restricting $k$ to be in
$\ZtwopiL$, we have
\begin{eqnarray*}
	a_k &=& \frac{1}{L} \sum_{j=-\infty}^{+\infty}
		\int_{t=0}^{t=L} e^{-ikt} f(jL+t) dt\\
	&=& \frac{1}{L} \sum_{j=-\infty}^{+\infty} \int_{u=jL}^{u=(j+1)L}
		e^{-ik(u-jL)} f(u) du \\
	&=& \frac{1}{L} \sum_{j=-\infty}^{+\infty}
		\int_{u=jL}^{u=(j+1)L} e^{-iku} f(u) du \\
	&=& \frac{1}{L} 
		\int_{u=-\infty}^{u=+\infty} e^{-iku} f(u) du \\
	&=& \frac{\fhat(k)}{L}.
\end{eqnarray*}

Now simply evaluate $S(t)$ at $t=0$.  We have
$$
	S(t) = \sum_{x \in \Z L} f(x+t)
$$
and
$$
	S(t) = \sum_{k \in \ZtwopiL} a_k e^{ikt}
	= \frac{1}{L} \sum_{k \in \ZtwopiL} \fhat(k) e^{ikt}
$$
so
$$
	S(0) = \sum_{x \in \Z L} f(x)
	= \frac{1}{L} \sum_{k \in \ZtwopiL} \fhat(k).
$$
\end{proof}
\end{prop}

%% ================================================================
\section{Poisson summation on the circle}

The Poisson summation formula was worked out on the (familiar) real line in
order to give a concrete example for the more abstract setting.  Recall from
section \ref{sec:intro} that in this paper we are taking Pontryagin duals as
given.  In particular, if $G$ is an abelian group, write $\Ghat$ for its dual.
In the previous section, we had $G=\R$ with closed subgroup $H=\Z L$ and
quotient $G/H = T_L$.  Now, the dual group for $G$ is $\Ghat = \R$, and the
duals of $H$ and $G/H$ are as before in this paper, namely:
$$
\begin{array}{llllll}
\textbf{Group} &&& \textbf{Dual} && \\
G   &=& \R   & \Ghat  &=& \R \\
H   &=& \Z L & \Hhat  &=& T_{2\pi/L} \\
G/H &=& T_L  & \GHhat &=& \ZtwopiL \\
\end{array}
$$

%Given this terminology, we can walk back through the proof of the Poisson
%summation formula on the real line, which will suggest how to correctly write
%the summation formula for the circle.
%\begin{itemize}
%\item We took $f:G \to \R$, with $\fhat:\Ghat \to \R$.
%\item We required
%	$$\sum_{k\in\GHhat} \fhat(k) < +\infty.$$
%\item We claimed
%	$$\sum_{x\in H} f(x) = \frac{1}{2\pi} \sum_{k \in \GHhat} \fhat(k).$$
%\item Given $f:G \to \R$ and a measure $\mu$ on $G$, we had
%	$\fhat:\Ghat \to \R$ via
%	$$
%		\fhat(k) = \int_G e^{-ikx}f(x) d\mu.
%	$$
%\item Given $f(x):G\to \R$, we summed over cosets to obtain $S(t):G/H \to \R$:
%	$$
%		S(t) = \sum_{x \in H} f(x+t).
%	$$
%\item This yielded Fourier coefficients
%	$$
%		a_k = \int_{G/H} e^{-ikt} S(t) d\nu
%	$$
%	for $k \in \GHhat$ where $\nu$ is a measure on $G/H$.
%\item We reconstructed $S(t)$ via
%	\begin{equation}
%		S(t) = \sum_{k \in \GHhat} a_k e^{ikt}.
%	\end{equation}
%\item The Fourier transform of $S(t)$ was
%	\begin{eqnarray*}
%		a_k &=& \int_{G/H} e^{-ikt} S(t) dt \\
%		&=& \int_{t\in G/H} e^{-ikt} \left[\int_{x \in H} f(x+t) dx\right] dt
%	\end{eqnarray*}
%	for $k \in \GHhat$.
%\item Using Fubini's theorem to interchange the order of integrals, we
%	obtained
%	\begin{eqnarray*}
%		a_k &=& \int_{x \in H} \left[\int_{G/H} e^{-ikt} f(x+t) dt \right] dx.
%	\end{eqnarray*}
%\item Using a change of variable, and restricting $k$ from %	$\Ghat$ to $\GHhat$, we found
%	$$
%		a_k = \frac{1}{2\pi} \fhat(k).
%	$$
%\item Evaluating $S(t)$ at $t=0$ gave the desired result.
%\end{itemize}

Given this terminology, we can walk back through the proof of the Poisson
summation formula on the real line, which will suggest how to correctly write
the summation formula for the circle.  
We took $f:G \to \R$, with $\fhat:\Ghat \to \R$.
We required
$$
	\sum_{k\in\GHhat} \fhat(k) < +\infty,
$$
and we claimed that
$$
	\sum_{x\in H} f(x) = \sum_{k \in \GHhat} \fhat(k).
$$
The key to the proof was that $e^{-ikx}$ was well-defined on the quotient space
$\GHhat$, and that $G$ was the union of translates of $G/H$.  Now we simply
replace $G$, $H$, et al. by the following, using the information about $\Delta
x$ and $\Delta k$ obtained in section \ref{sec:riemann_sums}:
$$
\begin{array}{llllll}
\textbf{Group} &&& \qquad \textbf{Dual} && \\
	\\
G   &=& \R/\Z L   & \qquad \Ghat  &=& \Z\frac{2\pi}{L} \\
    &=&     T_L   & \qquad        & &    \\
	\\
H   &=& \Z\frac{L}{N} / \Z L & \qquad \Hhat  &=& \Z_N \frac{2\pi}{L} \\
    &=& \Z_N\left(\frac{L}{N}\right)  & & & \\
	\\
G/H &=& \frac{\R/\Z L}{\Z\frac{L}{N}/\Z L} & \qquad
		\GHhat &=& \ZtwopiNL \\
    &=& \R / \Z\frac{L}{N} & && \\
    &=& T_{L/N}            & && \\
\end{array}
$$

\begin{prop}
\label{prop:circlepoi}
Let $L$ be a positive real number.  Let $f(x) \in L^1(T_L, dx)$ with Fourier
transform $\fhat \in L^1(\ZtwopiL, dk/2\pi)$.  Further suppose that
$$
	\sum_{k \in \ZtwopiNL} \fhat(k) < +\infty.
$$
Then
$$
	\sum_{x \in \ZNLN}f(x) = N \sum_{k \in \ZtwopiNL} \fhat(k).
$$
\begin{proof}
For $k \in \ZtwopiL$ we have
$$
	\fhat(k) = \int_{0}^{L} e^{-ikx}f(x) dx.
$$
Define
$$
	S(t) = \sum_{x \in \ZNLN} f(x+t).
$$
Since $S(t)$ is well-defined on the quotient $T_{L/N}$, it has Fourier
coefficients $$
	a_k = \frac{N}{L} \int_0^{L/N} e^{-ikt} S(t) dt
$$
for $k \in \ZtwopiNL$.  We may reconstruct $S(t)$ via
\begin{equation}
	S(t) = \sum_{k \in \ZtwopiNL} a_k e^{ikt}.
\end{equation}
But
\begin{eqnarray*}
	a_k &=& \frac{N}{L} \int_0^{L/N} e^{-ikt} S(t) dt
	= \frac{N}{L} \int_0^{L/N} e^{-ikt}
		\left[ \sum_{x \in \ZNLN} f(x+t) \right] dt.
\end{eqnarray*}
Since $f \in L^1(T_L, dx)$, by Fubini's theorem we may again interchange the
integral and sum to obtain
\begin{eqnarray*}
	a_k &=& \frac{N}{L} \sum_{x \in \ZNLN} \int_0^{L/N} e^{-ikt} f(x+t) dt.
\end{eqnarray*}
Since $x \in \ZNLN$, let $x = j L/N$ for $j \in \Z_N$.  Using change of
variable $u = jL/N+t$ in the integrand, and recalling that we are restricting
$k$ to be in $\ZtwopiNL$, we have
\begin{eqnarray*}
	a_k &=& \frac{N}{L} \sum_{j=0}^{N-1}
		\int_{t=0}^{t=L/N} e^{-ikt} f(jL/N+t) dt\\
	&=& \frac{N}{L} \sum_{j=0}^{N-1} \int_{u=jL/N}^{u=(j+1)L/N}
		e^{-ik(u-jL/N} f(u) du \\
	&=& \frac{N}{L} \sum_{j=0}^{N-1}
		\int_{u=jL/N}^{u=(j+1)L/N} e^{-iku} f(u) du \\
	&=& \frac{N}{L} 
		\int_{u=0}^{u=L} e^{-iku} f(u) du \\
	&=& N \fhat(k).
\end{eqnarray*}

Evaluating $S(t)$ at $t=0$, we have
$$
	S(t) = \sum_{x \in \ZNLN} f(x+t)
$$
and
$$
	S(t) = \sum_{k \in \ZtwopiNL} a_k e^{ikt}
	= N \sum_{k \in \ZtwopiNL} \fhat(k) e^{ikt}
$$
so
$$
	S(0) = \sum_{x \in \ZNLN} f(x)
	= N \sum_{k \in \ZtwopiNL} \fhat(k).
$$
\end{proof}
\end{prop}

%% ================================================================
\section{The Poisson summation formula for the DFT}
\label{sec:poiformula}

Now apply proposition \ref{prop:circlepoi} not to $f(x)$ but rather to
$e^{-ipx}f(x)/N$.

A little lemma shows that a phase shift in the time domain is a translation in
the frequency domain.  To see this, compute:
$$
	\frac{1}{L} \int_0^L  e^{-ikx} \left[e^{-ipx} f(x) \right] dx =
	\frac{1}{L} \int_0^L e^{-i(k+p)x} f(x) dx =
	\fhat(k+p).
$$

By the lemma and proposition \ref{prop:circlepoi}, we have
\begin{equation}
\label{eqn:poiform1}
	\frac{1}{N} \sum_{x \in \ZNLN} e^{-ipx} f(x)
	= \sum_{k \in \ZtwopiNL} \fhat(k+p).
\end{equation}
The left-hand side is precisely the $N$-point discrete Fourier transform, at
$p$, of $f(x)$ on $T_L$.  For brevity, call this $\fhat_D(p)$.  Then equation
\ref{eqn:poiform1} becomes
\begin{equation}
\label{eqn:poiform2}
\framebox[6.0in][c]{
	$\displaystyle
	\fhat_D(p) = \sum_{k \in \ZtwopiNL} \fhat(k+p).$
}
\end{equation}

We may summarize the Poisson formula in the following way:  The $N$-point DFT
folds energy from all frequencies which are congruent mod $N$ into the same
bin.  Thus the error in the DFT is precisely the amount of energy at
frequencies outside the range $k=-N/2, \ldots, N/2-1$.  The Riemann-Lebesgue
lemma \cite{bib:faris} assures us that for $f \in L_1$, the Fourier
coefficients $\fhat(p)$ eventually go to zero.  However, as yet we lack a bound
on how fast they go to zero.  In particular, consider a sinusoidal input with
frequency $aN+1$, for any integer $a$.  Then the $N$-point DFT will fold all
the input's energy into $\fhat_D(1)$.  Yet, the true Fourier coefficient will
have its energy at $\fhat(aN+1)$.  This high-frequency signal is highly
oscillatory, though, which is to say that its second derivative takes on large
values.  This example motivates the following error analysis.

%% ================================================================
\section{First-order error analysis}
\label{sec:poierr}

Compare the right-hand side of equation \ref{eqn:poiform2} to the continuous
Fourier transform $$
	\fhat(p) = \frac{1}{L} \int_0^L e^{-ipx} f(x) dx
$$
for $p \in \ZtwopiL$.  The error in the DFT with respect to the CFT is
$$
	\left| \fhat(p) - \fhat_D(p) \right|
	= \left|
		\fhat(p) - \sum_{k \in \ZtwopiNL} \fhat(k+p)
	\right|
	= \left|
		\sum_{k \in \ZtwopiNL, k \ne 0} \fhat(k+p)
	\right|.
$$

\begin{prop}
The error in the DFT depends linearly on the upper bound on the second
derivative of the input.  Furthermore, the convergence of the DFT to the true
Fourier transform is quadratic in the number of mesh points.

\begin{proof}
Suppose $f(x)$ has bounded second derivative, i.e. let $M_2$ be a positive real
number such that for all $x \in T_L$,
\begin{equation}
	\label{eqn:bounded_deriv}
	|f''(x)| \le M_2.
\end{equation}
To relate this bound on the second derivative to a bound on $\fhat(k+p)$ for $k
\ne 0$, use the derivative property of the Fourier transform.  Letting $q=k+p$
for brevity, this is
\begin{equation}
	\label{eqn:fpiq}
	\fprimehat(q) = iq \fhat(q).
\end{equation}
A straightforward computation involving the transform $\fprimehat(q)$ yields a
harmonic series, which is divergent.  Note, however, that the square harmonic
series converges.  Thus, a slight trick is called for.  Since we are assuming
$f \in C(T_L)$, we have $f \in L^2(T_L, dx)$.  Following the technique of
theorem 23.11 of \cite{bib:faris}, put $1=(k+p)/(k+p)$ and use the
Cauchy-Schwarz inequality on $\ell^2$ to obtain
\begin{eqnarray*}
	\left| \sum_{k\ne 0, k \in \ZtwopiNL} \fhat(k+p) \right|
	&=& \left| \sum_{k\ne 0, k \in \ZtwopiNL}
			\frac{1}{k+p} (k+p) \fhat(k+p) \right| \\
	&\le&
		\sqrt{\sum_{k\ne 0, k \in \ZtwopiNL} \frac{1}{(k+p)^2}}
		\sqrt{\sum_{k\ne 0, k \in \ZtwopiNL} |(k+p)\fhat(k+p)|^2} \\
\end{eqnarray*}
The use of the inequality is justified since we take $f \in C(T_L)$, hence $f
\in L^2(T_L, dx)$, and therefore $\fhat \in L^2(T_L, dk/2\pi)$.

For the first term, since $k \in \ZtwopiNL$, take $k = \ell2\pi N/L$.
Then
\begin{eqnarray*}
	\sum_{k\ne 0, k \in \ZtwopiNL} \frac{1}{k^2}
	&=& \sum_{\ell\ne 0, \ell \in \Z} \left(\frac{L}{\ell2\pi N}\right)^2
	= \frac{2L^2}{4\pi^2N^2} \sum_{\ell=1}^\infty \frac{1}{\ell^2}
	= \frac{L^2}{2\pi^2N^2} \frac{\pi^2}{6}
	= \frac{L^2}{12N^2}.
\end{eqnarray*}

For the other term in the error estimate, by equation \ref{eqn:fpiq},
\begin{eqnarray*}
	\sum_{k\ne 0, k \in \ZtwopiNL} |(k+p)\fhat(k+p)|^2 &=&
	\sum_{k\ne 0, k \in \ZtwopiNL} |\fprimehat(k+p)|^2 =
	\sum_{k\ne 0, k \in \ZtwopiNL} \left|\frac{\fdprimehat(k+p)}{k+p}\right|^2\\
\end{eqnarray*}

To relate $\fdprimehat(k+p)$ to the known bound on $f''(x)$, we may use the
following lemma.  Let $g(x) \in C(T_L)$.  By the compactness of $T_L$, $g(x)$
is bounded: say by a positive real constant $M$.  Then $ |\ghat(q)| \le M.  $
To see this, compute
\begin{eqnarray*}
	|\ghat(q)| &=& \left| \frac{1}{L} \int_0^L e^{-iqx}g(x) dx \right|
	\le \frac{1}{L} \int_0^L |e^{-iqx}g(x)| dx
	= \frac{1}{L} \int_0^L |e^{-iqx}||g(x)| dx \\
	&=& \frac{1}{L} \int_0^L |g(x)| dx
	\le \frac{1}{L} \int_0^L M dx
	= M.
\end{eqnarray*}

Apply the lemma to the second term in the error estimate:
\begin{eqnarray*}
	\sum_{k\ne 0, k \in \ZtwopiNL} \left|\frac{\fdprimehat(k+p)}{k+p}\right|^2
	&\le& \sum_{k\ne 0, k \in \ZtwopiNL} \left|\frac{M_2}{k+p}\right|^2
	= M_2^2 \sum_{k\ne 0, k \in \ZtwopiNL} \left|\frac{1}{k+p}\right|^2 \\
	&\le& M_2^2 \sum_{k\ne 0, k \in \ZtwopiNL} \frac{1}{k^2}
	= \frac{M_2^2L^2}{12N^2}.
\end{eqnarray*}

Combining these results yields
\begin{eqnarray*}
	\left| \sum_{k\ne 0, k \in \ZtwopiNL} \fhat(k+p) \right|
	&\le&
		\sqrt{\sum_{k\ne 0, k \in \ZtwopiNL} \frac{1}{(k+p)^2}}
		\sqrt{\sum_{k\ne 0, k \in \ZtwopiNL} |(k+p)\fhat(k+p)|^2} \\
	&\le& \frac{L}{N\sqrt{12}}
		\sqrt{\sum_{k\ne 0, k \in \ZtwopiNL} |(k+p)\fhat(k+p)|^2} \\
	&=& \frac{L}{N\sqrt{12}}
		\sqrt{\sum_{k\ne 0, k \in \ZtwopiNL} |\fprimehat(k+p)|^2} \\
	&=& \frac{L}{N\sqrt{12}}
		\sqrt{\sum_{k\ne 0, k \in \ZtwopiNL}
			\left|\frac{\fdprimehat(k+p)}{k+p}\right|^2} \\
	&\le& \frac{LM_2}{N\sqrt{12}}
		\sqrt{\sum_{k\ne 0, k \in \ZtwopiNL} \frac{1}{k^2}} \\
	&=& \frac{L^2M_2}{12N^2}.
\end{eqnarray*}
\end{proof}
\end{prop}

%% ================================================================
\section{Higher-order error analysis}

Now suppose that $f$ is $m+1$ times differentiable.  Then its $(m+1)$st
derivative is continuous on the compact set $T_L$, and hence is in $L^2(T_L,
dx)$.  Furthermore, $f^{(m+1)}(x)$ is bounded, say by $M_{m+1}$.

\begin{prop}
The error in the DFT depends linearly on the upper bound on the $(m+1)$st
derivative of the input.  Furthermore, the convergence of the DFT to the true
Fourier transform is inversely proportional to $N^{m+1}$ where $N$ is the
number of mesh points.

\begin{proof}
We may repeatedly apply equation \ref{eqn:fpiq} to obtain
\begin{eqnarray*}
	\left| \sum_{k\ne 0, k \in \ZtwopiNL} \fhat(k+p) \right|
	&\le& \frac{L}{N\sqrt{12}}
		\sqrt{\sum_{k\ne 0, k \in \ZtwopiNL} |\fprimehat(k+p)|^2} \\
	&\le& \frac{L}{N\sqrt{12}}
		\sqrt{\sum_{k\ne 0, k \in \ZtwopiNL}
		\left|\frac{\widehat{f^{(m+1)}}(k+p)}{(k+p)^{m}}\right|^2} \\
	&\le& \frac{LM_{m+1}}{N\sqrt{12}}
		\sqrt{\sum_{k\ne 0, k \in \ZtwopiNL} \frac{1}{k^{2m}}} \\
	&=& \frac{LM_{m+1}}{N\sqrt{12}}
		\sqrt{\sum_{\ell\ne 0, \ell \in \Z}
		\left(\frac{L}{\ell 2\pi N}\right)^{2m}} \\
	&=& \frac{L^{m+1}M_{m+1}}{N^{m+1}(2\pi)^m\sqrt{12}}
		\sqrt{\sum_{\ell\ne 0, \ell \in \Z}
		\frac{1}{\ell^{2m}}} \\
	&=& \frac{L^{m+1}M_{m+1}}{N^{m+1}(2\pi)^m}
		\sqrt{2 \sum_{\ell=1}^{\infty} \frac{1}{\ell^{2m}}}.
\end{eqnarray*}

The sum inside the square root is the Riemann zeta function at positive
even integers, which is well-known to be
$$
	\sum_{\ell=1}^{\infty} \frac{1}{\ell^{2m}}
	= \zeta(2m) = \frac{2^{2m-1} |B(2m)| \pi^{2m}}{(2m)!}
$$
where $B(2m)$ is the $2m$th Bernoulli number.  The first few values are
\cite{bib:mathworld} $\zeta(2)=\pi^2/6$, $\zeta(4)=\pi^4/90$,
$\zeta(6)=\pi^6/945$, and $\zeta(8)=\pi^8/9450$.  Detailed discussion of the
Riemann zeta function is beyond the scope of this paper; for the present
purpose, it suffices to know that, for even positive integers $2m$, $\zeta(2m)$
is bounded above by $\pi^{2m}$.  Also, $\sqrt{12} < 4$.  Thus
$$
	\sqrt{2 \sum_{\ell=1}^{\infty} \frac{1}{\ell^{2m}}}
	\le \sqrt{\frac{2^{2m}\pi^{4m}}{(2m)!}}
	= 2^m \pi^{2m} \sqrt{\frac{1}{(2m)!}}
$$
Thus an error bound in the case $f$ is $m+1$ times differentiable is
$$
	\left(\frac{L}{N}\right)^{m+1}
	\frac{M_{m+1} 2^m \pi^{2m}} {4 (2\pi)^m \sqrt{(2m)!}}
	= \left(\frac{L\pi}{N}\right)^{m+1} \frac{M_{m+1}}{4\pi \sqrt{{(2m)!}}}.
$$
\end{proof}
\end{prop}

%Here, we see that the DFT converges quite rapidly to the CFT as $N$ increases.

%To estimate the asymptotic error behavior as $m$ increases, we may use
%Stirling's approximation for the factorial function, namely,
%$$
%	n! \approx \left(\frac{n}{e}\right)^n \sqrt{2\pi n}.
%$$
%Then
%$$
%	\frac{1}{(2m)!} \approx \left(\frac{e}{2m}\right)^m \sqrt[4]{4\pi m}.
%	= \left(\frac{e}{2m}\right)^m \sqrt[4]{4\pi} m^{1/4}.
%$$
%We may now write the error bound as
%$$
%	\left(\frac{L\pi e}{2 N}\right)^{m+1} \frac{\sqrt[4]{4\pi} }{2\pi e}
%		\frac{M_{m+1}}{m^{m-1/4}}.
%$$

%We conclude that, given a bound on the $(m+1)$st derivative of $f(x)$,
%convergence of the DFT to the CFT is exponential in $N$.
