Transcript 1. dia

WAVELETS AND
FILTER BANKS
Ildikó László, PhD
ELTE UNIV., BUDAPEST,
HUNGARY
[email protected]
NEWER VERSION BUT STILL NEEDS
IMPROOVMENTS AND CORRECTIONS!!!
BIBLIOGRAPHY:
- G. Strang, T. Nguyen:Wavelets and
Filter Banks, Wellesly-Cambridge Press
- F. Schipp, W.R. Wade:Transforms on
Normed Fields
- Ingrid Daubechies: Ten Lectures on
Wavelets
- Stephane Mallat: A Wavelet tour of
Signal Processing
- Charles K. Chui: An Introduction to
Wavelets etc.
AT THE END YOU WILL HAVE TO BE
ABLE TO „GET” THINGS LIKE:
Some rows (10...) of A4:
50
100
150
200
250
50
100
150
200
250
50
100
150
200
250
Some rows (10...) of A3:
50
100
150
200
250
50
100
150
200
250
50
100
150
200
250
Some rows (22...) of A4:
50
100
150
200
250
50
100
150
200
250
50
100
150
200
250
Some rows (48...) of A3:
50
100
50
100
50
100
150
150
200
250
200
250
150 200
250
1.INTRODUCTION
- the classical Fourier analysis where a signal is
represented by its trigonometric Fourier
transform, is one of the most widely spread tools
in signal and image processing;
- at the end of the 19th century Du Bois-Reymond
constructed a continuous function with divergent
Fourier series;
- Hilbert: whether there exist any orthonormal
system for which the Fourier series with respect to
this system do not posses this singularity?
- Alfred Haar in 1909 constructed such an
orthonormal system for which the Haar-Fourier
series of continuous functions converge uniformly:

1
,



 (t )   1 ,

0 ,


1
0t 
2
1
 t 1
2
otherwise
- This was the first wavelet; (1910, Szeged, PhD)
-The basic goal of Fourier series is to take a
signal – considered as a function of time
variable t, and decompose it into its various
frequency components;
-The basic building blocks are sin(nt) and cos(nt),
which vibrate at a frequency of n times per 2
interval.
- Consider the following function:
f (t )  sin(t )  2 cos(3t )  0.3 sin(50t ).
- This function has three components that
vibrate at frequency 1 – the sin(t) part,
at frequency 3 – the 2cos(3t) part,
at frequency 50 – the 0.3sin(50t) part;
- we can express a function f(t) in terms
of the basis functions, sine and cosine:
f (t )   an cos(nt)  bn sin(nt)
n
-A trigonometric expansion is a sum of the
form:
f ( x)  a0   ak cos(kx)  bk sin(kx)
k
where the sum could be finite or infinite.
-One disadvantage of Fourier series is that its
building blocks, sines and cosines, are periodic
waves that continue forever.
- In many appl., given a signal f (t ) one is
interested in its frequency component locally
in time;
(similar to music notation, which tells the player
which notes –frequency inf. –to play at any given
moment)
- The standard Fourier transform,
1
Ff   
2

f (t )e
it
also gives the frequency content of
dt
f (t )
,.
- but inf. concerning time-localization of
e.g., high frequency bursts cannot be read off
easily from Ff
- time localization can be achieved by first
windowing the signal f .
f(t)
1
g(t)
0
(T
win
f )(t )   f ( s ) g ( s  t )e
 is
- which is the windowed Fourier transform;
- cuts off only a well-localized slice of f ;
-The wavelet transform provides a similar
time-frequency description, with a few
important differences;
ds
Fourier transform:
- represents the frequency comp. of a signal
- doesn‘t offer localization in time
Wavelet transform:
- cuts up a signal into frequency components
- studies each component with a resolution
matched to its scale
- offers localization in time
- A wavelet  is a function of zero average:


(
t
)
dt

0


- which is dilated by j
and translated by k;
- j is the scaling parameter;
- k is the translation coeff.; allows us to move
the time localization centre;
f is localized around k;
Wavelets - compactly supported small waves
(don‘t extend from –infty to +infty)
Each wavelet - is built up from the same
„MOTHER“ wavelet by
translation and dilation
 ( x)  2  (2 x  k )
j
k
wjk(x) =
j/2
j/2
j
2 w(2 x
j
- k)
- the wavelet transform of f at the scale j and
position k is computed by correlating f with
the wavelet:

WTf ( j, k ) 
 f (t )2
j
2
 (2 t  k )dt
*

- have only recently been used –
1988 Ingrid Daubechies
j
Wavelets:- basis functions
- linearly indep. functions;
Fourier and Wavelet transforms:
- representation of a signal
f(t) by basis functions
f (t )   F ( )e d
it
f (t )   b jk w jk (t )
j ,k
- The Fourier transform of a complex, two dimensional
function f(x,y) is given by:

F( fx, f y ) 
  f ( x, y) exp[i2 ( f x  f
x
y
y)]dxdy

- where f x and f y are referred to as frequencies;
- the inverse Fourier transform

f ( x, y) 
  F( f , f
x

y
) exp[i 2 ( f x x  f y y)]dfx df y
- that is, f(x,y) is a linear combination of elementary
functions, where the complex number F ( f f )
x, y
is a weighting factor;
– when dealing with linear
systems – this can be used for decomposing a complicated
input signal into more simple inputs;
- that is, the response of the system can be calculated as
the superposition of the responses given by the system to
each of these “elementary” functions of the form:
exp[i2 ( f x x  f y y)]
-
Examples;
Fourier transform theorems;
Linear Systems. Invariant Linear Systems;
Sampling theory;
Let us consider a rectangular lattice of samples of the
function g(x,y) as defined:
g_s(x,y)=comb(x/X) comb(y/Y)
y
x
Y
X
-The sampled function g consists of an array of 
functions, spaced at intervals of width X in the x
direction and Y in the y direction.
-The area under each  function is proportional to
the value of the g function at that particular point.
-The spectrum G_s of g_s can be found by convolving
the transform of the comb function with the transform
of g.
G(fx, fy ) fy
G(fx, fy )
fy
fx
fx
1/X
1/Y
1.2
1.0
X= -1
-1.2
1.2
1
T=
-1
-1.2
y=T*x=
2.2
0.2
-2.2
0.2
x = T^(-1) *y =
1
1
0
0
1
-1
0
0
0
0
1
1
0
0
1
-1
=
1.1
1.1
-1.1
-1.1
Compressed, reconstructed; y(1)=y(3)=0
• Convolution – example:
n = 0, 1, 2
k = 0, 1, 2
n=0 x(2) x(1) x(0)
y(0)=h(0)x(0)
h(0) h(1) h(2)
n=1 x(2) x(1) x(0)
y(1)=h(0)x(1)+
h(0) h(1) h(2)
h(1)x(0)
n=2 x(2) x(1) x(0)
h(0) h(1) h(2)
y(2)=h(0)x(2)+h(1)x(1)+h(2)x(0)
FILTERS
- A Filter – a linear time-invariant operator;
- can be characterized by its impulse
response function;
- acts on an input vector x;
the output y:
y(n)   h(k ) x(n  k )  h * x.
k
is a convolution of x and the impulse response
of the system;
- let us consider the input signal x(n) with the
pure frequency  :
in
x(n)  exp
- then the output in the time domain is:
y(n)   h(k ) exp
i ( n k )
in
 exp
k
 h(k ) exp
ik
k
- where the last term can be recognized as:
the Fourier transform of the impulse response
of the system:
ik
H ()   h(k ) exp
k
- the Fourier transform of the output signal:
Y ( )   y(n) exp
 in

n
 h(k ) exp
i ( n  k )
n
k
 X ( ) H ( )
in
exp

