\documentclass[12pt]{article}
\usepackage{fullpage,amsmath,amsfonts}

\newcommand{\x}{\boldsymbol{x}}
\renewcommand{\v}{\boldsymbol{v}}
\newcommand{\w}{\boldsymbol{w}}
\renewcommand{\a}{\boldsymbol{a}}
\renewcommand{\b}{\boldsymbol{b}}
\renewcommand{\c}{\boldsymbol{c}}
\newcommand{\e}{\mathrm{e}}
\renewcommand{\i}{\mathrm{i}}
\renewcommand{\Re}{\operatorname{Re}}
\renewcommand{\Im}{\operatorname{Im}}
\newcommand{\Range}{\operatorname{Range}}
\newcommand{\Ker}{\operatorname{Ker}}
\newcommand{\Span}{\operatorname{Span}}
\renewcommand{\dim}{\operatorname{dim}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\C}{\mathbb{C}}


\begin{document}

\title{Numerical Methods I}
\author{Lab Session 6}
\date{October 30, 2003}
\maketitle

\noindent
In \emph{trigonometric interpolation} we use a truncated Fourier
series to interpolate a periodic function $f$ on the interval
$[0,2\pi]$.  Let $N$ be an even number of equidistant nodes $x_j=jh$
where $h=2\pi/N$ and $j=0,\dots,N-1$.

The trigonometric interpolant of $f$ is the function
\begin{equation}
  \tilde f(x) = \sum_{k=-N/2}^{N/2-1} c_k \, \e^{\i k x} \,.
  \label{e.reconstruction}
\end{equation}
To determine the coefficients $c_k$, we impose the interpolation
condition $\tilde f(x_j)=f(x_j)$ for $j=0, \dots, N-1$, i.e.
\begin{equation}
  \sum_{k=-N/2}^{N/2-1} c_k \, \e^{\i k x_j} = f(x_j) \,.
\end{equation}
We multiply this equation by $\e^{-\i m x_j}$ and sum over all $j$, so
that
\begin{equation}
  \sum_{j=0}^{N-1} \e^{-\i m x_j}
  \sum_{k=-N/2}^{N/2-1} c_k \, \e^{\i k x_j} =
  \sum_{j=0}^{N-1} \e^{-\i m x_j} \, f(x_j) \,.
\end{equation}
Notice that
\begin{align}
  \sum_{j=0}^{N-1} \e^{-\i m x_j}
  \sum_{k=-N/2}^{N/2-1} c_k \, \e^{\i k x_j}
  & = \sum_{k=-N/2}^{N/2-1} c_k \sum_{j=0}^{N-1} \e^{\i (k-m) x_j} \notag \\
  & = \sum_{k=-N/2}^{N/2-1} c_k \sum_{j=0}^{N-1} \e^{\i (k-m) h j} \notag \\
  & = \sum_{k=-N/2}^{N/2-1} c_k \cdot
  \begin{cases}
    N & \text{if } k=m \\
    \dfrac{1 - (\e^{\i (k-m) 2\pi/n})^n}{1 - \e^{\i (k-m) 2\pi/n}} = 0
    & \text{if } k \neq m 
  \end{cases} \notag \\
  & = N \, c_m \,.
\end{align}
We therefore obtain that
\begin{equation}
  c_m = \frac1N \sum_{j=0}^{N-1} \e^{-\i m x_j} \, f(x_j) \,.
  \label{e.dft}
\end{equation}

\paragraph{Remark:} The normalization convention is arbitrary.  Most
FFT libraries, including \textsf{Octave}, have the factor $1/N$ in the
inverse transform \eqref{e.reconstruction} rather than in the forward
transform \eqref{e.dft}.  The convention above has the advantage that
the \emph{discrete Fourier transform} \eqref{e.dft} can be seen as the
Riemann sum approximation of the continuous Fourier transform.  If
both the forward and the inverse transform get a factor of $1/\sqrt
N$, the transform is unitary, which has certain theoretical
advantages, but is tedious in many instances.

\begin{enumerate}

\item Write an \textsf{Octave} program to compute the coefficients
$c_m$ for $m=-\tfrac{N}2,\dots,\tfrac{N}2-1$ for a given function $f$,
e.g.
\[
   f(x) = x \, (x-2 \, \pi) \, \e^{-x} \,.
\]

\item This operation can be done very efficiently via the \emph{Fast
Fourier Transform}, available in \textsl{Octave} as the function
\texttt{fft}.  Compare your result with that produced by \texttt{fft}
to find out how the coefficients $c_m$ are laid out.

\item Write an \textsf{Octave} function \texttt{tpolyval(c,x)}, which,
in analogy with the built-in function \texttt{polyval}, evaluates the
trigonometric interpolant $\tilde f$ via equation
\eqref{e.reconstruction}.  Use the coefficient layout as used by
\texttt{fft}.

\item Plot $f$ and $\tilde f$ for the given example.

\item What happens when $N=10$ and you interpolate $f(x)=\sin(6x)$?

\end{enumerate}
\end{document}
