An Introduction to Wavelets –

The Engineers’ Perspective

Faculty Expository Seminar, 10/17 and 10/24

John B. Little

Department of Mathematics and Computer

Science

College of the Holy Cross Worcester, MA 01610

email:  little@mathcs.holycross.edu

Abstract

• In these talks, we will present a rapid introduction to wavelets – a recently developed supplement or alternative to Fourier methods in signal and imaging processing, etc.

• 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.

Outline of Talks

1. Wavelets

2. Example – Haar Wavelets

3. Filter Banks

4. Example – Daubechies Wavelets

5. Some Generalizations – Multiwavelets

Wavelets

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

• signal analysis

• denoising

• data compression (e.g. FBI fingerprint record system)

• statistics

• finite element methods for PDE

Why Wavelets?

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

But sinusoids sin (\frac{nπt}{L}), cos (\frac{nπt}{L}) are not timelimited the way real-world signals are.

Wavelet expansions express signals as sums of 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.

Multi-Resolution Analysis

One framework for this approach: Consider continuous-time signals in L^{2}(R). A multiresolution analysis of L^{2} is a sequence of closed subspaces V_{j}, j ∈ Z satisfying:

1. V_{j} ⊂ V_{j+1} for all j

2. f(t) ∈ V_{j} ⇔ f(2t) ∈ V_{j+1}

3. ∩_{j∈Z}V_{j} = {0}

4. ∪_{j∈Z}V_{j} = L^{2}(R)

5. There exists ϕ ∈ V_{0} such that {ϕ(t − k) : k ∈ Z} is an orthonormal basis of V_{0}.

ϕ is called the scaling function.

MRA to Wavelets

Let W_{j} be the orthogonal complement of V_{j} in V_{j+1} , so V_{j+1} = V_{j} ⊕ W_{j}. Then

L^{2}(R) = V_{0} ⊕ W_{0} ⊕ W_{1} ⊕ W_{2} ⊕ ⋯

We assume that there is some ψ ∈ W_{0}, called the mother wavelet, such that {ψ(t−k) : k ∈ Z} is an orthonormal basis of W_{0}. From 2 in definition of MRA, W_{j} has an orthonormal basis {2^{j/2}ψ(2^{j}t − k) : k ∈ Z} formed from scaled, dilated, and translated versions of ψ.

General wavelet expansion of a signal f looks like this:

f(t) = _{k}^{∑} a_{0k}ϕ(t − k) + ∞ ; ∑ ; j=0 _{k}^{∑} b_{jk}2^{j/2}ψ(2^{j}t − k)

Wavelet Coefficients

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

a_{0k} = ∫_{R} f(t)ϕ(t − k) dt
b_{jk} = ∫_{R} f(t)2^{j/2}ψ(2^{j}t − k) dt

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

Dilation and Wavelet Equations

Since V_{0} and W_{0} are contained in V_{1}, and 2 in the definition of MRA implies that {ϕ(2t − k) : k ∈ Z} is a basis for V_{1}, we have

• the dilation equation:

ϕ(t) = √2 _{k∈Z}^{∑} h_{0}[k]ϕ(2t − k)

for some coefficients h_{0}[k], and

• the wavelet equation:

ψ(t) = √2 _{k∈Z}^{∑} h_{1}[k]ϕ(2t − k)

for some coefficients h_{1}[k].

Haar Wavelets

The simplest example of this setup is the Haar wavelet system. The Haar scaling function is the “box” function