What is the connection between:
WAVELETS, FILTERS and FILTER BANKS?
- it is the High_Pass that leads to w(t)
- the Low_Pass filter leads to a scaling
function:  (t ).
- the scaling function in continuous time
comes from an infinite repetition of the
lowpass filter, with rescaling at each repetition;
- the wavelet follows from (t )by just one
application of the highpass filter.
- averages come from the scaling functions;
- details come from the wavelets.
signal at level j (local averages)
+
=signal at level (j+1)
details at level j (local differences)
•LOW-PASS FILTER –or MOVING
AVERAGE
y(n)   h(k ) x(n  k ).
k
To build up the simplest lowpass filter,
we use the Haar filter coefficients:
h(0)=1/2 and h(1)=1/2.
- the output at time t=n is the average of the
input x(n) and that at time t=n-1 : x(n-1).
y(n)=1/2 x(n)+1/2 x(n-1)
averaging filter=1/2 (identity) + 1/2 (delay)
Every linear operator acting on a signal
x can be represented by a matrix: y=H*x
.
y(-1)
y(0)
y(1)
.
=
½
½
0
0
0
0
½
½
0
0
0
0
½
½
0
0 0
0 0
0 0
½0
. .
.
x(-1)
x(0)
x(1)
.
- with the simple input signal, with pure
frequency: x(n)  expin
in
i ( n 1)
1
1
y ( n) 
exp 
exp

2
2
 i
in
1
1
( 
exp ) exp .
2
2
- that is, the output can be written as:
y(n)=SOMETHING * x(n)
- this “something” is expressing the effect of
our filter, or system on the input signal;
- The frequency response or transfer function
of the system:
 i
1
1
H ( )  (  exp ).
2
2
- because this is a periodic function, we want:
i( )
H ()  H () exp
- it results:
 i
H ( )  (cos ) exp
2
2
.
.
- the first term is the amplitude, or magnitude
of the transfer function:
H 0( )  cos
2
- the second term is its phase:

( )  
2
Low-Pass Filter –H_0 (Frequency response
function, or transfer function)
IMPORTANT: our filters are CASUAL, that is
the output cannot arrive earlier than the input!!!
1
|H |
_0
pi
-pi
0
pi
The magnitude
-pi
0
The phase of H_0
High-PASS FILTER –or MOVING
DIFFERENCE
Ex. with the Haar coeff.: h(0)=1/2; h(1)=-1/2
The output at t=n:
y(n)=1/2 x(n)-1/2 x(n-1)
.
y(-1)
y(0)
y(1)
.
=
½ 0 0
-½ ½ 0
0 -½ ½
0 0 -½
0 0 0
0 0
0 0
0 0
½0
. .
.
x(-1)
x(0)
x(1)
.
The number of the coefficients is finite!
It results: the impulse response is also finite – FIR!
- a highpass filter takes differences, i.e.:
picks out the bumps in the signal;
difference filter=1/2(identity)-1/2(delay).
in
i ( n 1)
1
1
y ( n) 
exp  exp
2
2
i
in
1
1
 (  exp ) exp
2
2
in
 H1 ( ) exp .
- the quantity:
 i
1
1
H1 ( )  (  exp )
2
2
- is the highpass response, or the transfer
function of the highpass filter.
- it results:
 i
H1 ( )  (sin  )i exp
2
2
.
- the magnitude of this periodic highpass
filter is:
H1 ( )  sin  .
2
- the second term is the phase of the highpass
transfer function:
i ( )
exp
i

 i exp 2
 
 i
2

i
exp

for      0
for 0     .
High-Pass Filter – H
|H |
_1
_1
1
pi/2
0
-pi
0
pi
-pi
pi
The magnitude
and phase of the
High-Pass filter
- although there is a discontinuity in the
phase at:   0, we still say, that:
this is a linear phase filter;
H_0 (  )
H_1 ( )
C(  ) = 2 H_0( )
D(  ) = 2 H_1( )
FILTER BANK: LOW_PASS and HIGH_PASS
+DOWNSAMPLING
Input signal: x(n)
 a j 1k j 1k
