\documentstyle[psfig]{siam}
%\voffset=-0.75in
%\hoffset=-0.5in
%\textheight=10in
%\textwidth=7.05in

%
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\bear}{\begin{eqnarray}}
\newcommand{\bears}{\begin{eqnarray*}}
\newcommand{\eear}{\end{eqnarray}}
\newcommand{\eears}{\end{eqnarray*}}
\newcommand{\bull}{\item[$\bullet$]}
\newcommand{\eqdef}{\stackrel{\rm def}{=}}
\newcommand{\pf}{\noindent{\sc Proof: }}
\newcommand{\QED}{$\Box$}

\newcommand{\polys}{{\cal P}}
\newcommand{\trigs}{{\cal T}}
\newcommand{\splines}{{\cal S}}
\newcommand{\R}{{\bf R}}
\newcommand{\C}{{\bf R}}
\newcommand{\T}{{\bf T}}
%
\newcommand{\sA}{{\cal A}}
\newcommand{\sN}{{\cal N}}
\newcommand{\bfR}{{\bf R}}
\newcommand{\bfx}{{\bf x}}
\newcommand{\showlab}[1]{\label{#1}\mbox{ \tiny{\{#1\}}\ \ \ }}
%\newcommand{\showlab}[1]{\label{#1}}

%\newtheorem{corollary}{Corollary}
%\newtheorem{definition}{Definition}
%\newtheorem{exercise}{Exercise}
%\newtheorem{lemma}{Lemma}
%\newtheorem{theorem}{Theorem}

\begin{document}
%\addtolength{\baselineskip}{0.5\baselineskip}
%
% ------------------------------------------------------
\title{Math 669 Notes}
\author{J. G. Wade}
\date{\today}
\maketitle
% ---------------------------------------------------
\begin{abstract}
  These notes cover the period of the course around the end of October and
early November. We give a detailed derivation of the ``Discrete Fourier
Transform'' from the point of view of wanting to use a ``Gaussian quadrature''
type scheme to approximate the integrals which define the 
the Fourier series coefficients. Thus we again interpret what we're doing as
{\em (a)} approximating a function $f$ via the Projection Theorem, and then
{\em (b)} 
using a ``natural'' quadrature scheme (based on the roots of the underlying
orthogonal basis functions) 
to compute the inner product integrals arising from the Projection Theorem.
\end{abstract}
% --------------------------------------------------

%*****************************************************************************
\section{Discrete Fourier Transform}

  For $f\in L^2(-\pi,\pi)$ we have the Fourier series
\bear \showlab{f.fourier}
   f(x) &=& \sum_{k=-\infty}^\infty \hat{f}_k e^{ikx}, \\
   \hat{f}_k &=& {1\over 2\pi} \int_{-\pi}^\pi f(x) e^{-ikx}\,dx.  
   \showlab{fhat.k}
\eear
The problem we're dealing with in this section is that of deriving a
practically-implementable method for approximating the integral in the
formula for $\hat{f}_k$.
This combined with truncation of the series in (\ref{f.fourier}) 
(which is orthogonal projection, as we've discussed in detail)
will yield an effective and practical method for approximating $f$, 
assuming that $f$ is uniformly bounded. 

Denote by $\Omega$ the unit circle in the complex plane. 
We make the following change of variables: for $x\in [-\pi,\pi]$, we 
set $\omega\in \Omega$, $\omega\eqdef e^{ix}$. Then 
$dx = (1/i\omega)d\omega$, and the basis function
$e^{ikx}$ appearing in (\ref{f.fourier}) may be written as 
$ e^{ikx} = \omega^k$.
Then, with 
$
  g(\omega) \eqdef f(-i\log\omega),
$
we have that (\ref{f.fourier}) is equivalent to 
\bear \nonumber
   g(\omega) &=& \sum_{k=-\infty}^\infty  \hat{f}_k \omega^k \\
    \hat{f}_k &=& {1\over i2\pi}\int_\Omega 
                  g(\omega)\,\omega^{-k}\;{1\over\omega}d\omega.
   \showlab{fhat.g}
\eear

In analogy with 
othogonal polynomials on $[-1,1]$, we seek a ``Gaussian quadrature''
method based on, say,  $m$ points ($m$ depending on $n$), such that 
\[
   \sum_{j=0}^{m-1} 
   h(\rho^{(m)}_j)\gamma^{(m)}_j = \int_\Omega h(\omega)\;{1\over\omega}d\omega
\]
for $h$ a polynomial in $\omega$ of as high a degree as possible. 

\begin{lemma} \showlab{dft.alias}
  Let $m$ be a positive integer and let 
  $\rho^{(m)}_j$ be the ``$m$ roots of unity'', that is, 
  \[
     \rho^{(m)}_j \eqdef e^{i2\pi j/m}, \quad j=0,\ldots,m-1
  \]
  (so that $(\rho^{(m)}_j)^m=1$ for $j=0,\ldots,m-1$).
  Let $h$ be continuous on $\Omega$ and of the form
  \be \showlab{h.form}
    h(\omega) = \sum_{\sigma=-\infty}^\infty  \hat{h}_\sigma \omega^\sigma
  \ee
  where $\{\hat{h}_\sigma\}$ is absolutely summable. Then, 
  \be  \showlab{h.int.exact} 
   {1\over m}\sum_{j=0}^{m-1} 
   h(\rho^{(m)}_j) = \sum_{l=-\infty}^\infty \hat{h}_{lm}
  \ee
\end{lemma}
\pf
Setting 
\be \showlab{omega.m.def}
   \omega_m \eqdef e^{i 2\pi/m}, 
\ee
we have
\bears
   \sum_{j=0}^{m-1} h(\rho^{(m)}_j) 
    &=& \sum_{j=0}^{m-1}\sum_{\sigma=-\infty}^\infty  
                  \hat{h}_\sigma (\rho^{(m)}_j)^\sigma \\
    &=& \sum_{j=0}^{m-1}\sum_{l=-\infty}^\infty\sum_{s=0}^{m-1}
                  \hat{h}_{(lm+s)} (\rho^{(m)}_j)^{lm+s} \\
    &=& \sum_{j=0}^{m-1}\sum_{l=-\infty}^\infty\sum_{s=0}^{m-1}
                  \hat{h}_{(lm+s)} (\rho^{(m)}_j)^{s} \\
    &=& \sum_{j=0}^{m-1}\sum_{l=-\infty}^\infty\sum_{s=0}^{m-1}
                  \hat{h}_{(lm+s)} \omega_m^{js}
\eears
where in the second to last line we have used the fact that 
$(\rho^{(m)}_j)^m = 1$ and
in the last line we have used the fact that 
$(\rho^{(m)}_j)^s = (e^{i2\pi j/m})^s = (e^{i2\pi /m})^{js} = 
\omega_m^{js}$. Interchanging the sums we have
\bears
    \sum_{j=0}^{m-1} h(\rho^{(m)}_j) 
    &=& \sum_{l=-\infty}^\infty\sum_{s=0}^{m-1} \hat{h}_{(lm+s)} 
        \sum_{j=0}^{m-1} \omega_m^{js} \\
    &=& \sum_{l=-\infty}^\infty \left\{
         m\hat{h}_{lm} + \sum_{s=1}^{m-1} \hat{h}_{(lm+s)} 
         \sum_{j=0}^{m-1} (\omega_m^s)^j \right\} \\
    &=& \sum_{l=-\infty}^\infty \left\{
         m\hat{h}_{lm} 
         + \sum_{s=1}^{m-1} \hat{h}_{(lm+s)} 
           \Bigl[ {\omega_m^{sm} - 1 \over
                   \omega_m^s    - 1 }\Bigr]\;\right\} \\
    &=& \sum_{l=-\infty}^\infty 
         m\hat{h}_{lm}.
\eears
This yields (\ref{h.int.exact}). \QED

Now, 
\be \showlab{h0.1}
  {1\over i2\pi} \int_\Omega h(\omega) \;{1\over \omega}d\omega = \hat{h}_0.
\ee
This can also be verified by formally integrating (\ref{h.form})
term-by-term and using the fact that
\[
       {1\over i2\pi}\int_\Omega \omega^k \;{1\over\omega}d\omega 
   \;\;=\;\; {1\over  2\pi}\int_{-\pi}^{\pi} e^{ikx} 
   \;\;=\;\; \delta_{k,0}.
\]
(Compare (\ref{h0.1}) to the statement of the Cauchy Residue Theorem.)
In light of (\ref{h.int.exact}) and (\ref{h0.1}), we have that
\be  \showlab{quad.error}
     {1\over m}\sum_{j=0}^{m-1} h(\rho^{(m)}_j)
   \;-\;\; {1\over i2\pi} \int_\Omega h(\omega)\;{1\over \omega}d\omega \;=\;
    \sum_{l=-\infty\atop l\neq 0}^{\infty} \hat{h}_{lm}.
\ee

In our situation, $h(\omega) = f(-i\log(\omega)\bigr)\omega^{-k}$ for 
some $k$. Substituting this into the left side of (\ref{quad.error}) 
and using the fact that 
$f(-i\log(\rho^{(m)}_j)\,\bigr)
=f(-i\log(e^{i2\pi j/m})\,\bigr)
=f(2\pi j/m)$ and the fact that 
$(\rho^{(m)}_j)^{-k} = \omega_m^{-jk}$, 
 we obtain
\bears
   &&  {1\over m}\sum_{j=0}^{m-1} 
       f(-i\log(\rho^{(m)}_j)\,\bigr)(\rho^{(m)}_j)^{-k}
      \;-\; 
      {1\over i2\pi} \int_\Omega 
       f(-i\log(\omega)\,\bigr)\omega^{-k} \;{1\over \omega}d\omega \\
   &=& {1\over m}\sum_{j=0}^{m-1} 
       f(2\pi j/m)\omega_m^{-jk}
      \;-\; 
      {1\over 2\pi} \int_\Omega 
       f(x)e^{-ikx}\;dx \\
   &=& {1\over m}\sum_{j=0}^{m-1} 
       f(2\pi j/m)\omega_m^{-jk}
      \;-\; \hat{f}_k.
\eears
On the other hand, since we're identifying $h$ with 
$f(-i\log(\omega)\,\bigr)\omega^{-k}$, we have 
$\hat{h}_\sigma = \hat{f}_{\sigma+k}$, so the
right side of (\ref{quad.error}) becomes
\[
    \sum_{l=-\infty\atop l\neq 0}^{\infty} \hat{h}_{lm}
   = \sum_{l=-\infty\atop l\neq 0}^{\infty} \hat{f}_{lm+k}.
\]
Thus (\ref{quad.error}) is equivalent to 
\be  \showlab{alias.error}
     {1\over m}\sum_{j=0}^{m-1} 
       f(x^{(m)}_j)\omega_m^{-jk}
      \;=\; \hat{f}_k + \sum_{l=-\infty\atop l\neq 0}^{\infty} \hat{f}_{lm+k},
\ee
where 
\be  \showlab{x.m.def}
    x^{(m)}_j \eqdef {2\pi j\over m}
\ee
and $\omega_m$ is given by (\ref{omega.m.def}).
We summarize this discussion in the following theorem.
\begin{theorem} \showlab{thm:alias}
  For a given positive integer $m$, 
  with $x^{(m)}_j$ given by (\ref{x.m.def})
  and $\omega_m$ given by (\ref{omega.m.def}), 
  we have (\ref{alias.error}).
\end{theorem}

   The left side of (\ref{alias.error}) may be regarded as a 
numerical quadrature formula for approximating $\hat{f}_k$ --- in fact, 
it is nothing other than the ``Trapezoidal Rule''! With this interpretation, 
the second term on the right of (\ref{alias.error}) is an error term. It is
commonly referred to as the ``aliasing error'' since all of the coefficients
$\hat{f}_s, \; s=k\;{\rm mod}\; m$ are ``aliased to'' (added to) $\hat{f}_k$.
The aliasing error will be central in some of our subsequent discussions.

  Now, as stated at the top, the goal is to arrive at a practical and 
effective method for approximating $f$ by, say, $f_n$, by
approximately computing the orthogonal projection of $f$ upon $\trigs_n$.
We say ``by {\em approximately} computing the orthogonal projection''
because we're going to approximate $\hat{f}_k$ by (\ref{alias.error}).

\newpage
\noindent To summarize, \underline{Here's the method:}
\begin{quote}
\hrulefill\newline
For a given $n$, we set $m\eqdef 2n+1$, and
\be \showlab{fn.def}
    f_n(x) \eqdef \sum_{k=-n}^n \tilde{f}^{(n)}_k e^{ikx}
\ee
where
\be \showlab{fn.tilde}
  \tilde{f}^{(n)}_k \eqdef
  {1\over m}\sum_{j=0}^{m-1}  f(x^{(n)}_j)\omega_n^{-jk},
\ee
and $\omega_m$ and $x^{(m)}_j$ are 
given by (\ref{omega.m.def}) and (\ref{x.m.def}), respectively.
The $\tilde{f}^{(n)}_k$'s are related to the $\hat{f}_k$'s through
the ``aliasing error'' 
\be \showlab{alias.error2}
     \tilde{f}^{(n)}_k - \hat{f}_k  = 
      \sum_{l=-\infty\atop l\neq 0}^{\infty} \hat{f}_{lm+k}.
\ee
\hrulefill
\end{quote}

\bigskip
   

Note that if we set the notation
\bears
   \vec{f}^{(m)} &=& \{f(x^{(m)}_j)\}_{j=0}^{m-1} \\
   \tilde{f}^{(n)} &=& \{\tilde{f}^{(n)}_k\}_{k=-n}^{n}, 
\eears
then we may write 
\be \showlab{dft}
   \tilde{f}^{(n)} = F^{(m)}\vec{f}^{(m)}, 
\ee
where $F^{(m)}$ is the matrix given by $F^{(m)}_{jk} = \omega_m^{-jk}$. 
It is called the ``Fourier matrix'', and (\ref{dft}) is referred to as 
the ``Discrete Fourier Transform'' DFT (of degree $n$) of $f$. 
See \cite[\S 3.8]{A}.

 We have taken the point of view that we're applying the Projection 
Theorem and computing the corresponding inner products via Gaussian quadraure.
The community of mathematicians and computational scientiests who use
this method in differential equations (such as \cite{CHQZ}) 
refer to the Fourier projection approximations as a 
``spectral method'', and the method we've just derived as 
``pseudospectral method'', also as ``collocation''.
However, there is another useful way to view the DFT, and that is as an
{\em interpolation} scheme. Indeed, this is the point of view taken by 
Atkinson in \cite[\S 3.8]{A}. 
\begin{theorem}
  The approximation $f_n$ given by (\ref{fn.def}) satisfies
\[
    f_n(x^{(m)}_j) = f(x^{(m)}_j), \quad j=0,\ldots,m.
\]
\end{theorem}
\pf 
By the uniqueness of the interpolating trig. polynomial, there is 
exactly one $g=\in \trigs_n$ such that 
\[
    g(x^{(m)}_j) = f(x^{(m)}_j), \quad j=0,\ldots,m,
\]
and it suffices that we show 
\[
    f_n(x^{(m)}_j) = g(x^{(m)}_j), \quad j=0,\ldots,m.
\]
Now, suppose we apply our method to this $g$. Since 
$g\in \trigs_n$ it has the form
\[
   g(x) = \sum_{k=-n}^n \hat{g}_k e^{ikx},
\]
so that the aliasing error is zero.
Hence our method recovers $g$ exactly. \QED

%%%%%%%%%%% Figure ``.eps'' %%%%%%%%%%%%%
\begin{figure}
\centerline{\hbox{\psfig{figure=dft1.eps,height=3in,width=6in}}}
\caption{The ``true'' $f$ for which we're the DFT or ``pseudospectral''
 approximation is computed.
}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%% Figure ``.eps'' %%%%%%%%%%%%%
\begin{figure}
\centerline{\hbox{\psfig{figure=AB_tilde.eps,height=6in,width=6in}}}
\caption{The FFT-constructed cosine and sine coefficients.}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\newpage

\noindent\underline{{\tt ``dftdemo.m''}}
\begin{verbatim}
%Purpose: Demonstrate the "dft".
%
% First we set up a trig. polynomial of degree N, and a set of 
% x-values for evaluation for graphing, etc. The parameter "sob"
% refers to a Sobolev smoothness parameter.
  sob = 1/2;
  N = 2^10-1;
  decayrate = (1:N).^(sob+1/2);
  A0 = 10;  A=randn(N,1)./decayrate';    B=randn(N,1)./decayrate';

  M = 2*(N+1);
  k=(0:M-1);
  xk = 2*k*pi/M;

  fk = trigeval(A0,A,B);
  [A0_tilde,A_tilde,B_tilde] = triginterp(fk);

  clg; subplot(111);
  plot(xk,fk,'w'); 
  title(' The "true" f'); grid on;  xlabel('x-->'); ylabel('y-->'); 
  drawnow;
  print -deps dft1.eps

  subplot(211);  
     plot((1:N),A_tilde,'w'); 
     title('Freqency-domain (Cosine coefficients)');
     grid on;
  subplot(212);  
     plot((1:N),B_tilde,'w'); 
     title('Freqency-domain (Sine coefficients)');
     grid on;
  print -deps AB_tilde.eps;
\end{verbatim}

\newpage
\noindent\underline{{\tt ``triginterp.m''}}
\begin{verbatim}
  function [A0_tilde, A_tilde, B_tilde] = triginterp(fvals);
%
% Use the FFT to do a fastr trig interpolation. 
% The length of fvals should be 2^(r+1) for some r.
% The interpolation is done for this many evenly spaced x's
% on the interval [-pi,pi] (right-end not included in the x's).

 M = max(size(fvals)); N=M/2-1;
 ftilde = fft(fvals)/M;
 j=1:M/2;

% Convert the d's back to A's and B's.
 A0_tilde = ftilde(1);
 A_tilde  = real( ftilde(j+1) + ftilde(M-j+1) );
 B_tilde  = imag( ftilde(j+1) - ftilde(M-j+1) );

 
 A_tilde = A_tilde(1:N);
 B_tilde = B_tilde(1:N);
\end{verbatim}

\vfil
\noindent\underline{{\tt ``trigeval.m''}}
\begin{verbatim}
  function fval = trigeval(A0, Avec, Bvec);

% use the FFT to do a fast trig. evaluation. 
% The length of Avec and Bvec should be 2^r-1 for some r.

  Ncoef = max(size(Avec));
  M = 2*(Ncoef+1);
  fhat0 = A0;
  fhat_plus = (Avec + i*Bvec)/2;
  fhat = [ A0;
           fhat_plus;
           0;
           flipud(conj(fhat_plus))
         ];
  
  fval = M*real(ifft(fhat));
\end{verbatim}


\begin{thebibliography}{99}

\bibitem[A]{A}
      Kendall E. Atkinson, {\em An Introduction
      to Numerical Analysis}, $2^{nd}$ ed., Wiley, 1989.

\bibitem[CHQZ]{CHQZ}
  C. Canuto, M.~Y. Hussaini, A. Quarteroni, and T.~A. Zang,
  {\em Spectral Methods in Fluid Dynamics}, 
  Springer series in computational physics, 
  Springer, 1987.

\bibitem[vL]{vL}
  van Loan, Charles, 
  {\sl Computational Frameworks for the Fast Fourier Transform}, 
  SIAM, 1992.

\end{thebibliography}

\end{document}