\documentclass{slides}
\usepackage{epsfig}
\newcommand{\maintitle}[1]{\begin{center}\Large #1 \end{center}}
\newcommand{\slidetitle}[1]{\begin{par}\Large #1\end{par}}
\newtheorem{theorem}{Theorem}

\begin{document}

\maintitle{An Introduction to Wavelets -- The Engineers' Perspective}


\vfil
\begin{center}
{\large Faculty Expository Seminar, 10/17 and 10/24}

\vfil
{
John B. Little
}

\vfil

Department of Mathematics and Computer Science

{\small
College of the Holy Cross

Worcester, MA 01610

{\tt email: little@mathcs.holycross.edu}}
\end{center}



%%%%%%%%%%%%%%%%
\begin{slide}
\slidetitle{Abstract}

\begin{itemize}
\item In these talks, we will present
a rapid introduction to {\it wavelets} -- a
recently developed supplement or alternative 
to Fourier methods in signal and imaging processing,
etc.  

\item We will focus on the ``Fast Wavelet Transform'' 
approach, filter banks, etc. -- the engineers'
perspective on this subject rather than the analysts' 
perspective, though of course there will be a significant
analysis component too.
\end{itemize}
\end{slide}
\begin{slide}
\slidetitle{{Outline of Talks}}

\begin{enumerate}
\item Wavelets
\item Example -- Haar Wavelets
\item Filter Banks
\item Example -- Daubechies Wavelets
\item Some Generalizations -- Multiwavelets
\end{enumerate}
\end{slide}


\begin{slide}
\slidetitle{Wavelets}

In recent years, wavelets have become an attractive alternative
or complement to Fourier-based methods for applications such as

\begin{itemize}
\item signal analysis
\item denoising
\item data compression (e.g. FBI fingerprint record system)
\item statistics
\item finite element methods for PDE
\end{itemize}
\end{slide}

\begin{slide}
\slidetitle{Why Wavelets?}

Fourier expansions express signals as 
sums of sinusoidal 
components of different frequencies.

But sinusoids $\sin\left({n\pi t\over L}\right)$, 
$\cos\left({n\pi t\over L}\right)$ are not {\it time-limited}
the way real-world signals are.

Wavelet expansions express signals as sums of {\it time/distance
localized} components at different time/distance
scales (``detail levels'') rather than as sums of 
components of different frequencies.

This often yields better, more compact approximations 
for strictly time-limited signals.
\end{slide}

\begin{slide}
\slidetitle{Multi-Resolution Analysis}

One framework for this approach:

Consider {\it continuous-time} signals in $L^2({\bf R})$.
A {\it multiresolution analysis} of $L^2$ is a sequence
of closed subspaces $V_j$, $j \in {\bf Z}$ satisfying:

\begin{enumerate}
\item $V_j \subset V_{j+1}$ for all $j$
\item $f(t) \in V_j \Leftrightarrow f(2t) \in V_{j+1}$
\item $\cap_{j\in {\bf Z}} V_j = \{0\}$
\item $\overline{\cup_{j\in {\bf Z}} V_j} = L^2({\bf R})$
\item There exists $\phi \in V_0$ such that 
$\{\phi(t-k) : k \in {\bf Z}\}$ is an orthonormal basis
of $V_0$.
\end{enumerate}

$\phi$ is called the {\it scaling function}.
\end{slide}

\begin{slide}
\slidetitle{MRA to Wavelets}

Let $W_j$ be the orthogonal complement of $V_j$ in 
$V_{j+1}$, so $V_{j+1} = V_j \oplus W_j$.  Then
$$L^2({\bf R}) = \overline{V_0 \oplus W_0 \oplus W_1 \oplus W_2 \oplus \cdots}$$

We assume that
there is some $\psi \in W_0$, called the 
{\it mother wavelet}, such that 
$\{\psi(t-k) : k \in {\bf Z}\}$ is an orthonormal basis
of $W_0$.  
From 2 in definition of MRA, $W_j$ has an orthonormal basis
$\{2^{j/2}\psi(2^j t - k) : k \in {\bf Z}\}$
formed from scaled, dilated, and translated versions of $\psi$.

General wavelet expansion of a signal $f$ looks like this:
$$f(t) = \sum_k a_{0k} \phi(t-k) + \sum_{j=0}^\infty \sum_k b_{jk}2^{j/2}\psi(2^j t - k)$$
\end{slide}

\begin{slide}
\slidetitle{Wavelet Coefficients}

If appropriate orthogonality properties hold, the 
wavelet expansion coefficients of a signal
can be computed in 
a fashion similar to Fourier coefficients:

\begin{eqnarray*}
a_{0k} &=& \int_{\bf R} f(t) \overline{\phi(t - k)}\ dt\\
&&\\
b_{jk} &=& \int_{\bf R} f(t) \overline{2^{j/2} \psi(2^jt - k)}\ dt\\
\end{eqnarray*}

The larger $j$ is, the higher the level of detail encoded
by the $b_{jk}$.
\end{slide}

\begin{slide}
\slidetitle{Dilation and Wavelet Equations}

Since $V_0$ and $W_0$ are contained in $V_1$,
and 2 in the definition of MRA implies
that $\{\phi(2 t - k) : k \in {\bf Z}\}$ is a 
basis for $V_1$, we have
\begin{itemize}
\item the {\it dilation equation}:
$$\phi(t) = \sqrt{2} \sum_{k\in {\bf Z}} h_0[k] \phi(2t - k)$$
for some coefficients $h_0[k]$, and 
\item the {\it wavelet equation}:
$$\psi(t) = \sqrt{2} \sum_{k\in {\bf Z}} h_1[k] \phi(2t - k)$$
for some coefficients $h_1[k]$.
\end{itemize}
\end{slide}

\begin{slide}
\slidetitle{Haar Wavelets}

The simplest example of this setup is the Haar 
wavelet system.  

The Haar scaling function is the ``box'' function
$$\phi(t) = \cases{1 & if $0 \le t < 1$\cr
                   0 & otherwise\cr}$$
The dilation equation is 
$$\phi(t) = \phi(2t) + \phi(2t - 1),$$
that is, $h_0[0] = h_0[1] = {1\over\sqrt{2}}$ and all others $0$.

The Haar wavelet is 
$$\psi(t) = \cases{1 & if $0 \le t < 1/2$\cr
                   -1& if $1/2 \le t < 1$\cr
                   0 & otherwise\cr}$$
The wavelet equation is
$$\psi(t) = \phi(2t) - \phi(2t - 1)$$
so $h_1[0] = {1\over\sqrt{2}}, h_1[1] = -{1\over\sqrt{2}}$, all others 0.  
\end{slide}

% Haar Wavelets Maple Session slides here

\begin{slide}
\slidetitle{Analysis and Synthesis Equations}

Back to general setup.

Let $f \in V_{j+1} = V_j \oplus W_j$.  Then 
$$f(t) = \sum_k a_{j+1,k} 2^{(j+1)/2} \phi(2^{j+1}t - k),\ {\rm or}$$
$$f(t) = \sum_k a_{j,k} 2^{j/2} \phi(2^j t - k) + \sum_k b_{j,k} 2^{j/2} \psi(2^j t - k)$$
From the dilation and wavelet equations, get
\begin{itemize}  
\item the {\it analysis equations}
\begin{eqnarray*}
a_{j,k} &=& \sum_n h_0[n-2k] a_{j+1,n}\\
b_{j,k} &=& \sum_n h_1[n-2k] a_{j+1,n}\\
\end{eqnarray*}
\item and the {\it synthesis equation}
$$a_{j+1,n} = \sum_k (h_0[n-2k] a_{j,k} + h_1[n-2k] b_{j,k})$$ 
\end{itemize}
\end{slide}

\begin{slide}
\slidetitle{Wavelets and Filter Banks}

Key observation here is that these have {\it same form} as the outputs
of 2-channel {\it filter banks} used for multirate digital
signal processing tasks such as
sampling rate conversion, subband processing, etc.

The $h_0[n]$ and $h_1[n]$ can be understood as 
coefficients of two digital filters acting on ``signals''
composed of the wavelet expansion coefficients $a_{j,k}$, $b_{j,k}$
at detail level $j$, yielding the $a_{j+1,n}$ via the synthesis
relation.

Similarly, two related filters in the analysis
relation produce
$a_{j,k}$ and $b_{j,k}$ from $a_{j+1,n}$.


Let's take a short ``detour'' to ``recall''
some DSP fundamentals.
\end{slide}

\begin{slide}
\slidetitle{Digital Signals}

A {\em digital} signal is a quantized function of a 
{\em discrete} variable.  
\begin{itemize}
\item only have signal values at 
$n\Delta t$, $n$ an integer, $\Delta t$ the 
{\em sampling interval} 
\item signal values are elements of a discrete set.
\end{itemize}

For simplicity, we ignore quantization effects
and study {\em discrete-time, complex-valued} signals.

Can write a signal as a complex-valued function of 
an integer variable, $x[n]$, or equivalently, a 
list of complex numbers indexed by the integers:
\[
 \cdots, x[-2], x[-1], x[0], x[1], x[2], \cdots
\]
\end{slide}


\begin{slide}
\slidetitle{Filters}

Signal processing operations can be described 
mathematically by means of functions $\Gamma$:
\[
\Gamma : \{ {\rm signals} \} \to \{ {\rm signals} \}.
\]

We mainly consider operations with the following properties

\begin{enumerate}
\item (Linearity)  If $z[n] = c x[n] + y[n]$, then
\[
\Gamma(z)[n] = c \Gamma(x)[n] + \Gamma(y)[n]
\]
\item (Shift-invariance)  If $y[n] = x[n + k]$ for all 
$n$, $k$ fixed, and the outputs are
$w = \Gamma(x)$, $z = \Gamma(y)$, then
\[
z[n] = \Gamma(y)[n] = \Gamma(x)[n+k] = w[n+k]
\]
\end{enumerate}

A $\Gamma$ satisfying both of these properties is said
to be a (linear, shift-invariant) {\em filter}.
\end{slide}

\begin{slide}
\slidetitle{A ``Moving Average'' Filter}

For any input signal $x[n]$, define 
\begin{eqnarray*}
y[n] &=& \Gamma(x)[n]\\
     &=& \left(\sum_{k=0}^4 x[n-k]\right)/5
\end{eqnarray*}

If $x[n] = \cos(n\pi/2)$, 
\[
\cdots,1,0,-1,0,1,0,-1,0,1,\cdots
\]
for instance, then $\Gamma(x)[n] = $
\[
\cdots, 0,1/5,0,-1/5,0,1/5,0,-1/5,0,\cdots
\]

It is easy to see that this satisfies the linearity
and shift-invariance properties, so this defines
an filter.
\end{slide}
\begin{slide}
\slidetitle{Impulse Response of a Filter}

A nice result about filters:  Let $\delta$ be 
\[
\cdots, 0,0,1,0,0, \cdots
\]
($1$ at time $n = 0$) -- called a {\em unit
impulse} input.  $\Gamma(\delta)$ is called the 
{\em impulse response} of the filter.

{\em If we know $h = \Gamma(\delta)$, then we
also know the output of the filter on any other
signal $x$.}  Reason: We can write
\[
x[n] = \sum_{k=-\infty}^\infty x[k]\delta[n-k]
\]
Hence by linearity and shift-invariance
\[
y[n] = \Gamma(x)[n] = \sum_{k=-\infty}^\infty x[k] h[n-k].
\]
\end{slide}

\begin{slide}
\slidetitle{The Transfer Function}

Package $x[n]$ by its ``$z$-transform'', the series
\[
X(z) = \sum_{n=-\infty}^\infty x[n] z^{-n}
\]
Do same for the impulse response $h[n]$
and the output $y[n]$ for input $x[n]$:
\[
H(z) = \sum_{n=-\infty}^\infty h[n] z^{-n},\ 
Y(z) = \sum_{n=-\infty}^\infty y[n] z^{-n} 
\]

Coefficient of $z^{-n}$ in the product 
$H(z)X(z)$ is the discrete convolution:
\[
\sum_{k=-\infty}^\infty h[k]x[n-k],
\]
same as $y[n]$ (!)  

$\Rightarrow Y(z) = H(z)X(z)$.  $H(z)$ is called the 
{\em transfer function}.
\end{slide}

\begin{slide}
\slidetitle{Connection with Fourier Domain}

Note:  
$$H(e^{i\omega}) = \sum_{n=-\infty}^\infty h[n] e^{-i n\omega}$$
is the (discrete-time) Fourier transform 
of $h$.

So $H(z)$ also determines the 
{\it frequency response} characteristics of the filter on 
input signals.

For instance, with input a pure sinusoid $x[n] = e^{i n\omega_0}$ the output 
$y[n] = \Gamma(x)[n]$
\begin{eqnarray*}
     &=& \sum_{k=-\infty}^\infty h[k] x[n-k]\\
     &=& \sum_{k=-\infty}^\infty h[k] e^{i (n-k)\omega_0}\\
     &=& \left(\sum_{k=-\infty}^\infty h[k] e^{-ik\omega_0}\right) e^{in\omega_0}\\
     &=& H(e^{i\omega_0}) x[n]\\
\end{eqnarray*}    
\end{slide}
\begin{slide}
\slidetitle{Frequency Response}

The magnitude of the frequency response for the 
5-point moving average filter from before:

\psfig{figure=5Pt.eps,height=3in,width=6in}

This is an example of a {\it ``low-pass''} filter -- it
``passes'' sinusoids of small normalized frequency  
and ``stops'' sinusoids of large normalized frequency.
\end{slide}

\begin{slide}
\slidetitle{Back to Filter Banks}

The $h_0[n]$ and $h_1[n]$ are coefficients of
two filters appearing in the {\it synthesis filter
bank}:

Synthesis bank:

\begin{picture}(350,200)
\put(40,137){$a_{j,k}$}
\put(40,57){$b_{j,k}$}
\put(80,140){\vector(1,0){30}}
\put(80,60){\vector(1,0){30}}
\put(110,133){\fbox{$(\uparrow 2)$}}
\put(110,53){\fbox{$(\uparrow 2)$}}
\put(170,140){\vector(1,0){40}}
\put(170,60){\vector(1,0){40}}
\put(210,134){\fbox{$H_0$}}
\put(210,54){\fbox{$H_1$}}
\put(250,140){\line(1,0){40}}
\put(250,60){\line(1,0){40}}
\put(290,140){\vector(0,-1){30}}
\put(290,60){\vector(0,1){30}}
\put(290,100){\circle{20}}
\put(290,110){\line(0,-1){20}}
\put(280,100){\line(1,0){20}}
\put(300,100){\vector(1,0){40}}
\put(350,97){$a_{j+1,k}$}
\end{picture}

where $(\uparrow 2)$ is {\it upsampling}, or {\it ``interpolation''} by 2
(expanding with zeroes).

$H_0$ is the {\it synthesis low-pass filter}; $H_1$ is 
the {\it synthesis high-pass filter}.  For example, in Haar case,
$H_0$ is two-point averaging filter, $H_1$ is two-point difference filter.

\end{slide}

\begin{slide}
\slidetitle{Filter Banks, cont.}

Similarly, the analysis relations correspond to 
an analysis filter bank:

\begin{picture}(400,200)
\put(40,97){$a_{j+1,n}$}
\put(110,100){\vector(1,0){30}}
\put(140,100){\line(0,1){40}}
\put(140,100){\line(0,-1){40}}
\put(140,140){\vector(1,0){30}}
\put(140,60){\vector(1,0){30}}
\put(170,130){\fbox{$\widetilde{H_0}$}}
\put(170,50){\fbox{$\widetilde{H_1}$}}
\put(210,140){\vector(1,0){40}}
\put(210,60){\vector(1,0){40}}
\put(250,132){\fbox{$(\downarrow 2)$}}
\put(250,52){\fbox{$(\downarrow 2)$}}
\put(310,140){\vector(1,0){40}}
\put(310,60){\vector(1,0){40}}
\put(355,137){$a_{j,k}$}
\put(355,57){$b_{j,k}$}
\end{picture}

(where $\widetilde{H_i}(z) = H_i(1/z)$ and $(\downarrow 2)$ 
is {\it downsampling}, or {\it decimation} by 2)

\end{slide}

\begin{slide}
\slidetitle{``Fast Wavelet Transform''}

The filter bank approach to wavelet systems is 
fruitful for several reasons:

\begin{itemize}
\item {\it Iteration of the low-pass analysis branch}, starting
from $a_{j+1,k}$ for some large $j+1$ (or sampled signal values) gives 
a very fast way to compute wavelet expansions:

\begin{picture}(400,150)
\put(-40,47){$V_{j+1}$}
\put(10,50){\vector(1,0){15}}
\put(25,50){\line(0,1){25}}
\put(25,50){\line(0,-1){25}}
\put(25,25){\vector(1,0){15}}
\put(25,75){\vector(1,0){15}}
\put(40,15){\fbox{$\widetilde{H_0}$}}
\put(40,65){\fbox{$\widetilde{H_1}$}}
\put(80,25){\vector(1,0){15}}
\put(80,75){\vector(1,0){15}}
\put(95,17){\fbox{$(\downarrow 2)$}}
\put(95,67){\fbox{$(\downarrow 2)$}}
\put(155,25){\vector(1,0){15}}
\put(155,75){\vector(1,0){15}}
\put(170,17){$W_j$}
\put(170,67){$V_j$}
\put(195,75){\vector(1,0){15}}
\put(210,75){\line(0,1){25}}
\put(210,75){\line(0,-1){25}}
\put(210,100){\vector(1,0){15}}
\put(210,50){\vector(1,0){15}}
\put(225,90){\fbox{$\widetilde{H_0}$}}
\put(225,40){\fbox{$\widetilde{H_1}$}}
\put(265,100){\vector(1,0){15}}
\put(265,50){\vector(1,0){15}}
\put(280,93){\fbox{$(\downarrow 2)$}}
\put(280,43){\fbox{$(\downarrow 2)$}}
\put(340,100){\vector(1,0){15}}
\put(340,50){\vector(1,0){15}}
\put(355,92){$V_{j-1} \cdots$}
\put(355,42){$W_{j-1}$}
\end{picture}




\end{itemize}
\end{slide}

\begin{slide}
\slidetitle{Wavelet System Design}

Can {\it design wavelet systems} by designing the analysis 
(or synthesis) bank filters.

\begin{itemize}
\item Often want $H_i$ to be FIR filters $\Rightarrow$
$\phi, \psi$ have {\it compact support}.
\item Often want $\{\phi(t-k), \psi(2^jt - k)\}$
{\it orthogonal} in $L^2({\bf R})$ -- 
Sufficient: {\it ``double-shift orthogonality''}:
\begin{equation} 
\sum_{n \in {\bf Z}} h_i[n] h_j[n - 2k] = \delta[i - j]\delta[k] 
\end{equation} 
(condition also implies
``perfect reconstruction'' -- synthesis bank inverts analysis bank).
\item Often want good approximation properties -- for approximation
to order $p$, sufficient: $H_0(z)$ has {\it zero order $p \ge 1$ at} $z = -1$ ($\omega = \pi$).
(Condition also related to {\it smoothness} of $\phi,\psi$.)
\end{itemize}
\end{slide}

\begin{slide}
\slidetitle{The ``Cascade Algorithm''}

Recall the {\it dilation equation}:
$$\phi(t) = \sqrt{2} \sum_{k\in {\bf Z}} h_0[k] \phi(2t - k)$$
where the coefficients $h_0[k]$ in the synthesis
low pass filter appear.  

The dilation equation says that $\phi(t)$ is a fixed
point of a certain operator $T$ on $L^2$.  

Hence if we start from the coefficients in a filter $H_0$,
we can try to find $\phi$ by a {\it fixed point iteration}:
$$\phi^{(i+1)} = T(\phi^{(i)})$$

This is called the ``cascade algorithm'' for determining
the scaling function of a wavelet system.
\end{slide}
\begin{slide}
\slidetitle{Convergence of the Cascade Algorithm}

For {\it convergence} of the cascade algorithm, need an
eigenvalue condition:  $T$ has $\lambda = 1$ as a simple
eigenvalue, and all other $|\lambda| < 1$.  

$H_0(-1) = 0$ is a {\it necessary} condition for this to 
hold.
\end{slide}

%Daubechies Example Here
\begin{slide}
\slidetitle{To Multiwavelets}
\begin{itemize}
\item There are limitations on what can be accomplished
this way, though.
\item For instance, a symmetric
orthogonal FIR filter can have only two nonzero coefficients.
\item To get around these limitations, people have
introduced several generalizations of MRA, etc.
\item We'll consider one -- {\it multiwavelets}.
\end{itemize}
\end{slide}

\begin{slide}
\slidetitle{Multiwavelet Systems}

Basic idea here is to allow {\it more than one} scaling function,
mother wavelet:

\begin{itemize}
\item $V_j = {\rm Span} \{\phi_0(2^j t - k), \ldots, 
\phi_{r-1}(2^j t-k) : k \in {\bf Z}\}$
\item $W_j = {\rm Span} \{\psi_0(2^j t - k), \ldots, 
\psi_{r-1}(2^j t-k) : k \in {\bf Z}\}$
\item {\it dilation equation} becomes
$$\Phi(t) = \sqrt{2} \sum_{n\in {\bf Z}} C_0[n] \Phi(2t - n)$$
where $\Phi = (\phi_0,\ldots,\phi_r)^T$, and $C_0[n]$ is an $r\times r$
matrix for each $n$.  Similarly for {\it wavelet equation}.
\end{itemize}
\end{slide}

\begin{slide}
\slidetitle{Multiwavelets for Scalar Signals}

Get multiwavelet analysis and synthesis filter banks 
similar to the ordinary case.

Although one natural way to use multiwavelets is for
{\it vector-valued signals}, we will consider a different
use.  

For simplicity, take $r = 2$, and write
$$C_0[n] = \pmatrix{h_0[2n]&h_0[2n+1]\cr
                    h_1[2n]&h_1[2n+1]\cr},$$
then the entries in the $C_0[n]$ are split into coefficients
from {\it two scalar filters}, $h_i$, $i = 0,1$.  Similarly for
matrix filter $C_1[n]$ from the wavelet equation.

Group consecutive sampled values from a scalar signal into 
2-component vectors.
\end{slide}

\begin{slide}
\slidetitle{Multiwavelets for Scalar Signals, cont.}

\begin{itemize}
\item The 2-channel multiwavelet filter banks can be
regrouped into {\it 4-channel scalar filter banks}
(with interleaving of outputs).
\item Condition for {\it orthogonality} (analogous to (2) above)
becomes
\begin{equation} 
\sum_{n \in {\bf Z}} h_i[n] h_j[n - 4k] = \delta[i - j]\delta[k] 
\end{equation}
\item Natural analog of the ``zeroes at $z = -1$ ($\omega = \pi$)''
condition for approximation order is the {\it ``balanced'' property}: 
Selesnick, Lebrun and Vetterli, etc. 
\end{itemize}
\end{slide}

\begin{slide}
\slidetitle{Balanced Multiwavelet Systems}

{\em references:}

1. Selesnick, Multiwavelet Bases with Extra
Approximation Properties, 

2. Lebrun and Vetterli, Balanced Multiwavelets Theory and Design,

Both: {\sl IEEE Trans. Signal Proc.}, v. 46, 1998.

{\bf Def.} The $r = 2$ multiwavelet system with filters $h_0,h_1$ as above
is {\it balanced} of order $p$ if $(1 + z^{-1} + z^{-2} + z^{-3})^p$ divides
$H_0(z) + \alpha_p(z^4)H_1(z)$, where $\alpha_p$ is given by an explicit
formula yielding:
\begin{eqnarray*}
\alpha_1 &=& 1\\
\alpha_2 &=& \frac{3-z^{-1}}{2}\\
\alpha_3 &=& \frac{15-10z^{-1}+3z^{-2}}{8}, etc.\\
\end{eqnarray*}
\end{slide}

\begin{slide}
\slidetitle{Balanced Multiwavelets, cont.}

An orthogonal multiwavelet system which is balanced of order
$p$ has approximation order (at least) $p$:

The polynomial signals $1,t,t^2,\ldots,t^{p-1}$
can be {\it exactly} decomposed (on a finite interval) 
into linear combinations of the scaling functions
$\phi_i(t)$ and their integer translates.

The balancing order is also related to the order
of smoothness of the $\phi_i(t)$.

\end{slide}


\end{document}