k
Output at level j=1
Scaling coeff. ajk
C
2
D
2
Wavelet coefficients
at level j=1 : bjk
-THE OUTPUT OF THE LOW_PASS FILTER IS
NOT A FINAL OUTPUT
- the result: an average of the original signal
-THIS IS THAT OUTPUT WHICH CAN (WILL) BE
FILTERED AGAIN :
-GOES TO THE NEXT LEVEL:
- filtered by the lowpass:
-results in an even coarser average
of the original signal;
- filtered by the highpass:
- fine details of the previously
averaged signal;
-the output of the HIGH_PASS FILTER IS
FINAL OUTPUT
– WILL NOT BE FILTERED AGAIN;
- these are the WAVELET
COEFFICIENTS at the FINEST SCALE:
- HIGHEST FREQUENCY COMPONENTS!!!
N
=256
Some rows (139...) of A1:
N/2
L_P
128
H_P
50
100
150
200
250
A1
Some rows (48...) of A3:
32
L 0
H
0 1
0
15
30
45
0
1
60
A3
N/2
-the filter bank built up with the Haar
coefficients - IS ORTHOGONAL
- scaling functions come from the iteration of the
LOW_PASS FILTER
-they are RESCAILED at each ITERATINON
- this results in COARSER AND COARSER
AVERAGING OF THE INPUT SIGNAL
Some rows (10...) of A3:
15
30
45
60
L_P:1-32
65
Some rows (5...) of A4:
5
10
15
20
L_P:1-16
25
c(0)=c(1)= 2*h(0)
L=( 2) C = 2*
½
0
0
0
½
0
0
0
0 0
½½
0 0
0 0
0
0
½
0
0
0
½
0
0 ..
0 ..
0 ..
½. .
d(0)= 2/2, d(1)= - 2/2
B=( 2) D = 2*
-½
0
0
0
½
0
0
0
0 0 0
-½ ½ 0
0 0 -½
0 0 0
0 0 ..
0 0 ..
½ 0 ..
0 -½ . .
H1
Because of
causality the
upper triangle
consists of
zeros only!
THIS MATRIX REPRESENTS THE WHOLE
ANALYSIS BANK
( 2) C
( 2) D
=
L
B
= 1/ 2
1 1 0 0 0 0 .
0 0 1 1 0 0 .
.
-1 1 0 0 0 0 .
0 0 -1 1 0 0 .
.
.
.
.
.
.
.
L_P
H_P
THE ORTHOGONAL HAAR SYSTEM - FB
2 H0
v0 = ( 2) C x = L x
C
2
D
2
X
2 H1
v1 = ( 2) D x = B x
-The Analysis Matrix, or Analysis Bank built
up by the Haar coefficients is an orthogonal
Filter Bank;
-For an orthogonal system, the inverse matrix=
the transpose;
-That is: S=A-1=AT
THIS MATRIX REPRESENTS THE WHOLE
SYNTHESIS BANK
-1
L
B
=
T
L
T
B
= 1//2
1
1
0
0
0
0
1
1
N/2
.
.
.
.
.
-1
1
0
0
0
0
-1
1
.
.
.
.
.
N/2
WHEN THE MATRIX IS ORTHOGONAL:
THE INVERSE IS THE TRANSPOSE
DOWNSAMPLING
.
x(0)
( 2) x(1)
x(2)
.
=
.
x(0)
x(2)
x(4)
.
-
UPSAMPLING
.
x(0)
( 2) x(2)
x(4)
.
=
x(0)
0
x(2)
0
x(4)
- Repr. of the down- sampling op. by matrix
multiplication:
- as if every second row of a unitary
matrix would be missing!
100000
001000
000010
. . . . . .
.
x(0)
x(1)
x(2)
.
=
.
x(0)
x(2)
x(4)
.
DOWNSAMPLING y(n) = ( 2) x(n) = x(2n)
x(n)
UPSAMPLING
u(n) = ( 2) v(n)
v(n)
0 123 4 5 6 7
0 123 4 5 6 7
u(n)
y(n)
0 123 4 5 6 7
0 123 4 5 6 7
1
2
1 1 0 0 0 0 0 0
0 0 1 1 0 0 0 0
0 0 0 0 1 1 0 0
0 0 0 0 0 0 1 1
1-1 0 0 0 0 0 0
0 0 1 -1 0 0 0 0
0 0 0 0 1 -1 0 0
0 0 0 0 0 0 1 -1
x(-1)
x(0)
x(1)
x(2)
x(3)
x(4)
x(5)
x(6)
=
1
2
x(0)+x(-1)
x(2)+x(1)
x(4)+x(3)
x(6)+x(5)
x(0)-x(-1)
x(2)-x(1)
x(4)-x(3)
x(6)-x(5)
Low_Pass part
High_Pass part
v0 = ( 2) Cx = Lx
Lx = 1/ 2
.
x(0)+x(-1)
x(2)+x(1)
x(4)+x(3)
.
v1 = ( 2) Dx = Bx
Bx = 1/ 2
.
x(0)-x(-1)
x(2)-x(1)
x(4)-x(3)
.
RECONSTRUCTION
u0 =( 2)v0
x(0)+x(-1)
0
u0 = 1/ 2 x(2)+x(1)
0
x(4)+x(3)
.
- the Low_Pass part
of the synthesis matrix;
u1 =( 2)v1
u1 = 1/ 2
x(0)-x(-1)
0
x(2)-x(1)
0
x(4)-x(3)
.
- the High_Pass part
of the synthesis matrix;
FILTERING AFTER UPSAMPLING
RECOVERS THE INPUT DELAYED: x(n-1)
x(0)+x(-1)
0
F filters 1/ 2 x(2)+x(1) to give 1/2
0
Low_Pass
x(4)+x(3)
Of the
.
Synthesis
x(0)+x(-1)
x(0)+x(-1)
x(2)+x(1) = w0
x(2)+x(1)
x(4)+x(3)
.
x(0)-x(-1)
0
G filters 1/ 2 x(2)-x(1) to give
0
High_Pass
x(4)-x(3)
For the
.
Synthesis
w0
x(n-1)
w1
1
2
-x(0)+x(-1)
x(0)-x(-1)
-x(2) +x(1) = w1
x(2)-x(1)
x(4) -x(3)
.
w_0(n) = ½ (u_0(n) + u_0(n-1)
w_1(n) = ½ (u_1(n) - u_1(n-1)
x(0)+x(-1)
x(0)+x(-1)
1/2 x(2)+x(1) = w0
x(2)+x(1)
x(4)+x(3)
.
w0
x(n-1)
w1
-x(0)+x(-1)
x(0)-x(-1)
1/2 x(1) -x(0) = w1
x(2)-x(1)
x(3) -x(2)
.
w0(n) + w1(n) =1/2[x(0)+x(-1) + x(0)-x(-1) ]=x(n-1)
ANALYSIS
x(n)
SYNTHESIS
L
2
2
F
H
2
2
G
 =x(n-1)
Discret Wavelet Transform - DWT
Signal at level j
Signal at
level j-1
 a j 1k j 1k
k
L
2
H
2
Wavelet coefficients
At level j : bjk
Fast Wavelet Transform - FWT
Sign. at
j+1:aj+1,k
a
Signal at
lev.: j-1
j,k
L
H
L
2
H
2
2
2
Wavelet coeff.
at lev. j: bj,k
Wavelet coeff.
at j+1 : bj+1,k
N
A =
…
L
H
0
L
N/2
0
1
H
N/2
A2
A1
THE TRANSFORM OF A SIGNAL –
IN THE FREQUENCY DOMAIN
y(n) = h(k)x(n-k)
Y( ) = X( ) H( )
THE EFFECT OF THE DOWNSAMPLING
IN THE FREQUENCY DOMAIN
Y( ) = ½ {X( /2) + X( /2 + )} = V()
DOWNSAMPLING y(n) = ( 2) x(n) = x(2n)
x(n)
UPSAMPLING
u(n) = ( 2) v(n)
v(n)
0 123 4 5 6 7
0 123 4 5 6 7
u(n)
y(n)
0 123 4 5 6 7
0 123 4 5 6 7
- if our input signal is:
then the output is:
in
x(n)  exp
y(n)  x(2n)
- ? the effect of the downsampling in the
frequency domain is?
Y ( )   y(n) exp
 in
n
in
x
(
2
n
)
exp

Y ( ) 
n

 x(m) exp
i m 
2
"m"even


1
2
im
 (1  (1) ) x(m) exp
m
2
m

im
im(  ) 
2
1  x(m) exp 2   x(m) exp

2 m
m

- From where it results:
Y ( )  V ( )  1

X ( )  X (   ).
2
2
2
- that is, two input signals, that of frequency

and   
2
2
give the same output!
- This means that the recovery of the input
signal from this output wouldn’t be possible.
x(n)
0 123 4 5 6 7
u(n)
y(n)
0 123 4 5 6 7
- again: downsampling
in time domain;
0 123 4 5 6 7
X( )
-2
-
X( /2)
1
1
2

-2
-
X( /2+ )

2
Period: 4 
Period: 2 
Thus the result of this output, V( ) would be:
1
1/2
-2
-

That is, a constant!
We could not
recover the original
2 signal from this!
X( )
-2
1
-
2

X ( )1
2
-2
-
Aliasing
X (   )
2

2
- a high frequency component:   
2
overlaps the low frequency
component:  2 .
FIRST FILTER THEN DOWNSAMPLE
X( )
X( )
1
1
IS FILTERED
-2
-

2
-2
- -/2  /2 
2
  2
X( /2)
-2
-
X( /2+ )
1/2

2
-do not overlap
1
-2
X( )
- -/2  /2 
X( /2)

2
-2
-
X( /2+ )
1/2

2
V ( )  1 X ( )  X (   )
2
2
2

- as can be seen, if the signal is filtered before
downsampling, that is, is band limited, then
the alias doesn’t overlap the original signal;
- that is, the original signal can be recovered;
-RECONSTRUCTION
- the output of the filtered and downsampled
signal, Y ( ) is denoted by: V ( )
- when we analyze a signal, first comes filtering,
then downsampling;
- when we reconstruct a signal, first comes
upsampling, then filtering;
- in the time domain the upsampled signal
is denoted by u and takes the form:
 u ( 2k )  v ( k )
u  ( 2)v  
u (2k  1)  0.
- in the frequency domain this becomes simple,
as only the even indices enter the sum:
U ( ) 
 u (n) exp
 in
n
  v(k ) exp
 ik 2
k
  u (2k ) exp
 i 2 k
k
 V (2 ).
y(n)=v(n)=x(2n)
x(n)
0 123 4 5 6 7
u(2k)=v(k)
u(2k+1)=0
u(n)
y(n)
0 123 4 5 6 7
0 123 4 5 6 7
- downsampling and upsampling (and
filtering) in the time domain
UPSAMPLING after DOWNSAMPLING
- the analysis bank includes downsampling;
- the synthesis bank includes upsampling;
-the two operations:
v  ( 2) x then u  ( 2)v give u  ( 2)( 2) x.
- in the frequency domain:


1 X ( )  X (   )
2
2
2
U ( )  V (2 )  1 X ( )  X (   ).
2
V ( ) 
V( ) =½ {X( /2) + X( /2+ )}
X( )
1
-2
- -/2  /2 
X( /2)
2
-2

-
U( ) =V(2 ) =½ {X() + X( + )}
X( ) 1/2

-/2

X( + )
/2

X( /2+ )
1/2
2
AT RECONSTRUCTION:
u =( 2)v
U( ) =V(2 ) =½ {X( ) + X( + )}
IMAGING: ONE INPUT – TWO OUTPUTS!
 WAS
IF
1
-2
FILTERED BEFORE ( 2):
X( )
- -/2  /2 
X( /2)
2
-2
-
X( /2+ )
1/2

V( ) =½ {X( /2) + X( /2+)}
2
U( ) =V(2 ) =½ {X() + X( + )}
X( ) 1/2

-/2
X( +
/2
)

THE OUTPUTS V( ) AND U(  )FROM
DOWNSAMPLING AND UPSAMPLING STILL
CONTAIN THE ALIASING PART. BUT A NONOVERLAPPING ALIAS CAN BE FILTERED AWAY.
- SAMPLING OPERATIONS IN THE
z-DOMAIN:
- denoting: z
of x(n):
i
 exp
the Fourier transform
X ( z )   x(n) z .
n
2 X ( z
n
V ( z) 
1
1
2
U ( z)  V ( z )  1
2
1

)  X ( z 2 )

X ( z )  X ( z ).
2
SCALING FUNCTIONS AND WAVELETS
Dilation eq. for the scaling function:
Φ(t) =2 h_0(k) Φ(2t-k)
Wavelet eq.:
w(t) =2  h_1(k) Φ(2t-k)
Then, for the Haar case we get:
Φ(t) = Φ(2t) + Φ(2t-1)
w(t) = Φ(2t) - Φ(2t-1)
Scaling Function - Haar coeff.: h0 , h1
Φ(t-1)
Φ(t)
1
1
0
1
0
1
Φ(2t)
1
1/2
0 1/2
2
1
Φ(2t-1)
j=1
-Φ(2t-1)
1
0
1
j=0
k=1
0
1/2
1
1
Φ(t) = Φ(2t) + Φ(2t-1)
Φ(t-1)
j=0
k=1
1
0
1
1
1
2
Φ(2t-2)
1
0
(t  1)  (2t  2)  (2t  3)
2
1 1.5
1
2
1
0
Φ(2t-3)
2k=2;
2k+1=3
j=1 k=1
1 1.5 2
1
w(t)
1
1
0 1/2
w(2t)
w(2t-1)
1
1
1/2
0
0
w(t) = Φ(2t) - Φ(2t-1)
w
(
t
)
w
(
2
t
)

0

1/2
The Haar Wavelets
w
(
t
)
w
(
2
t

1
)

0

!!!That is, the Haar wavelet basis is an orthogonl basis.
1
1
2
2
Φ(2t-2)
1
1
0
0
1 1.5
1
w(t-1)
1
0
j=1 k=1
1 1.5 2
w(t  1)  (2t  2)  (2t  3)
2
1
Φ(2t-3)
Averages (lowpass filter) - scaling coefficient:
a j 1,k
1

(a j , 2 k  a j , 2 k 1 )
2
Differences (highpass filter) – wavelet coefficient:
b j 1,k
1

(a j , 2 k  a j , 2 k 1 )
2

(
t
)

(
t

k
)


(
k
)


1 if k=0
0 otherwise
w
(
t
)
w
(
t

k
)


(
k
)

Haar scaling functions and wavelets at the same
scaling level (same j) do not overlap. They are
orthogonal.
w
(
t
)
w
(
2
t

1
)

0

The rescaled Haar wavelets:
wjk(x) =
j/2
j
2 w(2 x
- k)
form an orthonormal basis:

w
(
t
)
w
(
t
)
dt


(
j

J
)

(
k

K
).
jk
JK


Generally (for the Haar):
N
(t )  2 h0 (k )(2t  k )
k 0
N
w(t )  2 h1 (k )(2t  k )
k 0
Different forms of the scaling function or –
DILATION EQUATION:
N
(t )  2 h0 (k )(2t  k )
k 0
-with the original coeff.
N
(t )  2  c(k )(2t  k )
k 0
- after downsampling
The wavelet equation: (wavelets from filters)
N
w(t )  2 h1 (k )(2t  k )
k 0
N
w(t )  2  d (k )(2t  k )
k 0
Analysis of a function


b jk 
f (t )w jk (t )dt


Synthesis of a function
f (t )   b jk w jk (t )
N
j ,k
(t )   a jk (k )(2t  k )
k 0
Perfect reconstruction:
- we will have causal filters only, i.e.:
- at the output
xˆ[n]  x[n  1]
- PR with delay;
- Perfect reconstruction means, that the synthesis
bank is the inverse of the analysis bank;
1
T
- in the Haar ex.:
SA A
the synthesis bank is the transpose of the
analysis bank;
- Two channel FB – Haar ex.:
H 0 ( z) and H1 ( z)
are generally low_pass and high_pass filters, but
they are not ideal;
1 H 0 ( )
H1 ( )



2
0
H1 ( )

2


- aliasing
- Downsampling operation can produce aliasing;
- if X ( )  1, that is the input has all frequencies,
than: after downsampling, at the output of the
low_pass channel we have:
V ()  {H0 ( 2 ) X ( 2 )  H0 ( 2   ) X ( 2   )}
1
2




- in the z domain:
1
2
1
2
1
2
1
2
V ( z)  12 {H 0 ( z ) X ( z )  H 0 ( z ) X ( z )}
- where:
X ( z )  X (   )
1
2
X ( z )  X ( 2 )
- the goal is, to design the synthesis filters in such a
way, that there is no aliasing;
- only an overall delay to be allowed;
i
z  exp
X ( z )   x(n) z .
n

n

1 X ( )  X (   )
2
2
2
U ( )  V (2 )  1 X ( )  X (   ).
2
V ( ) 
V( ) =½ {X( /2) + X( /2+ )}
X( )
1
-2
- -/2  /2 
X( /2)
2
-2
X( ) 1/2

-/2

-

X( + )
/2

X( /2+ )
1/2
2
A Filter Bank is a set of filters, linked by sampling
operators and delays.
L
x(n)
H
2
2
2
2
F

=x(n-1)
G
Without ( 2) and ( 2) and delay, PERFECT
RECONSTRUCTION would mean:
F0 (z)H0 (z) + F1(z)H 1(z) = I (the identity matrix)
Without ( 2) and ( 2) – with delay: l
F0(z)H 0(z) + F1 (z)H1(z) =I z
-l
We expect an overall delay z -l as every filter is
causal.
Let’s follow the signal through the filter:
Y( ) = ½ {X( /2) + X( /2 + )} = V()
After upsampling:
U( ) =V(2 ) =½ {X( ) + X( + )}
In the z domain:
V( z ) = ½[X(z1/2) + X(-z1/2)]
After upsampling:
U(z)= ½[X(z) + X(-z) ]
So, at the output of the Lowpass filter we have:
½[H0(z)X(z) + H0(-z)X(-z)]
-the „-“ sign comes from:  +  which = -z
After filtering by the lowpass filter of the synthesis
we have:
½F0(z)[H 0(z)X(z) + H0 (-z)X(-z)]
For the highpass part we get a similar rel.:
½F1 (z)[H 1(z)X(z) + H 1(-z)X(-z)]
We find the output by adding up these two terms.
F0 ( z) H 0 ( z)  F1 ( z) H1 ( z)X ( z) 
1
2 F0 ( z ) H 0 ( z )  F1 ( z ) H1 ( z )X ( z ).
1
2
- For perfect reconstruction with a time delay l,
the output must be: z l X (z ).
l
-that is, the “distortion term” must be z ,
the alias term must be zero.
The sampling operators introduced ALIASING:
These are the last terms.
ALIAS CANCELLATION CONDITION:
½F0 (z) H 0(-z)X(-z) + ½F1 (z) H 1(-z)X(-z) =0
From where it results:
F0 (z) H 0(-z) + F1 (z) H 1(-z) =0
1.
If the alias cancellation eq. is fulfilled, then
at the output we have:
½F0 (z)H0(z)X(z) + ½F1 (z)H1(z)X(z)
To get PR, the output has to be: x(n-l) which in the
-l
z is: X(z) = z X(z)
It results:
½[F0 (z)H0 (z)+F1 (z)H 1(z)]X(z) =z -l X(z)
Finally, the NO DISTORTION CONDITION:
F0 (z)H0 (z)+F1 (z)H 1(z) =2z -l
2.
The „extra 2“ is the result of that, that we built up
the FB from 2 filters.
These two conditions, in matrix form:
H0 (z)
F0 (z)
F1 (z)
H0 (-z)
=
H 1(z)
2z -l
0
H1 (-z)
The question is: how to design filters that meet
these conditions?
There are many ways!
One of them: to determine some of the filters
from the others!
From the Alias Cancellation cond. we choose:
F0 (z) = H1 (-z) and F1 (z) = - H 0(-z)
Define the product filter: P0 (z) = F0(z) H0(z)
and
P1 (z) = F1 (z) H 1(z)
P0 (z) it is: lowpass filter product,
P1 (z)
highpass filter product
- from the alias cancellation condition results:
P1 ( z)  H0 ( z) H1 ( z)  H0 ( z) F0 ( z)  P0 ( z)
Then, the No Distortion cond.:
F0 (z)H0 (z)+F1 (z)H 1(z) = F0 (z)H0 (z)-F 0(z)H0(-z) =
l
 P0 ( z)  P0 (z)  2z .
(cond. on odd powers)!!!
P0 (z) – P0 (-z) =2z -l
2.
Example:
H0
a, b, c
H 1 p, q, r, s, t
p, -q, r, -s, t
F0 (z) = H 1(-z)
-a, b, -c
F1 (z) = - H 0(-z)
STEPS: - Design a lowpass prod. Filter
which satisfies eq. 2.
- Factor P0 into F0 and H0 .
Then find F1 and H1.
F0 (z) H 0(-z) + F1 (z) H 1(-z) =0
1.
F0 (z)H0 (z)+F1 (z)H 1(z) =2z -l
2.
P0 (z) = F0 (z) H 0(z)
P1 (z) = F1 (z) H1 (z)
F0 (z) = H1 (-z) and F1 (z) = - H 0(-z)
P0 (z) – P0 (-z) =2z -l
2. – (4.9)
- There are many ways to design P0 in step 1.
and there are many ways to factor it in step 2.
-The length of P0 determines the sum of the
length of F0 and H 0.
IMPORTANT: EQ. 2. IS A CONDITION ON
THE ODD POWERS!
Those odd powers MUST HAVE coefficient zero
except z -l .
Which has coefficient one.
Eq. 2. can be made a little more convenient.
The left side is an odd function, so l is odd!
Normalize P0 (z) by z l to center it:
l
The normalized product filter is P(z) = z P0 (z).
l
P(-z)=(-z) P0 (-z) = -zl P0 (-z)
THE PERFECT RECONSTRUCTION
CONDITION THEN:
P(z) + P(-z) =2
Or:
3.
P() + P(+ ) =2 4.
- cond. on the even power
FILTER BANKS – PERFECT RECONSTRUCTION
From 3. and 4. it results that:
P(z) must be „halfband filter“!
That is, all even powers in P(z) are zero,
except the constant term, z0 which is 1!
The odd powers cancel, when
P(z) combines with P(-z).
In practice, many times we first determine:
H0 and H 1
Ex.:
Order flip
H0
a, b, c, d
Alternating flip
H 1 d, -c, b, -a
d,c,b,a
F0 (z) = H 1(-z)
Alternating sign
-a, b, -c, d
F1 (z) = - H 0(-z)
Example:
1 4
P0 ( z)  (1  z ) Q( z)
We have to choose the polynomial Q (z )
in such a way to satisfy PR condition.
That is, the odd powers must have coefficient
l
zero, except z has coefficient one.
A possible choice for
Q(z) :
1
Q( z)  1  4z  z
2
Results:
1 4
1
2
P0 ( z )  (1  z ) (1  4 z  z )
1
2
3
4
6
P 0 ( z )  (1  9 z  16z  9 z  z ).
16
- it results:
P0 ( z)  P0 ( z)  2z
3
The central term is:
z
3
z
l
The centering operation gives P(z):
P( z)  z P0 ( z)
3
1
3
1
3
P( z )  ( z  9 z  16  9 z  z ).
16
This is halfband: the only term with even
0
power is z with coefficient 1. So, the PR
condition is verified. The odd powers cancel
in the sum.
- that is:
P( z )  P (  z )  2
- in the frequency domain:
P( )  P(   )  2
- halfband condition;
P0 ( z)  P0 ( z)  2z
3
- condition on the odd powers;
P( z )  P (  z )  2
- condition on the even powers;
- How to factor P0 ( z ) ?
- There are a variety of factorizations in:
P0 ( z)  H 0 ( z) F0 ( z)
- the polinomial P0 ( z )– prewious ex., has 6 roots.
- the 2 roots from Q(z) are:
z  2 3
- the other 4 roots are at: z=-1.
- note:
2 3 
Im
-4th order
zero;
1
2 3
2 3
Re
1
z=-1
Ex.:
1
2
F0 (z)  1/ 4(1 2z  z ).
1
2
3
4
H0 (z)  1/ 4(1 2z  6z  2z  z )
1 2
(1  z )
Then H0 :
2 3
Im
-2nd order
zero;
2 3
-2nd
order
zero;
Re
z=-1
1
- filter length =3
¼ {1,2,1}
-1
1
- filter length =5
¼ {-1,2,6,2,-1}
- symmetric filters, (linear phase);
2 3
-2nd order
zero;
-1
2 3
-2nd order
zero;
-1
1
1
2 3
Orthogonal filters; (min. phase/max. phase)
- filter length =4
1
4 2
1 
- filter length =4
3 ,3  3 ,3  3 ,1  3
1
4 2
1 

3 ,3  3 ,3  3 ,1  3

- in this case one filter is the transpose of the other;
- in the time domain:
- in the z domain:
f 0 [n]  h0 [3  n]
3
1
F0 ( z)  z H0 ( z )
- the Haar filter bank: H 0 ( z ) 
1
2
H1 ( z ) 
1
2
(1  z 1 )
1
(1  z )
- from the alias cancellation cond.:
F0 ( z )  H1 ( z ) 
1
2
F1 ( z )   H 0 ( z ) 
(1  z 1 )
1
2
1
(1  z )
- the product filter:
P0 ( z)  F0 ( z) H0 ( z)  12 (1  z 1 )2
- the perfect reconstruction cond.:
P0 ( z )  P0 ( z ) 
1
2
1
2
1
2
(1  2 z  z )  (1  2 z  z ) 
2z
1
1
2
1
P( z)  z P0 ( z)  (1  z)(1  z )
1
1
2
Im
-2nd order
zero;
Re
z=-1
1
- zeros of P(z):
1 z  0
1  z 1  0
Orthonormal Filter Banks
C
x(n)
D
2
2
2
2
C
T
T
D

=x(n-1)
ORTHOGONALITY CONDITION or
CONDITION „O“
L
B
c(3) c(2) c(1) c(0)
0
0
0 c(3) c(2) c(1)
.....
=
d(3) d(2) d(1) d(0) 0
0
0
d(3) d(2) d(1)
.....
0 0.
c(0) 0 .
0
d(0) 0
LT BT
=
c(3) 0
c(2) 0
c(1) c(3)
c(0) c(2)
. c(1)
. c(0)
.
.
.
.
.
.
d(3)
d(2)
d(1)
d(0)
0
0
0
0
d(3)
d(2)
d(1)
d(0)
PR condition – for orthogonal filter banks.
L
B
T
L
T
B
LL = I :
T
LB = 0 :
T
BB = I :



0
0
I
=
It results:
T
I
c(n)c(n-k) = (k)
c(n)d(n-k) = 0
d(n)d(n-k) =  (k)
This condition above, is called Cond. O or
double shift orthogonality.
2
2
2
2
c(0) +c(1) + c(2) + c(3) = 1
and
c(0)c(2) + c(1)c(3) = 0
- Cond. O, imposes two constraints on four
coeff.
- For the high_pass coeff. similar relations
hold.
2
2
2
2
d(0) +d(1) + d(2) + d(3) = 1
and
d(0)d(2) + d(1)d(3) = 0
- From these cond. results:
THE FILTER LENGTH CANNOT
BE ODD!!! (Example!!!)
If the filter length would be odd (N+1),
then:
L
B
c(4) c(3) c(2) c(1) c(0) 0 0
0 0 c(4) c(3) c(2) c(1) c(0)
0 0 0
0 c(4) c(3) c(2)
0 0 0
0
0
0 c(4)
= .....
d(4) d(3) d(2) d(1) d(0) 0
0
0 d(3) d(2) d(1) d(0)
.....
0 0 0 0
0 0 0 0
c(1) c(0) 0 0
c(3) c(2) c(1) c(0)
0
0
0
0
For the orthogonality cond.:
c(4)*0+c(3)*0+c(2)*0+c(1)*0+c(0)*c(4)=0
which would mean: c(0)*c(4)=0 !!!
That is: N cannot be even, so the length of
the Filter cannot be odd!!!
The PR, or halfband cond., or Cond. O
For the filters in the frequency and z
Domain:
2
2
|C(z)| + |C(-z)| = 2
2
2



+
|C( )| + |C( )| =2
-where :
ze
i
2
2
|D(z)| + |D(-z)| = 2
2
2



|D( )| + |D( + )| =2
The highpass filter coeff. can be expressed
by the lowpass filter coeff.:
d(k)=(-1)k c(N-k)
Which is the alternating flip.
L
B
c(3) c(2) c(1)
0
0 c(3)
.....
=
-c(0) c(1) -c(2)
0
0
-c(0)
.....
c(0)
0
0 0.
c(2) c(1) c(0) 0 .
c(3) 0
0
c(1) -c(2) c(3) 0
- In reality because of time reversal there is a
change in sign: c(0), -c(1)...
The alt. flip cond. for the filters:
-N
-1
D(z)=(-z) C(-z ) where z=e
i
It results, that the main problem is, to design
the lowpass filter C!!!
It would be ideal to have:
- symmetric
- orthogonal filters
But: this is possible only in the Haar case!
Let us denote: P() = |C()|
2
Power Spectral Response
it is an autocorrelation function:
> 0 = C( ) C( )
That is, C() or C(z) is a spectral component of a
HALFBAND FILTER!
It is real and nonnegative.
A similar relation holds for P(1 ) =|D( )|
2
N
P( )   p(n)e
N
in
N
 (  c ( k )e
ik
N
)( c(k )e )
0
ik
0
This is a real and nonnegative filter: p(n)=p(-n)
Or:
p(n)  c(k )c(k  n)
That is, p(n) is the autocorrelation of the coeff. c(k).
That is, P is symmetric, it is the autocorrelation
filter.
-C is the start of an orthogonal filter bank if and
only if the autocorrelation filter, P is a
HALFBAND FILTER!
- in the z domain the halfband filter cond.:
P(z)+P(-z)=2
- from where it results:
 1 if
p(m)   c(k )c(k  2m)  
0 if
m0
m  0.
Example:
This is never negative!

P( ) = 1 + cos
The PR cond. is satisfied
Or, in the z domain:
z-1 + z
P(z) = 1+
2
-i
i
P( ) =|C( )| = (1+ e )(1+e ) /2
=1+cos
2
The requirement P( ) > 0 is crucial!
That is:
P( z )  P( z ) 
1  cos  1  cos  2
Is a halfband filter.
It is positive and =0 when
When C ( )  0 as well.
 
or z=-1;
- for C ( ) it results:
i
C()  (1  e ) / 2.
- the coefficients come from the familiar
averaging filter:
c(0)  c(1)  1/ 2.
Another example: the famous Daubechies
Orthogonal filter.
2
P( ) = ( 1+ cos  ) (1- ½ cos )=
(1  2 cos  cos2  )(1  1 2 cos )  1  3 2 cos  1 cos3 
2
Has a double zero at  = 
- If we would keep the first term only – the square
of the previous ex., - they would lead to the hat
function.
- This would not yield a halfband filter.
- In the z-domain we are squaring:
1
1  2 ( z  z)
2
which produces the even power: z or 2
1
term.
- This is not halfband!
Daubechies extra factor must be included to cancel
this term.
Q( )  1 
Q( z )  1 
1
1
2
cos
1
4
( z  z)
- Thanks to this factor, we got a halfband filter.
2
In the z-domain z is missing.
1
1
P ( z )  [1  1 2( z  z )] [1  1 4 ( z  z )] 
1  1 4 ( z
1
2
1

1
 z)  ( z  z) * 1 1 4 ( z  z)
2
In the z domain:
P(z) =1/16[ -z 3 + 9z1 +16 + 9z-1-z-3 ]
From where:
-1
-2
-3
C(z)=1/4 2[(1+ 3)+(3+ 3)z +(3- 3)z +(1- 3)z ]

-2nd order
zero;
-1
2 3
-2nd order
zero;
-1
1
1
2 3
Orthogonal filters; (min. phase/max. phase)
- filter length =4
1
4 2
1 
- filter length =4
3 ,3  3 ,3  3 ,1  3
1
4 2
1 

3 ,3  3 ,3  3 ,1  3

- These are the 4 coeff. of the famous Daubechies
filter D4.
1.- The halfband filter has P(z)+P(-z)=2.
The factor C(z) goes into an orthonormal filter
bank. D(z) comes from C(z) by an alternating flip.
2. – The response C(z) has a double zero at z=-1.
In the frequency domain C( ) has a double
Zero at =  .
The response is flat at  because of (1+cos ) 2.
2
1
|P( )|
P(/2) =1
|C( )|
Spectral Factorization:
- Two questions arise immediately:
- can every polynomial with P ( )  0
be factored into C ( ) ?
2
- how can this spectral factorization be done?
- The answer to the first question is: YES,
Féjer-Riesz Theorem;
- The answer to the second question is not so easy;- there are different methods;
- Method A: - zeros of a polynomial
- with real symmetric coeff. p(n), we have
P( z )  P(1 / z )
- If
zi is _ a _ root, _ so _ is 1/ zi .
-When z i is inside the unit circle, 1 zi is outside.
- The roots on the unit circle must have even
multiplicity.
- real coeff. ensure that the complex conjugate z
is a root when z is a root;
- the complex roots of the unit circle come
4 at a time;
- to get real filters, z and
z
must stay together;
- to get orthogonal filters z and 1 zi
separately:
must go
- this is the spectral factorization C(z)C( 1 zi );
MAXFLAT (DAUBECHIES) FILTERS
- these filters (and wavelets) are orthogonal
- the frequency responses have maximum flatness
at  = 0 and  = 
- the highpass coefficients come from the lowpass
by alternating flip.
Condition „O“ for a maxflat filter:
P  C ( )
2
is a normalized halfband filter;
- that is: p(0)=1 and p(2)=p(4)=...=0
Condition Ap :
C() has a zero of order p at:
 ,
that is:
( p 1)



C( )  C ( )  C ( )  ...  C
( )  0.
- The eq.: C ( )  0 says that
n
c
(
n
)(

1
)
 0.

- That is, the odd-numbered coefficients have
the same sum as the even-numbered coeff.:
Cond. A1 on c(n):
 c(n)   c(n).
odd _ n
even_ n
- Cond. Ap on c(n):
2 p 1
 (1)
n
n c ( n)  0
k
n 0
for k=0,1,...,p-1.
k
- factor n comes from the k-th derivative of
 c(n)e
n
- (1) comes from substituting
in
  .
- Because: P ( )  C ( ) , P also vanishes
when    .
2
.
The odd sum must be 1, since the only nonzero
even-numbered coeff. is p(0)=1:

p(n) 
odd _ n
 p(n)  p(0)  1.
even_ n
The sum over all n is P(0)=2. This yields:
C (0)  2
so the DC term at
is not reversed in sign.
0
The p zeros at  mean that C ( ) has a
i p
factor (1  e
)
Condition Ap on
C ( ) :
 i
p
1 e 
 R( ).
C ( )  
 2 
P
(

)
From where it results,
that
has
a
p
factor:  1  cos  .

2

Example:

p
1+cos

P( )=2(
)
2
k
k
p+k-1
1-cos

(
)(
).
k
2
How to design Maxflat Filters:
-The coeff. of the H_P Filter come from
the L_P Filter coeff. by Alternating-Flip.
The degree of the product Filter P is:
2N=4p-2
1 z
P0 ( z )  
 2
1
2p

 Q2 p  2 ( z )

- Let us take the binomial series:
truncated after p terms:
(1  y)
p
 2 p  2  p 1
p( p  1) 2
 y
B p ( y )  1  py 
y  ...  
2
 p 1 
 (1  y )
p
 O( y ).
p
- the remainder O, has order
first term to be dropped.
y
p
because this is the
-we combine B(y) with the factor: (1  y)
that has p zeros at y=1.
- the variable y on [0,1] will correspond to the
frequency  on [0,  ].
p
~
p
p
p
p
P ( y )  2(1  y ) B p ( y )  2(1  y ) [(1  y )  O( y )]
 2  O( y ).
p
- this product has exactly the flatness we want at
y=0.
- It is a polinomial of degree 2p-1, with 2p coeff.,
and satisfies p cond. at each endpoint:
- it’s first p-1 derivatives are zero at
~
y=0 and y=1, except P(0)  2.
- if we go from ordinary polinomials in y to
trigonometric polinomials in  the change
from 0  y  1 into 0     is:
1  cos

y


2

1  cos
1  y 
2

Then, we get:
 1  cos 
P( )  2

2


 p  k  1 1  cos 


 .

k
2

k 0 

p p 1
k
1
zz
- In the z-domain: 1  2 y  cos 
.
2
p
k
1 p p 1
1 k
 p  k  1 1  z   1  z 
1 z  1 z 
  
 .

P( z )  2
 
 
 2   2  k 0  k
 2   2 
Some columns (163,...) of S: Some basis vectors: wavelets
2
0
-2
0
50
100
150
200
250
300
0
50
100
150
200
250
300
0
50
100
150
200
250
300
0
50
100
150
200
250
300
2
0
-2
2
0
-2
2
0
-2
Rows 2,4,6,129,131,133 of the Analysis Matrix: Lev 1
2
0
-2
0
2
2
4
6
8
10
12
14
16
18
20
2
4
6
8
10
12
14
16
18
20
2
4
6
8
10
12
14
16
18
20
2
4
6
8
10
12
14
16
18
20
2
4
6
8
10
12
14
16
18
20
2
4
6
8
10
12
14
16
18
20
0
-2
0
2
0
-2
20
0
-2
0
2
0
-2
0
2
0
-2
0
Daubechies – orthogonal wavelets
Wrap around – 3. (row:127);4. (row:128)
Biorthogonal PR:
Theorem: In a biorthogonal linear-phase Filter
Bank with two channels, the filter lengths are
all odd or all even! The analysis filter can be
a) – both symmetric of odd length;
b) - one symmetric and the other antisymmetric
of even length.
How to design orthogonal or biorthogonal
filters?
- to have real valued filters: z and z MUST
stay together;
-1
- for biorthogonal filters: z and z MUST
stay together;
-1
- for orthogonal filter: z and z MUST go into
different filters; this gives C(z)C(z )
-for an orthogonal FB, the filter length cannot
be odd!!!
The space L2
Definition: For an interval a  t  b,
2
the space L ([a, b]) is the set of all square
integrable functions defined on a  t  b.
In other words,


