# Wavelet Analysis The Continuous Wavelet Transform

Document Sample

```					         Wavelet Analysis
The Continuous Wavelet Transform

The equation for a general probing function fm(t) that slides
along the signal (x(t)) is:

X(t, m) =  x( ) f(t -  )m d


where f(t)m is a function family with m specifying the family number.

In the Continuous Wavelet Transform (CWT) the probing function
has a finite time duration and its family members vary in scale:
        1      t - b
W(a,b) =  x(t)          *       dt

a     a 

where b acts to translate the function across x(t) just as t does in
the Short Term Fourier Transform, and the variable a acts to vary
the time scale of the probing function, ψ.
Wavelet Families
1                                                  1
a = 0.5                                            a=1
0.5                                                0.5

Wavelet

Wavelet
0                                                  0

-0.5                                               -0.5

-1                                                 -1

The Wavelet                                        -1      -0.5       0        0.5   1                -1     -0.5       0        0.5   1
shown is the                                                      Time (sec)                                        Time (sec)

‘Morlet Wavelet’
described by the                                  1                                                  1
a=2                                                a=4
equation:                                       0.5                                                0.5

                
Wavelet

Wavelet
 (t) = e - t2 cos 
0                                                  0
2
ln 2
t                 -0.5                                               -0.5

-1                                                 -1

-1      -0.5       0        0.5   1                -1     -0.5       0        0.5   1
Time (sec)
Time (sec)

• If b = 0, and a = 1, then the Wavelet is in its natural form,
the "mother wavelet.” If a > 1, the wavelet is stretched
along the time axis and if a < 1, the wavelet is contracted.

• There is a built-in tradeoff between time and frequency
resolution which is key to the CWT and makes it well
suited to analyzing signals with rapidly varying high-
frequency components superimposed on slowly varying
low frequency components.
Δωψ(a) Δtψ(a) = Δωψ Δtψ = constant

• This trade-off is illustrated in the next figure which shows
the time-frequency boundaries for the Mexican hat
wavelet for various values of a.
20

18
Time-frequency boundaries
16
of the Mexican Hat Wavelet
for various values of a. The
14                       area of each of these boxes

is constant.
12

10

8

6

4

2
0.2   0.4   0.6   0.8       1        1.2   1.4   1.6   1.8
Time (sec)
Example 7-1 Write a program to construct the CWT of a signal
consisting of two sequential sine waves of 10 and 40 Hz. (i.e. the signal
shown in Figure 6-1). Plot the Wavelet Coefficients as a function of a
and b. Use the Morlet Wavelet.
decr_a = .5;                                      % Decrement for a
a_init = 4;                                       % Initial a
wo = pi * sqrt(2/log2(2));                        % Wavelet freq. scale factor
ti = ((1:N1/2)/fs)*10;                      % Time vector for Wavelet
%
% Calculate Continuous Wavelet Transform
for i = 1:resol_level
a(i) = a_init/(1+i*decr_a);                   % Set scale
t = abs(ti/a(i));                             % Scale t vector for Wavelet
mor = (exp(-t.^2).* cos(wo*t))/sqrt(a(i));        % Morlet Wavelet
Wavelet = [fliplr(mor) mor];                  % Make symmetrical
ip = conv(x,Wavelet);                         % Convolve Wavelet & signal
ex = fix((length(ip)-N)/2);                   % Calculate extra points /2
CW_Trans(:,i) =ip(ex+1:N+ex,1);               % Remove extra points
end

This program generates the following graph.
Example 7-1 Figure
Continuous Wavelet Transform of
a signal that changes in step-like
manner between 10 and 40 Hz
Discrete Wavelet Transform
• The CWT has one serious problem: it is highly
redundant. The CWT provides an oversampling of the
original waveform: many more coefficients are generated
than are actually needed uniquely to specify the signal.

• If the application calls for recovery of the original signal,
all of the coefficients will be required and the
computational effort could be excessive.

• In applications that require bilateral transformations, a
transform that produces the minimum number of
coefficients required is preferred.

• The “Discrete Wavelet Transform” (DWT) achieves this
parsimony by restricting the variation in translation and
scale, usually to powers of 2.
Discrete Wavelet Transform (continued)

• In the DWT, a new concept is introduced termed the
"Scaling function," a function that facilitates
computation of the DWT. To implement the DWT
efficiently, the finest resolution is computed first. The
computation then proceeds to courser resolutions, but
rather than start over on the original waveform, the
computation uses a smoothed version of the fine
resolution waveform. This smoothed version is
obtained with the help of the Scaling function.

 (t) =    
n = -
2 c(n)  (2t - n)
Discrete Wavelet Transform (continued)

• In the DWT, the Wavelet itself can be defined from the
Scaling function:

 (t) =    
n = -
2 d(n)  (2t - n)

where d(n) is a series of scalars that are related to
the waveform x(t) and that define the discrete Wavelet
in terms of the Scaling function. While the DWT can
be implemented using the above equations, it is
usually implemented using filter bank techniques.
Filter Banks

x(n)
Ho()         ylp(n)   Simple filter bank
consisting of only two
yhp(n)   filters applied to the
H1()                  same waveform. The
filters have highpass and
H1()    lowpass spectral
|H()|                          characteristics. Filter
H0()    outputs consist of two
subbands: a lowpass
subband ylp(n) and a
Frequency (Hz) fs/2   highpass subband yhp(n).
Filter Banks (Continued)

Analysis                Synthesis
Filters                 Filters
x(n)              ylp(n)                            x’(n)
H0()                   G0()        
yhp(n)
H 1()                  G1)

• In some analysis applications only the subband signals are of
interest, and reconstruction is not needed, but in many Wavelet
applications, some operation is performed on the subband signals,
ylp(n) and yhp(n), before reconstruction of the output signal as
shown.
Filter Banks ~ Downsampling

• The Filter Bank approach would normally require extra
points to represent x(n), but if the Analysis filters are
correctly chosen, it is possible to limit the total filtered
points to be the same as the original data. Filtered signal
samples are reduced by one-half by eliminating every
other point.   (This operation is known as “downsampling” and is
illustrated schematically by the symbol ↓2.)

•    After the data are operated on, they are “upsampled”
(indicated by the symbol ↑2) to replace the missing
points with zeros. Subsequent filtering will replace the
zeros with correct data.
Example of a Wavelet application that uses three Filter
Banks and includes the downsampling and upsampling
operations. Downsampled amplitudes are sometimes scaled by /2, a
normalization that can simplify the filter calculations when matrix methods
are used.
Analysis
Filters            H 0()             
2       ylp (n)
x(n)
H 0()          
2
H 2()

2        yhp2(n)

H 1()          
2        yhp1 (n)

Synthesis
Filters
ylp(n)
x’(n)
2      G 0()                     2        G 0()        

2      G 1()                      2         G 1()
yhp2(n)                             yhp1(n)
Filter Bank Filter Design
• Filters must meet a number of criteria to be able to
recover the original signal after passing through the two
Filter Banks. Recovery is complicated by the
downsampling.
• Assuming only a high and lowpass filter, the criterion for
aliasing cancellation is (Strang and Nguyen, 1997):

G0(z) H0(-z) + G1(z)H1(-z) = 0

where H0(z) is the Transfer Function of the Analysis lowpass filter,
H1(z) is the Transfer Function of the Analysis highpass filter, G0(z) is
the Transfer Function of the Synthesis lowpass filter, and G1(z) is
the Transfer Function of the Synthesis highpass filter.
Filter Bank Filter Design (continued)

• The necessity to be able to recover the original
waveform from the subband waveforms places another
important requirement on the filters which is satisfied
when:

G0(z)H0(z) + G1(z)H1(z) = 2z-N

where the Transfer Functions are the same as in the previous
equation. N is the number of filter coefficients (i.e., the filter order);
hence z-N is just the delay of the filter.
Daubechies Filters
• In many analyses, it is desirable to have subband signals
that are orthogonal, placing added constraints on the
filters. Although no filter exists with all the desirable
properties, a number of filters have been developed that
have most of the desirable properties. A filters family
developed by Daubechies (and named after her) is
popular and has many of the desired features.
• The coefficients of the lowpass filter, h0(n), for the 4-
coefficient Daubechies filter are given as:

h(n) = [(1 + /3), (3 + /3), (3 - /3), (1 - /3)]/8

• The other filters (highpass and synthesis) can be
constructed from the lowpass filter given the constraints.
Analysis Filters

• To generate a highpass filter whose output is orthogonal
to that of the lowpass filter, the highpass filter frequency
characteristics must have a specific relationship to those
of the lowpass filter:
H1(z) = -z-NH0(-z-1)

• This criterion can be implemented by applying the
“alternating flip” algorithm to the coefficients of h0(n):
h1(n) = [h0(N), -h0(N-1), h0(N-2), -h0(N-3), ...]

• where N is the number of coefficients in h0(n). This algorithm is
easy to implement in MATLAB.
Synthesis Filters

• Synthesis Filters are fairly constrained by anti-alaising
and orthogonality requirements. The conditions can be
met by making G0(z) = H1(-z) and G1(z) = -H0(-z).
Synthesis Filter Transfer Functions are related to the
Analysis Transfer Functions by:
G0(z) = H1(z) = z-NH0(z-1)
G1(z) = -H0(-z) = z-NH1(z-1)
• Which can be implemented using the “order flip”
algorithm:
g0(n) = [h0(N), h0(N-1), h0(N-2), ....]
g1(n) = [h1(N), h1(N-1), h1(N-2), ...]
Example 7-3 Construct an Analysis Filter Bank containing L
decompositions; that is, a lowpass filter and L highpass filters.
Decompose and recover the signal using appropriate Filter Banks.
The Analysis Filter Bank is developed in the routine ‘analyze.’
•     function an = analyze(x,h0,L)
•     lf = length(h0);                          % Filter length
•     lx = length(x);                           % Data length
•     an = x;                                   % Initialize output
•     % Calculate High pass coefficients from low pass coefficients
•     for i = 0:(lf-1)
•        h1(i+1) = (-1)^i * h0(lf-i);           % Alternating flip algorithm
•     end
•     %
•     % Calculate filter outputs for all levels
•     for i = 1:L
•        a_ext = an;
•        lpf = conv(a_ext,h0);                  % Low pass FIR filter
•        hpf = conv(a_ext,h1);                  % High pass FIR filter
•        lpf = lpf(1:lx);                       % Remove extra points
•        hpf = hpf(1:lx);
•        lpfd = lpf(1:2:end);                   % Downsample
•        hpfd = hpf(1:2:end);
•        an(1:lx) = [lpfd hpfd];                % Low pass output at beginning of array, but now
occupies only                                          % half the data points as last pass
•        lx = lx/2;                             % Decrease array size by two for next level
•        ….. Then plot……
Example 7-3 (continued) The Synthesis Filter Bank is implemented in
the routine ‘synthesize.’

•     function y = synthesize(a,h0,L)
•     lf = length(h0);                           % Filter length
•     lx = length(a);                            % Data length
•     lseg = lx/(2^L);                           % Length of first low and highpass segments
y = a;                                     % Initialize output
•     g0 = h0(lf:-1:1);                           Lowpass coefficients using order flip,
•     %
•     % Calculate High pass coefficients, h1(n), from low pass coefficients
•     for i = 0:(lf-1)
•        h1(i+1) = (-1)^i * h0(lf-i);            % Alternating flip Eq. 7-18
•     end
•     g1 = h1(lf:-1:1);                          % Highpass filter coefficients using order
•     % Calculate filter outputs for all levels
•     for i = 1:L
•        lpx = y(1:lseg);                                     % Get Low pass segment
•        hpx = y(lseg+1:2*lseg);                              % Get High pass outputs
•        up_lpx = zeros(1,2*lseg);               % Initialize for upsampling
•        up_lpx(1:2:2*lseg) = lpx;               % Up sample low pass (every odd point)
•        up_hpx = zeros(1,2*lseg);               % Repeat for high pass
•        up_hpx(1:2:2*lseg) = hpx;
–     syn = conv(up_lpx,g0) + conv(up_hpx,g1);       % Filter and combine
•        y(1:2*lseg) = syn(1:(2*lseg));          % Remove extra points from end
•        lseg = lseg * 2;                        % Double segment lengths for next pass
•     end
Example 7-3 (continued) After passing through the
Analysis and Synthesis Filter Banks the original signal
(upper trace) shows little change (lower trace) except
for an initial artifact and a delay induced by the filters.

0

-2

-4

-6

-8

-10

-12
-0.2   0   0.2   0.4        0.6   0.8   1   1.2
Time (sec)
Low Pass Outputs                                   High Pass Outputs
5                                                 1

0                                                 0
Example 7-
3 (cont)         -5                                                -1
0          200            400          600        0               200            400          600
10                                                 1
These plots
show the         0                                                 0
output of
the various     -10
0          100            200          300
-1
0               100            200          300
filters. Only   10                                                 1
the lowest
0                                                 0
lowpass
subband is      -10                                                -1
0          50             100          150        0               50             100          150
used in         20                                                 1
synthesis.
0                                                 0
Filter
frequency       -20                                                -1
0     20           40           60     80         0          20           40           60     80
plots are        2                                                 2
Low Pass Filter                             High Pass Filter
also shown.      1                                                 1

0                                                 0
0   100      200        300      400   500        0       100       200        300      400   500
Frequency (Hz.)                                       Frequency (Hz.)
Example 7-4 Decompose the previous signal using a four-level Filter
Bank. Set all data points in the two highest resolutions segments that
are below some threshold to zero, then reconstruction the signal with
the Synthesis filters. Modify the Analysis and Synthesis filter routines to
use circular convolution to reduce the time delay.

an = analyze1(x,h0,4);            % Decompose the signal: Analysis Filters
% Set the threshold
threshold = 5*var(an(N/4:N));     % Set threshold based on variance
for i = (N/4:N)                   % Examine the two highest resolution
% highpass subbands
if an(i) < threshold
an(i) = 0;
end
end
sy = synthesize1(an,h0,4);        % Reconstruct original signal
Example 7-5 Add a small step change to the second derivative of a
waveform. Apply a three-level Analysis Filter Bank and examine the high
frequency subbands for evidence of this subtle discontinuity.
Low Pass Outputs                      x 10
-3               High Pass Outputs
2.5                                                                            4                                               5

2
2
1.5
0
0
1

-2                                              -5
0.5                                                                                 0          200        400         600           0               200         400           600
0
4                                            0.01

-0.5                                                                           2                                               0
-1
0                                            -0.01
-1.5

-2                                           -0.02
-2                                                                                 0          100        200         300           0               100         200           300
5                                            0.02
-2.5
0     0.1   0.2   0.3   0.4      0.5       0.6   0.7   0.8   0.9   1
Time (sec)
0
0
-0.02
The discontinuity at 0.5
-5                                           -0.04
sec. is not visible in the                                               0          50        100
Low Pass Filter         150           0          50
High Pass Filter        100           150
waveform, but is clearly                                           1.5                                            1.5

identified in the highpass                                          1                                               1

subbands after                                                     0.5                                            0.5

application of the                                                  0                                               0
0   100      200   300      400   500           0         100     200     300     400     500
Analysis Filter Bank                                                               Frequency (Hz.)                                       Frequency (Hz.)
LP
Wavelet Tree
(Logarithmic Tree)        LP

Wavelet Packets                                     LP
HP

A Wavelet packet tree                                    HP
(lower tree) does a
complete decomposition
generating both lowpass                             HP

and highpass subbands at
LP
every level. The outputs of   Wavelet Packet Tree
this decomposition can be     (Balenced Tree)            LP

HP
used in feature detection.                          LP
LP

HP

HP

LP

LP

HP
HP
LP

HP

HP
Balanced Tree            10                                        5

Decomposition             0                                        0

The signal at the
-10                                       -5
upper left has been            0   50                100   150          0   50                100   150
Time (sec)                               Time (sec)
lowpass filtered 3      0.2                                        2

times. The signal in
0                                        0
the upper right has
been lowpass filtered   -0.2                                      -2
twice, then highpass           0   50
Time (sec)
100   150          0   50
Time (sec)
100   150
0.1                                      0.1
filtered. The rest of
the plots follow the      0                                        0
tree structure
sequentially.           -0.1                                     -0.1
0   50                100   150          0   50                100   150
0.2             Time (sec)               0.5             Time (sec)
The three different
sinusoids contained
0                                        0
in the original
waveform can be         -0.2                                     -0.5
0   50                100   150          0   50                100   150
identified.
Time (sec)                               Time (sec)

```
DOCUMENT INFO
Shared By:
Categories:
Tags:
Stats:
 views: 44 posted: 7/27/2012 language: English pages: 26
How are you planning on using Docstoc?