\documentclass[10pt,letterpaper]{article}
\usepackage{graphicx,amssymb}
\usepackage[colorlinks=true]{hyperref}

\topmargin -0.5in
\textheight 9in
\oddsidemargin=-0.05in
\evensidemargin=-0.05in
\textwidth 6.5in
\pagestyle{empty}

\newcommand{\bi}{\begin{itemize}}
\newcommand{\ei}{\end{itemize}}
\newcommand\ben{\begin{enumerate}}
\newcommand\een{\end{enumerate}}
\newcommand\vg{\vspace{2ex}}
\newcommand{\bc}{\begin{center}}
\newcommand{\ec}{\end{center}}
\newcommand{\ie}{{\it i.e.\ }}
\newcommand{\eg}{{\it e.g.\ }}
\newcommand{\eps}{\varepsilon}
\newcommand{\mbf}[1]{{\mathbf #1}}
\newcommand{\tbox}[1]{{\mbox{\tiny #1}}}


\begin{document}
\title{Math 56 Compu \& Expt Math, Spring 2013: Homework 5}
\date{\vspace{-2ex} due 10am Thursday May 2nd}
\maketitle

{\em This week less theory, but more fun real-world signal applications}


\ben

\item Spectral differentiation, i.e.\ approximating the derivative of
a function from $N$ samples of its value.%
\footnote{You should compare
this to finite-differencing, which is the other numerical
differentiation method you have seen.}
\ben
\item Adapt your code from \#3(c) from the previous homework to
compute the {\em derivative} of the interpolant on the fine grid,
generated from the same sets of $N$ samples of $e^{\sin x}$.
Plot the convergence vs $N$ of the max error from the analytically-known
derivative of $f$ on the fine grid.
[Hint: what are the Fourier coefficients of the derivative?]
\item Repeat for $|\sin x|^3$. What type and order of convergence results?
Is it as predicted by \#1(f) from HW4?
\een

\item Write your own FFT code for $N=2^n$, using the Cooley-Tukey
recursive algorithm, 
including a driver code that computes the 2-norm of the error between
the output of your FFT and  Matlab's native {\tt fft}.
Measure the ratio of your runtime to Matlab's on a vector of length $2^{16}$.
(Again, marvel at the speed of the FFTW library that Matlab uses.)

\item The power of the DFT to extract a periodic signal
in the presence of noise.
\ben
\item
Download and read in the audio file {\tt signoise.wav}.
This signal $\mbf{f}$ is known to contain one frequency component on top of
(iid, i.e\ uncorrelated) noise.
Can you see any periodicity in the graph of $\mbf{f}$ visually?
Use the FFT to answer: what are the {\em two}
dominant Fourier mode indices $m$ ?
% can't be your own, since not length 2^n
How are their indices related? (Explain.)
How are their coefficients related?
Explain using a property of the input signal.
%\item The desired component is in the form $A \sin(2\pi q t + \phi)$
%where $t$ is time. What is the amplitude $A$ ?
\item Use the information that the signal was sampled at 44100 Hz
to compute the true (physical) dominant frequency in Hz.
\item [BONUS]
Listen to the signal; can you hear any tone above the noise?
Discuss by comparing the amplitude of the desired component
to the typical size of the signal $\mbf{f}$.
\een


\item Discrete convolution theory.
\ben
\item
Consider the signal $f = [1,1,1,0,0,\ldots]^T$,
where the dots indicate continuation by zeros to the length of the vector,
and the signal $g = [1,1,1,1,0,\ldots]^T$.
Assuming that the length of both the signals is $N=6$, work out
by hand the (periodic) convolution $f\ast g$.
\item
Do the same except assuming that the length of both is $N=5$.
\item
The answer to one of (a) or (b) has the same nonzero vector elements as
the {\em non-periodic} convolution $\sum_i f_i g_{j-i}$.
Which? Use this to answer: to what length $N$
signals of lengths $N_1$ and $N_2$ should be zero-padded so
that the $N$-periodic convolution correctly computes the non-periodic one.
\item Prove commutativity: for all $f,g \in\mathbb{C}^N$, it holds that
$f\ast g = g \ast f$.
\item
For any $f,g\in\mathbb{C}^N$, what
is the sum of the signal $f\ast g$ expressed in terms of
simple properties of $f$ and $g$? [Hint: use (a) to make a conjecture
then prove it]
\een


\item Using convolution to apply a ``reverb effect'' to an audio signal
(i.e.\ digital signal processing).

Use {\tt audacity}, or any audio software of your choice,
to record (or take from a music track)
an interesting sound of duration at most a few seconds.
Export it as a WAV file (mono, 44100 Hz sample rate)
and read it into Matlab via {\tt wavread}.

Download and read in {\tt impulseresponse.wav}, which is
the sound of a clap in an echoey space, i.e.\ the response of the
space to an impulse (delta function).

Use {\tt fft} and {\tt ifft} to (periodically) convolve the two signals.
Make sure to zero-pad the vectors to a length $N$ long enough so
that the result is correct as a non-periodic convolution.

Finally, use e.g.\ {\tt wavwrite(y,44100,'output.wav');} to write out the
result as a
WAV file, and include in your HW. It should be playable by any audio
player---listen and enjoy!


\item 2D image de-blurring application, when aperture function known.
\ben
\item
Download and {\tt textread} into Matlab the data {\tt blurry.txt}
and reshape it into a 512-by-512 array.
View in monochrome with {\tt imagesc} and {\tt colormap(gray(256)); axis equal},
and note that the text is unreadable!
Similarly read in {\tt aperture.txt} that was the aperture image used
for blurring; it is a single block of nonzero pixels (but appears in four
corners due to periodic wrapping).
\item
Deconvolve the blurry image using the 2D version of the
technique from lecture.
[Hint: the 2D FFT is just like the 1D FFT except it takes a $N\times N$
array to another $N\times N$ array. See {\tt fft2} and {\tt ifft2}.]
Include the resulting (crisp!) image in your HW.
\item Add iid Gaussian white noise of standard deviation $3\times 10^{-5}$
to the original blurry image---this models a {\em tiny}
amount of measurement error (e.g.\ from a camera)---then
repeat (b). Explain why the result looks as it does.%
\footnote{This extreme sensitivity to noise is why I didn't have you
read in the images in TIFF format, even though it is easy in Matlab:
even the rounding error to 8-bit accuracy in the TIFF format introduces
way too much noise!}
[Hint: look at
the DFT of the aperture image.]
\item [BONUS] Invent an improved algorithm that reduces noise in the
deconvolved image while still preserving resolution, and show your
result. [Hint: see above hint!]
\een

\een
\end{document}

% can return to predicting Taylor convergence rates for atan
% when do pi digits.the Fourier seires



% Newton iter for f(x) = 1/x - y which computes 1/y using only mult & add
% (use for arb-prec FFT - based mult)

%Read on SVD, ? - no, skip this

%For {\tt A1} what perturbation of b results in the least
%change in the answer x? This is hard, as is 