2
L ([a, b])   f : [a, b]  C;  f (t ) dt  .
a


b
2
-funcions that are discontinuous are allowed as
members of this space;
- all the ex. used are or cont., or discontinuous
at a finite set of points;
b
- the condition:

f (t ) dt  
2
a
physically means that the total energy of the
signal is finite;
2
- the space L ([a, b]) is infinite dimensional;
- ex.:if a=0 and b=1, then the set of functions
{1, t, t2, t3, ...}is linearly independent and
2
belongs to L ([0,1])
- the function f(t)=1/t is an ex. of a function
2
that does not belong to L ([0,1]) since:
1
(
1
/
t
)
dt


.

2
0
L2 Inner Product:
Definition: The
is defined as
f ,g
b
2
L
L2
inner product on L
  f (t ) g (t )dt for
a
2
([a, b])
f , g  L ([a, b]).
2
- remember: for two vectors, X=(x1,x2,x3), and
Y=(y1,y2,y3) in R3, the standard (Euclidian)
inner product of X and Y is defined as
X , Y  x1 y1 x2 y2  x3 y3.
Multiresolution Analysis
- recall that: the sampling theorem approximately
reconstructs a signal f from samples taken
uniformly at intervals of length T.
- if the signal is band limited and its Nyquist
frequency is less than 1/T , then the reconstruction
is perfect;
- the smaller T is, the better we can approximate
or resolve the signal;
- the size of T measures our resolution of the
signal f.
- A typical FFT analysis of the samples taken
from f works at one resolution T.
- if the signal has bursts where it varies rapidly
with periods where it is slowly varying, this
single resolution analysis does not work well;
- to solve this problem:
- replace the space of band limited functions
with one tailored to the signal;
- analyze the signal using the scaled versions
2
of the same space, at resolutions T/2, T / 2 and
so on;
- Hence the term multiresolution analysis.
DEFINITION: Let V j , j= ..., -2,-1,0,1,2,3,...be
a sequence of subspaces of functions in L2 ( R).
The collection of subspaces V j , j  Z  is
called a multiresolution analysis with scaling
function  if the following conditions hold:
1.
2.
3.
4.
V j  V j 1
(nested)
(density)
Vj  L2 (R)
(separation) V j  0
(scaling) The function f(x) belongs to V j
j
if and only if the function f (2 x)belongs
to V0 .
5.
(orthonormal basis) The function 
belongs to V0 and the set  ( x  k ), k  Z 
is an orthonormal basis (using the L2
inner product) for V0 .
- there may be several choices of 
corresponding to a system of approximation
spaces;
- different choices for  may yield different
multiresolution analysis;
-  don’t have to be orthonormal;
- all that is needed is a  for which the set
 ( x  k ), k  Z is a basis.


