Haar Wavelet Analysis

A First Course in Wavelets with Fourier Analysis

Albert Boggess Francis J. Narcowich
Prentice-Hall, Inc., 2001
 Outlines
Why Wavelet
Haar Wavelets
The Haar Scaling Function
Basic Properties of the Haar Scaling Function
The Haar Wavelet
Haar Decomposition and Reconstruction
Algorithms
Decomposition
Reconstruction
Filters and Diagrams
Summary
 4.1 Why Wavelet
Wavelets were first applied in geophysics to
analyze data from seismic surveys.
Seismic survey

geophones

seismic trace

Sesimic trace
Direct wave (along the surface)
Subsequent waves (rock layers below ground)
 Fourier Transform (FT) is not a good tool
gives no direct information about when an
oscillation occurred.
Short-time FT : equal time interval, high-
frequency bursts occur are hard to detect.
Wavelets can keep track of time and
frequency information. They can be used to
zoom in on the short bursts, or to zoom
out to detect long, slow oscillations
frequency

frequency + time (equal time intervals)

frequency + time
 4.2 Haar Wavelets
4.2.1 The Haar Scaling Function
Wavelet functions
Scaling function (father wavelet)
Wavelet (mother wavelet)
These two functions generate a family of functions
that can be used to break up or reconstruct a signal
The Haar Scaling Function
Translation
Dilation
 Using Haar blocks to approximate a signal

Figure 2

High-frequency noise shows up as tall, thin blocks.

Needs an algorithm that eliminates the noise and not
distribute the rest of the signal.
Disadvantages of Harr wavelet: discontinuous and
does not approximate continuous signals very well.
Chap 6
Dubieties 3

Daubechies 4

Daubechies 8
4.2.2 Basic Properties of the Haar Scaling Function
The Haar Scaling function is defined as

1, if 0 x 1
(x)

0, elsewhere

(x-k) : same graph but translated by to the right (if

k>0) by k units
Let V0 be the space of all functions of the form

a (x k) a
kZ
k k R
 V0 consists of all piecewise constant functions whose
discontinuities are contained in the set of integers
V0 has compact support.

Figure 5

Figure 6

Typical element in V0
f ( x) 2 ( x) 3 ( x 1) 3 ( x 2) 2 ( x 3)
has discontinuities at x=0,1,3, and 4
 Let V1 be the space of piecewise constant functions of
finite support with discontinuities at the half integers
ak ( 2 x k ) ak R
kZ
(2 x)

f ( x) 4 (2 x) 2 (2 x 1) 2 (2 x 2) (2 x 3) V1
has discontinuities at x=0,1/2,3/2, and 2
 Suppose j is any nonnegative integer. The space of step
functions at level j, denoted by Vj , , is defined to be the
space spanned by the set

{ , (2 j x 1), (2 j x), (2 j x 1), (2 j x 2),}

over the real numbers.

Vj is the space of piecewise constant functions of finite
support whose discontinuities are contained in the set
{...,1 / 2 j ,0,1 / 2 j ,2 / 2 j ,3 / 2 j ,...}

V0 V1 ... V j 1 V j V j 1... means no

information is lost as the resolution gets finer. Vj contains
all relevant information up to a resolution scale order 2-j
 A function f(x) belongs to V0 iff f(2jx) belongs to Vj
Proof :
If a function f belongs to Vo , then
f(x) is a linear combinations of { ( x k ), k Z }
Therefore, f(2 j x) is a linear combinations of { ( 2 j x k ), k Z },
which means that f(2 j x) is a member of Vj .
The proof of the converse is similar
 A function f(x) belongs to Vj iff f(2-jx) belongs to V0

Proof :
If a function f belongs to Vj , then
f(x) is a linear combinations of { (2 j x k ), k Z }
Therefore, f(2 j x) is a linear combinations of { (2 j 2 j x k ), k Z },
which means that f(2 j x) is a member of Vo .
How to decompose a signal into its Vj-components

When j is large, the graph of (2jx) is similar to one of

the spikes of a signal that we may wish to filter out.
One way is to construct an orthonormal basis for Vj
using the L2 inner product
k 1
( x k ) dx 1dx 1
2
(x k ) L2
2
- k

