\documentclass{slides}
\usepackage[landscape]{geometry}
\parindent=0pt \parskip=0pt

\usepackage{amssymb,amsmath,amsthm}
%\usepackage{psfig}
\usepackage{graphics}
\input{../tex/isomath}
\input{./syms}

%% ================================================================
\begin{document}

%% ----------------------------------------------------------------
\begin{slide}
\begin{center}
\textbf{Poisson summation and the discrete Fourier transform} \\
John Kerl \\
University of Arizona, Math 523B \\
April 14, 2006
\end{center}

\end{slide}
%% ----------------------------------------------------------------
\begin{slide}
\begin{center}
\end{center}
In Euclidean space, given a vector, \ldots

\setlength{\unitlength}{1.0in}
\begin{picture}(7.0,3.2)
%\put(0.0,0.0){\framebox(7.0,3.2)}
\thicklines
\put(0.2,0.2){\vector(4,3){4.0}}
\end{picture}

\end{slide}
%% ----------------------------------------------------------------
\begin{slide}
\begin{center}
\end{center}
\ldots we can put down a coordinate frame (say an orthonormal one).  We can
\textbf{decompose} the vector into its components, or \textbf{reconstruct} the
vector from its components.

\setlength{\unitlength}{1.0in}
\begin{picture}(7.0,3.2)
%\put(0.0,0.0){\framebox(7.0,3.2)}
\put(0.0,1.5){3}
\put(2.0,-0.1){4}
\thicklines
\put(0.2,0.2){\vector(4,3){4.0}}
\put(0.2,0.2){\vector(1,0){1.0}}
\put(0.2,0.2){\vector(0,1){1.0}}
\thinlines
\put(0.2,0.2){\vector(1,0){4.0}}
\put(0.2,0.2){\vector(0,1){3.0}}
%
\put(4.2,0.2){\line(0,1){0.1}}
\put(4.2,0.4){\line(0,1){0.1}}
\put(4.2,0.6){\line(0,1){0.1}}
\put(4.2,0.8){\line(0,1){0.1}}
\put(4.2,1.0){\line(0,1){0.1}}
\put(4.2,1.2){\line(0,1){0.1}}
\put(4.2,1.4){\line(0,1){0.1}}
\put(4.2,1.6){\line(0,1){0.1}}
\put(4.2,1.8){\line(0,1){0.1}}
\put(4.2,2.0){\line(0,1){0.1}}
\put(4.2,2.2){\line(0,1){0.1}}
\put(4.2,2.4){\line(0,1){0.1}}
\put(4.2,2.6){\line(0,1){0.1}}
\put(4.2,2.8){\line(0,1){0.1}}
\put(4.2,3.0){\line(0,1){0.1}}
%
\put(0.2,3.2){\line(1,0){0.1}}
\put(0.4,3.2){\line(1,0){0.1}}
\put(0.6,3.2){\line(1,0){0.1}}
\put(0.8,3.2){\line(1,0){0.1}}
\put(1.0,3.2){\line(1,0){0.1}}
\put(1.2,3.2){\line(1,0){0.1}}
\put(1.4,3.2){\line(1,0){0.1}}
\put(1.6,3.2){\line(1,0){0.1}}
\put(1.8,3.2){\line(1,0){0.1}}
\put(2.0,3.2){\line(1,0){0.1}}
\put(2.2,3.2){\line(1,0){0.1}}
\put(2.4,3.2){\line(1,0){0.1}}
\put(2.6,3.2){\line(1,0){0.1}}
\put(2.8,3.2){\line(1,0){0.1}}
\put(3.0,3.2){\line(1,0){0.1}}
\put(3.2,3.2){\line(1,0){0.1}}
\put(3.4,3.2){\line(1,0){0.1}}
\put(3.6,3.2){\line(1,0){0.1}}
\put(3.8,3.2){\line(1,0){0.1}}
\put(4.0,3.2){\line(1,0){0.1}}
\end{picture}
\end{slide}
%% ----------------------------------------------------------------
\begin{slide}
\begin{center}
\end{center}
Decomposition (\emph{how much in each direction}) in $\R^n$:  let $\{\vecu_1,
\ldots, \vecu_n\}$ be an orthonormal basis.  Our vector is $\vecv$:
$$
	\vecv = \sum_{j=1}^n c_j \vecu_j
$$