-the most useful class of scaling functions are
those that have compact or finite support;
- the Haar scaling function is a good example
of a compactly supported function;
- the scaling functions associated with
Daubechie’s wavelets are not only compactly
supported, but also continuous;
- the two sharp spikes: noise
- the signal can be approximated using Haar
building blocks;
Haar Decomposition and Reconstruction
Algorithms
- Decomposition
Definition:
Suppose, j is any nonnegative integer. The space
of step functions at level j, denoted by V j is
defined to be the space spanned by the set:
 (2 x 1), (2 x), (2 x 1), (2 x  2)
j
j
over the real numbers.
j
j
V j is the space of piecewise constant
functions of finite support whose
discontinuities are contained in the set:
...,1/ 2 ,0,1/ 2 ,2 / 2 ,3 / 2 ,...
j
j
j
j
- a function in V0 is a piecewise constant
function with discontinuities contained in the
set of integers;
...,1/ 2 ,0,1/ 2 ,2 / 2 ,3 / 2 ,...
0
0
0
0
-any function in V0 is also contained in V1
which consists of piecewise const. fun. whose
discont. are contained in a set of half integers
...,1/ 2,0,1/ 2,1,3 / 2,...
- the same applies for
V1 and V2 : V1  V2 .
V0  V1  V2    V j 1 V j  V j 1    .
- this containment is strict;
- for ex. the fun.  (2 x) belongs to V1 but
does not belong to V0 - since  (2 x) is discont.
at x=1/2.
y
a2
a0
a4
a1
-1
0
1
2
a1
3
4
x
a3
Graph of typical element in V0 .
A fun. in V0 may not have disc. if for ex. a1
 a2