if j k , ( x j ) and ( x k ) have disjoint supports

( x j ), ( x k ) L2
( x j ) ( x k )dx 0, j k
-

so the set { ( x j ), k Z } is an orthonormal basis for V0

 Theorem:
j
The set of functions {2 (2 x j ); k Z } is
2 j

an orthonormal basis for Vj

Note : ( (2 j x)) 2 dx 1 / 2 j
-
 4.2.4 The Haar Wavelet
We want to isolate the spikes that belong to Vj, but that are not
members of Vj-1
The way is to decompose Vj as an orthonormal sum of Vj-1 and
its complement.
Start with V1, assume the orthonormal complement of Vo is
generated by translates of some functions , we need:

is a members of V1 and so can be expressed as

(x) l a l (2x - l ) for some choices of a l R
is orthogonal to V0 . This is equivalent to

(x) (x - k) 0 for all integers k

 Harr wavelet ( x) (2 x) (2 x 1)

(x) (2x) - (2(x - 1/2)) (2x) - (2x - 1) (1st : V1 )

1/ 2 1
k 0, (x) (x - k)dx 1dx - 1dx 1/2 - 1/2 0
- 0 1/ 2

k 0, (x) (x - k)dx 0, since supports of (x) and (x) do not overlap
-

(2nd : V0 )
 Any function f1 ak (2 x k ) V1 is orthogonal to V0
k

(linear combinations of ( x l ), l Z ) if and only if

a1 -a0 , a3 -a2 ,.... In this case,
f1 a2 k ( (2 x 2k ) (2 x 2k 1)) a2 k (x k )
kZ kZ

In other words, a function in V1 is orthogonal to V0 if and only if

it is of the form k
a2 k ( x k )
So let W0 be the space of all functions of the form
a (x k ), a
kZ
k k R , (assume only a finite number of a k are nonzero).

W0 is the orthogonal complement of V0 in V1

V1 V0 W0
 Theorem 4.8 (extend to Vj)

Let Wj be the space of functions of the form ak ( 2 j x k ) ak R

k

Wj is the orthogonal complement of Vj in Vj1 and Vj1 Vj Wj

Proof : To establish this theorem, we must show two facts :

1. Each function in Wj is orthogonal to every function in Vj
2. Any function in Vj1 that is orthogonal to Vj must belong to Wj .
Part 1.
Suppose that g kZ ak (2 j x k ) belongs to Wj
and suppose that f belongs to Vj . We must show

g, f L2
g(x)f(x) dx 0
-
Since f(x) belongs to Vj , the function f(2 j x) belongs to V0. So

0
-
kZ k
a ( x k ) f(2 j
x) dx (becasue is orthogonal to V0 )

2 j
-
kZ k
a (2 j
y k ) f( y ) d y (by letting y 2 -j
x)

2
j
g( y ) f( y ) d y .
-

Therefore g is orthogonal to any f Vj .

The second requirement can be derived as V1 case.
 Decomposing Vj

Vj W j Vj1 W j1 W j 2 Vj 2
...
W j1 W j 2 W0 2 V0
So each f in Vj can be decomposed uniquely as a sum
f w j 1 w j 2 w0 f 0
where each wl belongs to Wl , 0 l j 1 and f 0 belongs to V0
Intuitively, wl represents the " spikes" of f of width 1/2l 1
that cannot be represented as linear combinations of
spikes of other widths.
 Theorem:

The space L2 (R) can be decomposed as an infinite orthogonal direct sum

L2 (R) V0 W0 W1 W2 .
In particular, each f L2 (R) can be written uniquely as

f f 0 w j where f 0 belongs to V0 and w j belongs to Wj
j 0
 4.3 Haar Decomposition and Reconstruction
Algorithms

Noise filtering problem

Approximate f by a step function f j V j (for j suitably large)
Decompose f j into its components
f j f 0 w1 w j 1, w l Wl

Example :
Suppose that spikes of width less than 0.01 represents noise
2-6 0.01 2-7 , any w j with j 6 represents noise, set these
commponents equal to zero to filter out this noise.
 Implementation