How do we find the coefficients?  Since the basis vectors are
orthonormal, $\langle \vecu_i, \vecu_j \rangle = \delta_{ij}$.  Just dot each
basis vector with $\vecv$:
$$
	\langle \vecu_i, \vecv \rangle = \langle \vecu_i, \sum_{j=1}^n c_j \vecu_j \rangle
	= \sum_{j=1}^n c_j \langle \vecu_i, \vecu_j \rangle
	= c_i.
$$

\end{slide}
%% ----------------------------------------------------------------
\begin{slide}
Recall that in $\C^n$, we must conjugate one argument (by our convention, the
first) in order to get positive definiteness of the inner product:
$$
	\langle \vecv, \vecw \rangle = 
	\sum_{j=1}^n \overline{v}_j w_j.
$$

\end{slide}
%% ----------------------------------------------------------------
\begin{slide}
Consider the space of periodic continuous functions from the interval $[0,L]$
(or alternatively from the circle $T_L$) to the complex numbers.  (Continuous
on compact domain $\implies$ both $L_1$ and $L_2$.) Try to form an orthonormal
basis using $e^{ikx}$ as basis functions.  Recall that for functions from $\R$
to $\C$, $k$ can take any real value.  Here, what values can $k$ take?  For
example, look at $k=3$ and $k=3.2$:

\includegraphics{goodkbadk.eps}

\ldots so, periodicity requires $k = 2\pi\ell/L$ for $\ell \in \Z$.

\end{slide}
%% ----------------------------------------------------------------
\begin{slide}
Inner product of two functions $f(x)$ and $g(x)$ is
$$
	\langle f(x), g(x) \rangle = \frac{1}{L} \int_0^L \overline{f}(x) g(x) dx.
$$
The $e^{ikx}$ functions are readily seen to be orthonormal; it may be shown
that they span our space as well.

Decomposition (\emph{how much energy at each frequency}) of a function in terms
of the basis functions:
$$
	\fhat(k) = 
	\langle e^{ikx}, f(x) \rangle = \frac{1}{L} \int_0^L e^{-ikx} f(x) dx.
$$
This is a \bfidx{Fourier transform}.  Often called a \bfidx{Fourier series}
when the domain is $T_L$:  it maps from $L^2(T_L, dx)$ to $\ell^2(2\pi\Z/L)$.
\end{slide}
%% ----------------------------------------------------------------
\begin{slide}
The previous integral is over infinitely (uncountably many) $x$ in $T_L$.
\textbf{Approximate} it using a Riemann sum on a uniform mesh of $N$ points:
\begin{eqnarray*}
	\fhat(k) &=& \frac{1}{L} \int_0^L e^{-ikx} f(x) dx
	\approx \frac{1}{L} \sum_{j=1}^N e^{-ikx} f(x_j) \Delta x \\
	&=& \frac{1}{L} \sum_{j=1}^N e^{-ikjL/N} f(jL/N) \frac{L}{N}
	= \frac{1}{N} \sum_{j=1}^N e^{-ikjL/N} f(jL/N).
\end{eqnarray*}
This latter is called the \textbf{discrete Fourier transform}:
$$
	\fhat_D(k) = \frac{1}{N} \sum_{j=1}^N e^{-ikjL/N} f(jL/N).
$$
We have $\Delta x = L/N$ and $\Delta k = 2\pi/L$ so $\Delta x \Delta k = 2\pi/N$.
\end{slide}
%% ----------------------------------------------------------------
\begin{slide}
What are the distinct values of $k$ now?  Central point of this paper (to be
quantifed below in terms of Poisson summation):  when you sample \textbf{only}
at $N$ points, you \textbf{cannot} distinguish a frequency of $k$ from a
frequency of $N+k$.  For example, with $N = 8$, here are $\sin(1 \cdot 2\pi
x/L)$ and $\sin(9 \cdot 2\pi x/L)$:

\includegraphics{fig_9_1.eps}

\end{slide}
%% ----------------------------------------------------------------
\begin{slide}

Now look just at the 8 sampled points.  What is the true frequency?

\includegraphics{fig_9_1_bare.eps}

Result:  There are $N$ different $x$'s  with $x_j = jL/N$; there are $N$
different $k$'s  with $k = 2\pi\ell/L$.  Both $j$ and $\ell$ are taken mod $N$.

\end{slide}
%% ----------------------------------------------------------------
\begin{slide}
\begin{center}
\textbf{Duality}
\end{center}