4
1
1/2
Graph of  (2 x)
-The graph of the
function:
f ( x)  4 (2 x)  2 (2 x  1) 
2 (2 x  2)   (2 x  3)
1
-1
2
3
- V j contains all relevant information up to a
j
resolution scale of order 2 .
- as j gets larger, the resolution gets finer;
- One way to decompose a signal efficiently is to
construct an orthonormal basis for V j .

- V0 is generated by
2
have unit norm in L ;
that is:

and its translates and
k 1
 ( x  k ) L    ( x  k ) dx   1dx  1
2
2
2

k
- if j is different from k then:
 ( x  j ), ( x  k )


L2
   ( x  j ) ( x  k )dx  0

- so the set  ( x  k ), k  Z
basis for V0 .
jk
 is an orthonormal
- the idea is to decompose V j as an orthogonal
sum of V j 1 and its complement;
- if j=1, the orthogonal complement of V0 in
V1 is generated by the translates of some
function ;
1. -  is a member of V1 and so  can be
expressed as:  ( x)  l al (2 x  l ) for some
choice of al  R.
(And only a finite
number of the al are nonzero.)
2. -  is orthogonal to V0. This is equivalent
to  ( x) ( x  k )dx  0 for all integers k.
- the first requirement means that  is built
from blocks of width ½.
- the second requirement with k=0 implies:


(
x
)

(
x
)
dx

0
.


- the smallest  satisfying both of these
requirements is the function whose graph is:
- This graph consists of two
blocks of width one-half and
can be written as:
1
1
1/2
 ( x)   (2 x)   (2 x  1)
-1

- satisfying the first requirement.
- in addition:
12
1

(
x
)

(
x
)
dx

1
dx

1
dx

0
.




0
12

- thus,  is orthogonal to .
- If k  0, then the support of  (x ) and the
support of  ( x  k ) do not overlap, so:

(
x
)

(
x

k
)
dx

0
.

- therefore,  belongs to V1 and is orthogonal toV0 .
MULTIRESOLUTION:
- wavelets produce a natural multiresolution of every image;
- coarser and coarser approx. of
func. f(t);
- averaging for bigger and bigger int.
THE SMOOTH SIGNAL at one level +
DETILS at the same level:
COMBINE INTO A MULTIRESOLUTION at the next FINER LEVEL
coarse aver. at lev. j + details at lev. j
(wave. coef.)
= signal at the finer level j-1

(x)

)
(x
a
 j 1k j 1k
k
)
(x
w
(x)
b

(x)

)
(x
a
 jk jk  jk jk
k
k
Where:
a jk   f ( x) jk (x)dx

b jk 


f (t )w jk (t )dt
V2
W2
x=a3k3k(t) =a2k2k(t) + b2kw2k(t) =
=a1k1k(t) + b1k w1k(t)+ b2kw2k(t) =
= a0k0k(t) + b0kw0k (t)+ b1kw1k(t) +b2kw2k(t)
V0
W0
W1
W2
Some Applications:
- ECG Analysis (Laszlo, Schipp);
- Construction of Haar-like systems
(Laszlo, Schipp);
- Storage, compression, reconstruction...
(László);
- Fusion (Kozaitis, László) etc.
- Edge detection from a noisy image
(Kozaitis)
h1(x) = (-1)n h0(N - n),
f0(n) = h0(N - n), Orthogonal case
f1(n) = (-1)n+1 h0(n),
Rows 2,4,6,129,131,133 of the Analysis Matrix: Lev 1
2
0
-2
0
2
2
4
6
8
10
12
14
16
18
20
2
4
6
8
10
12
14
16
18
20
2
4
6
8
10
12
14
16
18
20
2
4
6
8
10
12
14
16
18
20
2
4
6
8
10
12
14
16
18
20
2
4
6
8
10
12
14
16
18
20
0
-2
0
2
0
-2
20
0
-2
0
2
0
-2
0
2
0
-2
0
Some rows (101...) of A: Some basis vectors: wavelets
1
0
-1
0
50
100
150
200
250
300
0
50
100
150
200
250
300
0
50
100
150
200
250
300
0
50
100
150
200
250
300
1
0
-1
1
0
-1
1
0
-1
Some columns (163,...) of S: Some basis vectors: wavelets
2
0
-2
0
50
100
150
200
250
300
0
50
100
150
200
250
300
0
50
100
150
200
250
300
0
50
100
150
200
250
300
2
0
-2
2
0
-2
2
0
-2
Part1: Level 4 WT
Part1: Level 3 WT
50
50
100
100
150
150
200
200
250
250
50
100 150 200 250
50
Part1: Level 2 WT
Part1: Level 1 WT
50
50
100
100
150
150
200
200
250
250
50
100 150 200 250
100 150 200 250
50
100 150 200 250
input image
horizontal edge
vertical edge
WT
WT
bh1
b1v
Wch
Wcv
bh2
b2v
Wh
threshold
Wv
bh3
b3v
edge location filter
X
bh4
X
b4v
Mv
M
Mh
modulus
edge result
inverted result
1. HW
Design a Filter Bank (FB) with two channels,
Lowpass (LP) and Highpass (HP) using the Haar
coefficients: - for LP: h(0)=1/2 and h(1)=1/2
but because of downsampling you put in: 2/2!
- for HP: h(0)=1/2 and h(1)= - 2/2
Build up 8 levels!!!
Show:
- by putting on the same figure 5 neighbouring
rows, show the wavelets at least at 3 different
Levels! Mandatory: the finest level and the
coarsest level
- the first (10-12) rows and columns of the (colored
picture) analysis matrix A=A8*A7*...A1
- represent (colored picture) the coefficients;
show that the lowpass and highpass coefficients are
at the right place –at least for two different levels
- put in a signal in the time domain (for ex. at the
finest level) and show that your FB performs perfect
reconstruction; (first one signal at the finest level, then
2 and 3 at the same level; then show this at least at one
more level.)
- show that the wavelets and the scaling functions as
well are the same in the synthesis matrix as in the
analysis matrix;
Extra: - put in 3 signals into the wavelet domain at
3 consecutive places
(this is done in the same way as when you put in in the
time domain, the difference is that then you will
multiply by the synthesis matrix and not by the
analysis matrix)
and show the effect of the downsampling;
-once at the finest level, then on another two levels;
clc
clear
x = zeros(size(1:256));
x(150) = 1;
x = x';
[row col]=size(x);
r=1/sqrt(2);
S=A';
% synthesis matrix
figure
subplot(5,1,1)
plot(A(1,:),'m');
figure
imagesc(A6(1:10,1:10))
colorbar
16
32
64
Even coarser
Coarser level
128
Finest scale
256
A
signal
HW_2.
Maxflat Filter design:
Write a program to compute and plot the
magnitude of the frequency response for the
maxflat filters for p= 3, 5, 7, 9, 11.
Verify the filters are halfband (take the IFT).
 1  cos 
H (e )  

2


i
 p  1  k  1  cos 


.

k 
2

k 0 
p p 1
k
HW_2: Part II. Orthogonal Filters
Given eq.:
C ( z) 

1
4
1
4 2
1  z  1  3  1  3 z 
1 2
1


1  3  3  3 z  3  3 z  1  3 z .
2
1
2
3
You have the filter coeff. That produce a wavelet.
Graph both its impulse, and magnitude of the
frequency response.
Generate the orthogonal high-pass filter using
the alternating-flip method and plot its impulse
and magnitude of the frequency response as well
as the Power Spectral Response.
Help: IFFT_C=real(fftshift(ifft(C)))
3. HW
Design a biorthogonal Filter Bank (FB) with two
channels, Lowpass (LP) and Highpass (HP) using
the coefficients of the polinomial:
P=[3 0 -25 0 150 256 150 0 -25 0 3];
Except of the effect of the downsampling, show
all those, you had to in the first HW. (This part is
considered as the second homework.)
Then, draw a picture in black-and-white to see the
Horizontal, vertical and diagonal components in the
Wavelet domain.
Then, take a nice picture you‘d like and show its
Transforms in the wavelet domain
At each (4 )level!!!
Show at the reconstruction as well at each level.
So, this time you need to build up 4 levels only!!!
Show the reconstruction of the picture.
clear
clc
format long e
p=[-1 0 9 16 9 0 -1]
Roots = roots(p)
% Po(z) = Fo(z)*Ho(z)
% Filter bank 1:BIORTHOGONAL The roots c and 1/c
are in the same polinomial
% This gives a linear phase filter bank (bi-orthogonal).
% Roots are 3.7321, .2679, -1 in one filter and -1 -1 -1 in
the other
%a=[-1 3.7321 .2679];
b=[-1 -1 -1];
Ho_B=poly(b) % 1 3 3 1
Fo_B=real(deconv(p,Ho_B)) % -1
3
3
-1
% Filter bank 2: ORTHOGONAL The roots c and 1/c
are in different polynomials. This gives an
% orthogonal filter bank.
% Roots .2679, -1, -1 in one filter and -1 -1 3.7327 in the
other
a2=[-1 -1 .2679];
b2=[-1 -1 3.7327];
Ho_O = poly(b2) % 1 -1.7327 -6.4654 -3.7327
Fo_O=real(deconv(p,Ho_O)) % -1 -1.7327 -0.46764
0.2544
P=[3 0 -25 0 150 256 150 0 -25 0 3]; % the coefficients of
the polinomial
Roots_P = roots(P)
z=[Roots_P(1,1) Roots_P(2,1) Roots_P(3,1) Roots_P(4,1)
Roots_P(10,1)];
Fo_B=poly(z)
Ho_B =real(deconv(P,Fo_B))
% Create synthesis matrix for each level
S_Lev1 = pinv(A1);
inputimage_orig =
imread('c:\MATLAB6p1\Project\f3.bmp','bmp');
inputimage=double(inputimage_orig(140:395,170:425));
%for the box2: (140:395,120:375);
[r,c]=size(inputimage);
figure
imagesc(inputimage)
title('original image')
colormap(gray)
The name of the picture used, ex.: f3.bmp
Biorthogonal case
h1(n) = (-1)n f0(n),
f1(n) = (-1)n+1 h0(n)
Perfect Reconstruction cond.: - PR
Alias cancellation:
F0(z)H0(-z)+F1(z)H1(-z)=0
No distortion:
F0(z)H0(z)+F1(z)H1(z)=2z-l
Perfect Reconstruction cond.:
Where:
P(z)+ P(-z)=2
(PR cond. In the z domain)
z= eiω
Next slides: - ex. of scaling functions
and wavelets
j=1, k=0
j=1, k=5
j=2, k=0
j=2, k=5
Command: - sort
- coef. in the matrix
- ex.: using a compression 16:1 more
then 61.000 coef. could be set to zero – from the
total 65.436 !!!