Step 1 : Approximate the original signal f by a step
function of the form
f j(x) al ( 2 j x l )
lZ
j : small enough to capture the essential features of the signal.
l : range depends on the domain of the signal.
Sample the signal at x ...,-1,-1/2 j ,0,1/2 j ,..., which leads
to al f (l / 2 j ) for l Z .
 Step 2
Decompose (2 j x l ) into its Wl components for l j.
(2x) ( ( x) (x))/2 (4.4)
(2x 1 ) ((x) (x))/2 (4.5)

x 2 j 1 x
The following relations holds for all x R :
(2 j x) ( (2 j 1 x) (2 j 1 x))/2 (4.6)
(2 j x 1 ) ((2 j 1 x) (2 j 1 x))/2 (4.7)
 Example 4.11

j 2, decompose f into W1 ,W0 ,V0 components.

f(x) 2 (4x) 2 (4x 1) (4x 2) (4x 3)
( 4 x) (x)(2 ( x))/2
2
( 4 x 1 ) ( ( 2 x)x))/2
(2
( 4 x 2 ) ( 4( x 1/ 2 )) (x(2( 1/) 2 ) ( x2( 1/))2/2)
( 4 x 3 ) ( 4( x 1/ 2 ) 1) ( (2(x 1/ 2 )(
2( x 1/ ))2 /2)
f(x) [x) (2 ( x)2 ] ([ x)x)
2 (2 ]
[x(2 )1 ( x2 )1 ] / 2 ([ x2 )x 1 (2 ) 1 ] / 2
x(2 )1 2( x)2
W1 -components:x(2 )1
V1 -components: 2 ( 2 x)
V1 V0 ,W0
f(x) x(2 )x) 1 ((x)
 General decomposition scheme
Step 1 :
Divide the sum f j(x) k ak ( 2 j x k ) into even and odd terms :
f j(x) a2 k ( 2 j x 2k ) a2 k 1 ( 2 j x 2k 1 ) (4.9)
kZ kZ

(2 j x) ( (2 j 1 x) (2 j 1 x))/2 (4.6) (2x 1 ) ((2 j 1 x) (2 j 1 x))/2 (4.7)

2 j x 2k 2 j ( x k 21 j )
( 2 j x 2k ) ( (2 j 1 x k) ( 2 j 1 x k )) / 2 (4.10)
( 2 j x 2k 1 ) ( (2 j 1 x k) ( 2 j 1 x k )) / 2 (4.11)
f j(x) a2 k ( ( 2 j 1 x k ) ( 2 j 1 x k )) / 2 a2 k 1( ( 2 j 1 x k ) ( 2 j 1 x k )) / 2
kZ kZ

a2 k a2 k 1 a a2 k 1
( )( 2 j 1 x k ) ( 2 k ) ( 2 j 1 x k )
kZ 2 2
w j 1 f j 1
Wj-1-component Vj-1-component
 Theorem 4.12 (Haar Decomposition)

f j(x) akj ( 2 j x k ) Vj
kZ

Then f j can be decomposed as

f j w j 1 f j 1 , where
w j 1 bkj 1( 2 j 1 x k ) W j 1
kZ

f j 1 akj 1 ( 2 j 1 x k ) V j 1
kZ

a j
a j
a j
a j
with bkj 1 2 k 2 k 1
akj 1 2 k 2 k 1
2 2

f j w j 1 f j 1 w j 1 w j 2 f j 2 ... w j 1 w j 2 w0 f 0
 Example 4.13
Figure 13, j 8, 0 x 1,
decompose f into W1 , W0 , V0 components.
ak8 f( k / 28 ),0 k 28 1
28 1
f 8 ( x ) f( k / 28 )( 28 k)
k 0

V8-component V7-component V6-component V4-component

W7-component
 4.3.2 Reconstruction
Having decomposed f into V0 and Wj' components for
0 j' j, then what ?
Noise filtering : The W j' - components corresponding to the
unwanted frequencies can be thrown out.
Data compression : The W j' - components that are small
can be thrown out, only keep larger bkj' .
bkj' have been modified, need a reconstruction algorithm
to rebuild f(x) alj ( 2 j x l )
lZ

Once this is done, f is a step function of height alj over

the interval l / 2 j x ( l 1 ) / 2 j
 General reconstruction scheme

f(x) f 0 ( x ) w0 ( x ) w j 1 (x) wl Wl where

f 0 ( x ) ak0 ( x k ) V0 and
kZ

wl bkl ( 2l x k ) Wl for 0 l j-1

k

Goal :
To rewrite f ( x ) l alj ( 2 j x l ) and
find a algorithm for the computation of constants alj
 General reconstruction scheme
(x) ( 2 x) ( 2 x 1 ) (4.12)
(x) ( 2 x) ( 2 x 1 ) (4.13)
( 2 j 1 x) ( 2 j x) ( 2 j x 1 ) (4.14)
( 2 j 1 x) ( 2 j x) ( 2 j x 1 ) (4.15)
Using 4.12 with x replaced by x-k , we have
f 0 ( x) ak0 ( x k ) (ak0 (2 x 2k ) ak0 (2 x 2k 1))
kZ kZ

so f 0 ( x) al1 (2 x l ) (4.16)
kZ

ak0 if l 2k is even
where al1
k if l 2k 1 is odd
0
a
Similarly , w0 k ak0 ( x k ) can be written as
w0 ( x) bl1 (2 x l ) (4.17)
lZ

bk0 if l 2k is even
where bl1
bk if l 2k 1 is odd
0

f 0 ( x) w0 ( x) al1 (2 x l ) ( f1 ( x))
lZ

ak0 bk0 if l 2k is even

where al1 al1 bl1
a
k
0
bk if l 2k 1 is odd
0
Add w1 k bk1 (2 x k ) to the sum
f 0 ( x) wo ( x) w1 ( x) al2 (2 2 x l ) ( f 2 )
lZ

a1k bk1 , if l 2k is even

where al2
a
k
1
bk , if l 2k 1 is odd
1

al0 , bl0 coefficients determine the al1 - coefficients.

al1 , bl1 coefficients determine the al2 - coefficients,
and so on in a recursive manner.
 Theorem 4.14 (Haar Reconstruction)

Suppose f f 0 w0 w j 1 with
f 0 ( x ) ak0 ( x k ) V0 and
kZ

w j' ( x ) bkj' ( 2 j' x k ) W j' for 0 j' j .

kZ

Then f ( x ) alj ( 2 j x l ) V j
lZ

where the alj' are determined recursively for j' 1,

then j' 2, and so on until j' j, by the algorithm
a j' 1
b j' 1
, if l 2k is even
al j' 1
j' k k
j' 1
a
k bk , if l 2k 1 is odd
 Example 4.15

sample signal 80% compression 90% compression

Figure 19, j 8, 0 x 1,
ak8 f (k / 28 ),0 k 28 1
f f 0 w0 w1 w2 w7 with
f 0 ( x) a00(x) V0 and w j ' ( x) bkj
'
(2 j ' x k ) W j '
kZ

for 0 j ' 7.
After ordering the bkj by size,
set the smallest 80% equal to zero (80% compression),
then reconstruct the signal.
 4.3.3 Filters and Diagrams
Decomposition algorithm
Discrete filters (convolution operators) H and L
a2jk a2jk 1 a2jk a2jk 1
via their impulse responses, bk
j 1

2
j 1
ak
2
n 1
which are the sequences h and : [ y * z ]k y j z k j
j 0

1 1 1 1
h ( 0 0 ), ( 0 0 )
2 2 2 2
k=-1,0
k=-1,0

If {x k } 2 , then H(x): h x and L(x) x.

The resulting sequences are thus
1 1 1 1
H(x)k (h x)k xk - xk 1 L(x)k ( x)k xk xk 1
2 2 2 2
If we keep only even subscripts, then
1 1 1 1
H(x)2 k (h x)2 k x2 k - x2 k 1 and L(x)2 k ( x)2 k x2 k x2 k 1 .
2 2 2 2
This operation of discarding the odd coefficients in a sequence
is called downsampling, denoted by D.
Go from level - j scaling coefficients akj to get the level j - 1
sacling and wavelet coefficients.decompose f into V0 and W j '
components for 0 j' j bkj 1 DH (a j ) k and akj 1 DL(a j ) k .
a j
a j
a j
a j
(bkj 1 2 k 2 k 1
akj 1 2 k 2 k 1
Theorem 4.12)
2 2

downsampling
operator
 Reconstruction

~ ~
Discrete filters (convolution operators) H and L
via their impulse responses, a j' 1
b j' 1
, if l 2k is even
~ ~ al j' 1
j' k k
j' 1
which are the sequences h and : a
k bk , if l 2k 1 is odd
~
h ( 0 1 - 1 0 ), ( 0 1 1 0 )

k=0,1 k=0,1
~ ~
For a sequence {x k }, we have (h x) k xk -xk-1 and ( x) k xk xk-1.
If x and y are sequences in which the odd entries are all 0, then
~ x 2k l 2k is even ~ y 2k l 2k is even
(h x)l ( y )l
- x 2k l 2k 1 is odd y 2k l 2k 1 is odd
 ~ ~
Adding the two sequences h x and y then gives us
~ ~ x2k y2k l 2k is even
(h x) l ( y )l
y2k -x2k l 2k 1 is odd
This is almost the pattern for the reconstructin formula given in
j ' 1 j ' 1
ak bk , if l 2k is even
Theorem 4.14. al j '1
j'
j ' 1
a
k bk , if l 2k 1 is odd
Although we have assumed that x2k 1 ' s and y2k 1 ' s are 0,
the x2k ' s and y2k ' s are ours to choose,
so we set x2k bkj 1 and y2k akj 1 ; that is,
x ( 0 b j-1
-1 0 b
0
j-1
0 b
1
j-1 j-1
0 b
2 0 ) and similarly for y.
The sequences x and y are called upsampling of the
sequences b j 1 and a j 1
We use U to denote the upsampling operator, so x Ub j 1 and y Ua j 1
~ ~
a j L Ua j 1 HUb j 1

upsampling
operator
 Summary
Given signal y f(t),
Let and be the Haar scaling function and wavelet
Step1. Sample
continuous signal : choose top level j J and 2 J is larger than
the Nyquist rate.
Let akJ f (k / 2 J ) k : a finite interval determined by the duration
of the signal.
ex. 0 t 1, 0 k 2 J -1 (or 1 k 2 J )
discrete signal : skip this step.
highest - level approximation to f
f J (x) akJ (2 J x k)
kZ
Step 2 Decomposition
Decompose f J wJ 1 w j 1 w j w0 f 0 where
w j 1 bl j 1( 2 j 1 x l ) f j 1 alj 1 ( 2 j 1 x l )
lZ lZ

The coefficients bl j 1 , alj 1 are determined from

the a j recursively by the algorithm
alj 1 DL( a j )k (4.18) bl j 1 DH ( a j )k (4.19)
where H and L are the high - and low pass filters.
When j J , akJ 1 and bkJ 1 are determined by akJ , which are the sampled
signal values from Step 1. Then j becomes J-1 and akJ 2 and bkJ 2 are
determined from akJ 1. Then j becomes J-2 , and so on, until the level
reaached is either satisfactory for some purpose or there are two few
coefficients to continue. Unless otherwised stated, will continue until
the j 0 level is reached.
Step 3 Processing. After decomposition, the signal in the form
J 1
f J ( x) w j f 0 (4.20)
j 0
J 1
( bkj 1 (2 j 1 x l )) ak0 ( x k ) (4.21)
j 0 kZ kZ

The signal can now be filtered by modifying the wavelet

coefficients bkj .
a. Filter out all high frequencies : all the bkj would be set to
zero for j above a certain value.
b. only a certain segment of a signal corresponding to particular
values of k is to be filtered.
c. data compression : the bkj that are below a certain threshold
would be set to zero.
Step 4. Reconstruction
Take the modified signal, f J
J 1
f J ( x) ( bkj 1 (2 j 1 x l )) ak0 ( x k ) (4.21)
j 0 kZ kZ

and reconstruct it as
f J akJ ( 2 J x k ).
kZ

This is accomplished by the reconstruction algorithm

~ j 1 ~ j 1
a L Ua HUb
j
(4.22) for j 1,...,J .
When j 1, the a1k are computed from the ak0 and bk0 .
When j 2, the ak2 are computed from the a1k and bk1 and so forth.
The range of k is determined by the time duration of the signal.
When j J (the top level), akJ represents the approximate value
of the processed signal at time x k/2 J .