ϕ(t) = \{ \begin{matrix}1 & if 0 ≤ t < 1 \\ 0 & otherwise \end{matrix}

The dilation equation is

ϕ(t) = ϕ(2t) + ϕ(2t − 1),

that is, h_{0}[0] = h_{0}[1] = \frac{1}{√2} and all others 0. The Haar wavelet is

ψ(t) = \{ \begin{matrix}1 & if 0 ≤ t < 1/2 \\ −1 & if 1/2 ≤ t < 1 \\ 0 & otherwise \end{matrix}

The wavelet equation is

ψ(t) = ϕ(2t) − ϕ(2t − 1)

so h_{1}[0] = \frac{1}{√2}, h_{1}[1] = −\frac{1}{√2}, all others 0.

Analysis and Synthesis Equations

Back to general setup. Let f ∈ V_{j+1} = V_{j} ⊕ W_{j}. Then

f(t) = _{k}^{∑} a_{j+1,k} 2^{(j+1)/2} ϕ(2^{j+1} t − k), or

f(t) = _{k}^{∑} a_{j,k}2^{j/2}ϕ(2^{j}t−k)+_{k}^{∑} b_{j,k}2^{j/2}ψ(2^{j}t−k) From the dilation and wavelet equations, get

• the analysis equations

a_{j,k} = _{n}^{∑} h_{0}[n − 2k]a_{j+1,n}
b_{j,k} = _{n}^{∑} h_{1}[n − 2k]a_{j+1,n}

• and the synthesis equation a_{j+1,n} = _{k}^{∑}(h_{0}[n − 2k]a_{j,k} + h_{1}[n − 2k]b_{j,k})

Wavelets and Filter Banks

Key observation here is that these have same form as the outputs of 2-channel 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.

Digital Signals

A digital signal is a quantized function of a discrete variable.

• only have signal values at nΔt, n an integer, Δt the sampling interval

• signal values are elements of a discrete set.

For simplicity, we ignore quantization effects and study 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:

⋯, x[−2], x[−1], x[0], x[1], x[2], ⋯

Filters

Signal processing operations can be described mathematically by means of functions Γ:

Γ : {signals } → {signals }.

We mainly consider operations with the following properties

1. (Linearity) If z[n] = cx[n] + y[n], then

Γ(z)[n] = cΓ(x)[n] + Γ(y)[n]

2. (Shift-invariance) If y[n] = x[n + k] for all n, k fixed, and the outputs are w = Γ(x), z = Γ(y), then

z[n] = Γ(y)[n] = Γ(x)[n + k] = w[n + k]

A Γ satisfying both of these properties is said to be a (linear, shift-invariant) filter.

A “Moving Average” Filter

For any input signal x[n], define

y[n] = Γ(x)[n]
= ( 4 ; ∑ ; k=0 x[n − k])/5

If x[n] = cos (nπ/2),

⋯, 1, 0, −1, 0, 1, 0, −1, 0, 1, ⋯

for instance, then Γ(x)[n] =

⋯, 0, 1/5, 0, −1/5, 0, 1/5, 0, −1/5, 0, ⋯

It is easy to see that this satisfies the linearity and shift-invariance properties, so this defines an filter.

Impulse Response of a Filter

A nice result about filters: Let δ be

⋯, 0, 0, 1, 0, 0, ⋯

(1 at time n = 0) – called a unit impulse input. Γ(δ) is called the impulse response of the filter.

If we know h = Γ(δ), then we also know the output of the filter on any other signal x. Reason: We can write

x[n] = ∞ ; ∑ ; k=−∞ x[k]δ[n − k]

Hence by linearity and shift-invariance

y[n] = Γ(x)[n] = ∞ ; ∑ ; k=−∞ x[k]h[n − k].

The Transfer Function

Package x[n] by its “z-transform”, the series

X(z) = ∞ ; ∑ ; n=−∞ x[n]z^{−n}

Do same for the impulse response h[n] and the output y[n] for input x[n]:

H(z) = ∞ ; ∑ ; n=−∞ h[n]z^{−n}, Y (z) = ∞ ; ∑ ; n=−∞ y[n]z^{−n}

Coefficient of z^{−n} in the product H(z)X(z) is the discrete convolution:

∞ ; ∑ ; k=−∞ h[k]x[n − k],

same as y[n] (!)

⇒ Y (z) = H(z)X(z). H(z) is called the transfer function.

Connection with Fourier Domain

Note:

H(e^{iω}) = ∞ ; ∑ ; n=−∞ h[n]e^{−inω}

is the (discrete-time) Fourier transform of h. So H(z) also determines the frequency response characteristics of the filter on input signals. For instance, with input a pure sinusoid x[n] = e^{inω_{0}} the output y[n] = Γ(x)[n]

= ∞ ; ∑ ; k=−∞ h[k]x[n − k]
= ∞ ; ∑ ; k=−∞ h[k]e^{i(n−k)ω_{0}}
= ( ∞ ; ∑ ; k=−∞ h[k]e^{−ikω_{0}} )e^{inω_{0}}
= H(e^{iω_{0}})x[n]

Frequency Response

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

This is an example of a “low-pass” filter – it “passes” sinusoids of small normalized frequency and “stops” sinusoids of large normalized frequency.

Back to Filter Banks

The h_{0}[n] and h_{1}[n] are coefficients of two filters appearing in the synthesis filter bank:

Synthesis bank:

[figure: a_{j,k}; b_{j,k}; (↑ 2); (↑ 2); H_{0}; H_{1}; a_{j+1,k}]

where (↑ 2) is upsampling, or “interpolation” by 2 (expanding with zeroes).

H_{0} is the synthesis low-pass filter; H_{1} is the synthesis high-pass filter. For example, in Haar case, H_{0} is two-point averaging filter, H_{1} is two-point difference filter.

Filter Banks, cont.

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

[figure: a_{j+1,n}; H̃_{0}; H̃_{1}; (↓ 2); (↓ 2); a_{j,k}; b_{j,k}]

(where H̃_{i}(z) = H_{i}(1/z) and (↓ 2) is downsampling, or decimation by 2)

“Fast Wavelet Transform”

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

• 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:

[figure: V_{j+1}; H̃_{0}; H̃_{1}; (↓ 2); (↓ 2); W_{j}; V_{j}; H̃_{0}; H̃_{1}; (↓ 2); (↓ 2); V_{j−1} ⋯; W_{j−1}]

Wavelet System Design

Can design wavelet systems by designing the analysis (or synthesis) bank filters.

• Often want H_{i} to be FIR filters ⇒ ϕ, ψ have compact support.

• Often want {ϕ(t − k), ψ(2^{j}t − k)} orthogonal in L^{2}(R) – Sufficient: “double-shift orthogonality”:

_{n∈Z}^{∑} h_{i}[n]h_{j}[n − 2k] = δ[i − j]δ[k]  (1)

(condition also implies “perfect reconstruction” – synthesis bank inverts analysis bank).

• Often want good approximation properties – for approximation to order p, sufficient: H_{0}(z) has zero order p ≥ 1 at z = −1 (ω = π). (Condition also related to smoothness of ϕ, ψ.)

The “Cascade Algorithm”

Recall the dilation equation:

ϕ(t) = √2 _{k∈Z}^{∑} h_{0}[k]ϕ(2t − k)

where the coefficients h_{0}[k] in the synthesis low pass filter appear.

The dilation equation says that ϕ(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 ϕ by a fixed point iteration:

ϕ^{(i+1)} = T (ϕ^{(i)})

This is called the “cascade algorithm” for determining the scaling function of a wavelet system.

Convergence of the Cascade Algorithm

For convergence of the cascade algorithm, need an eigenvalue condition: T has λ = 1 as a simple eigenvalue, and all other |λ| < 1.

H_{0}(−1) = 0 is a necessary condition for this to hold.

To Multiwavelets

• There are limitations on what can be accomplished this way, though.

• For instance, a symmetric orthogonal FIR filter can have only two nonzero coefficients.

• To get around these limitations, people have introduced several generalizations of MRA, etc.

• We’ll consider one – multiwavelets.

Multiwavelet Systems

Basic idea here is to allow more than one scaling function, mother wavelet:

• V_{j} = Span {ϕ_{0}(2^{j}t − k), …, ϕ_{r−1}(2^{j}t − k) : k ∈ Z}

• W_{j} = Span {ψ_{0}(2^{j}t − k), …, ψ_{r−1}(2^{j}t − k) : k ∈ Z}

• dilation equation becomes

Φ(t) = √2 _{n∈Z}^{∑} C_{0}[n]Φ(2t − n)

where Φ = (ϕ_{0}, …, ϕ_{r})^{T}, and C_{0}[n] is an r×r matrix for each n. Similarly for wavelet equation.

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 vector-valued signals, we will consider a different use.

For simplicity, take r = 2, and write

C_{0}[n] = ( \begin{matrix}h_{0}[2n] & h_{0}[2n + 1] \\ h_{1}[2n] & h_{1}[2n + 1] \end{matrix} ) ,

then the entries in the C_{0}[n] are split into coefficients from 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.

Multiwavelets for Scalar Signals, cont.

• The 2-channel multiwavelet filter banks can be regrouped into 4-channel scalar filter banks (with interleaving of outputs).

• Condition for orthogonality (analogous to (2) above) becomes

_{n∈Z}^{∑} h_{i}[n]h_{j}[n − 4k] = δ[i − j]δ[k]  (2)

• Natural analog of the “zeroes at z = −1 (ω = π)” condition for approximation order is the “balanced” property: Selesnick, Lebrun and Vetterli, etc.

Balanced Multiwavelet Systems

references:

1. Selesnick, Multiwavelet Bases with Extra Approximation Properties,

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

Both: IEEE Trans. Signal Proc., v. 46, 1998.

Def. The r = 2 multiwavelet system with filters h_{0}, h_{1} as above is balanced of order p if (1+z^{−1}+z^{−2}+z^{−3})^{p} divides H_{0}(z)+α_{p}(z^{4})H_{1}(z), where α_{p} is given by an explicit formula yielding:

α_{1} = 1
α_{2} = \frac{3 − z^{−1}}{2}
α_{3} = \frac{15 − 10z^{−1} + 3z^{−2}}{8} , etc.

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}, …, t^{p−1} can be exactly decomposed (on a finite interval) into linear combinations of the scaling functions ϕ_{i}(t) and their integer translates.

The balancing order is also related to the order of smoothness of the ϕ_{i}(t).