\begin{itemize}
\item
Functions from $\R$ (Fourier transform) to $\C$:  $x \in \R$; $k \in \R$.  The
dual of $\R$ is $\R$.

\item
Functions from $T_L$ (Fourier series) to $\C$: dual of $T_L$ is $\Z 2\pi/L$.

\item
Functions from $\Z_N \left(\frac{L}{N}\right)$ to $\C$ (discrete Fourier
transform): dual of $\Z_N L/N$ is $\Z_N 2\pi/L$.
\end{itemize}

Pontryagin duality is an interesting subject on its own; we only touch on it
here.
\end{slide}
%% ----------------------------------------------------------------
\begin{slide}
\begin{center}
\textbf{Poisson summation formula}
\end{center}
Intuition:  As seen above, the $N$-point DFT \textbf{folds} energy from
\textbf{all} frequencies which are congruent mod $N$ into the \textbf{same}
bin.

The formula:
$$
	\fhat_D(p) = \sum_{k \in \ZtwopiNL} \fhat(k+p).
$$

Example:  Take $N=8$.  Then the discrete Fourier transform for frequency 1, in
terms of the true transform values, is
$$\fhat_D(1) = \ldots + \fhat(-7) + \fhat(1) + \fhat(9) + \fhat(17) + \ldots.$$
\end{slide}
%% ----------------------------------------------------------------
\begin{slide}

\textbf{Proposition.}  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).
$$

\emph{Proof.}  Mimic Faris' proof for Poisson summation on the real line.
See term paper for full details.
\hfill\qedsymbol

%\end{slide}
%%% ----------------------------------------------------------------
%\begin{slide}
%
%\emph{Proof.}  For $k \in \ZtwopiL$ we have
%$$
%	\fhat(k) = \int_{0}^{L} e^{-ikx}f(x) dx.
%$$
%
%\end{slide}
%%% ----------------------------------------------------------------
%\begin{slide}
%
%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*}
%
%\end{slide}
%%% ----------------------------------------------------------------
%\begin{slide}
%
%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).
%$$
%$\hfill\qedsymbol$

\end{slide}
%% ----------------------------------------------------------------
\begin{slide}

\textbf{Lemma.}
A phase shift in the time domain is a translation in
the frequency domain.

\emph{Proof.}
$$
	\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).
$$
\hfill\qedsymbol

By the lemma, and applying the proposition not to $f(x)$ but rather to
$e^{-ipx}f(x)/N$, we have
\begin{equation*}
	\frac{1}{N} \sum_{x \in \ZNLN} e^{-ipx} f(x)
	= \fhat_D(p)
	= \sum_{k \in \ZtwopiNL} \fhat(k+p).
\end{equation*}
This proves the Poisson summation formula.
\end{slide}
%% ----------------------------------------------------------------
\begin{slide}
\begin{center}
\textbf{First-order error analysis}
\end{center}

Intuition: The error in the DFT is the amount of energy
at frequencies outside the range $k=0, 1, \ldots, N-1$.

By the Riemann-Lebesgue lemma, for $f \in L_1$, the Fourier coefficients
$\fhat(p)$ eventually go to zero.  But how \textbf{fast}?

In particular, we can construct 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 \textbf{oscillatory},
though, which is to say that its second derivative takes on large values.

\end{slide}
%% ----------------------------------------------------------------
\begin{slide}

Quantitatively, the error in the DFT with respect to the CFT is
\begin{eqnarray*}
	\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|.
\end{eqnarray*}

\end{slide}
%% ----------------------------------------------------------------
\begin{slide}

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

\emph{Proof.} Let $M_2$ be a positive real number such that for all $x \in
T_L$,
$$
	\label{eqn:bounded_deriv}
	|f''(x)| \le M_2.
$$
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:
$$
	\fprimehat(k+p) = i(k+p) \fhat(k+p).
$$
A straightforward computation involving the transform $\fprimehat(k+p)$ 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 Faris' notes, 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,
\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.

\end{slide}
%% ----------------------------------------------------------------
\begin{slide}

\textbf{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 for all $q \in \Z2\pi/L$, $|\ghat(q)|$ is
also bounded by $M$.

\emph{Proof.}
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 = \frac{ML}{L} = M.
\end{eqnarray*}
\hfill\qedsymbol

\end{slide}
%% ----------------------------------------------------------------
\begin{slide}

Apply the lemma to the 2nd 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*}
\hfill\qedsymbol

\end{slide}
%% ----------------------------------------------------------------
\end{document}
