The Story of Wavelets: Theory and Engineering Applications
Download
Report
Transcript The Story of Wavelets: Theory and Engineering Applications
Presents
Robi Polikar
Dept. of Electrical & Computer Engineering
Rowan University
The Story of Wavelets
Technical Overview
But…We cannot do that with Fourier Transform….
Time - frequency representation and the STFT
Continuous wavelet transform
Multiresolution analysis and discrete wavelet transform (DWT)
Application Overview
Conventional Applications: Data compression, denoising, solution of PDEs,
biomedical signal analysis.
Unconventional applications
Yes…We can do that with wavelets too…
Historical Overview
1807 ~ 1940s: The reign of the Fourier Transform
1940s ~ 1970s: STFT and Subband Coding
1980s & 1990s: The Wavelet Transform and MRA
What is a Transform
and Why Do we Need One ?
Transform: A mathematical operation that takes a function or sequence
and maps it into another one
Transforms are good things because…
The transform of a function may give additional /hidden information
about the original function, which may not be available /obvious
otherwise
The transform of an equation may be easier to solve than the original
equation (recall your fond memories of Laplace transforms in DFQs)
The transform of a function/sequence may require less storage, hence
provide data compression / reduction
An operation may be easier to apply on the transformed function, rather
than the original function (recall other fond memories on convolution).
December, 21, 1807
“An arbitrary function, continuous or with
discontinuities, defined in a finite interval by an
arbitrarily capricious graph can always be
expressed as a sum of sinusoids”
J.B.J. Fourier
Jean B. Joseph Fourier
(1768-1830)
Complex Function weight i Simple Function i
i
Complex function representation through simple building blocks
Basis functions
Using only a few blocks Compressed representation
Using sinusoids as building blocks Fourier transform
Frequency domain representation of the function
F ( ) f (t )e
jt
dt
1
jt
f (t ) F ( )e d
2
How Does FT Work Anyway?
Recall that FT uses complex exponentials (sinusoids) as building
blocks.
jt
e
cost j sint
For each frequency of complex exponential, the sinusoid at that
frequency is compared to the signal.
If the signal consists of that frequency, the correlation is high
large FT coefficients.
F ( ) f (t )e jt dt
f (t )
1
jt
F
(
)
e
d
2
If the signal does not have any spectral component at a frequency,
the correlation at that frequency is low / zero, small / zero FT
coefficient.
FT At Work
x1 (t ) cos(2 5 t )
x2 (t ) cos(2 25 t )
x3 (t ) cos(2 50 t )
FT At Work
x1 (t )
F
X 1 ( )
x2 (t )
F
X 2 ( )
x3 (t )
F
X 3 ( )
FT At Work
x4 (t ) cos(2 5 t )
cos(2 25 t )
cos(2 50 t )
x4 (t )
F
X 4 ( )
FT At Work
Complex exponentials
(sinusoids) as basis
functions:
F ( )
f (t ) e
jt
dt
1
f (t )
2
F
An ultrasonic A-scan using 1.5 MHz transducer, sampled at 10 MHz
F ( ) e
jt
dt
Stationary and Non-stationary
Signals
FT identifies all spectral components present in the signal, however it
does not provide any information regarding the temporal (time)
localization of these components. Why?
Stationary signals consist of spectral components that do not change in
time
all spectral components exist at all times
no need to know any time information
FT works well for stationary signals
However, non-stationary signals consists of time varying spectral
components
How do we find out which spectral component appears when?
FT only provides what spectral components exist , not where in time
they are located.
Need some other ways to determine time localization of spectral
components
Stationary and Non-stationary
Signals
Stationary signals’ spectral characteristics do not change with time
x4 (t ) cos(2 5 t )
cos(2 25 t )
cos(2 50 t )
Non-stationary signals have time varying spectra
x5 (t ) [ x1 x2 x3 ]
Concatenation
Stationary vs. Non-Stationary
X4(ω)
Perfect knowledge of what
frequencies exist, but no
information about where
X5(ω)
these frequencies are
located in time
Shortcomings of the FT
Sinusoids and exponentials
Stretch into infinity in time,
no time localization
Instantaneous in frequency,
perfect spectral localization
Global analysis does not allow analysis of non-stationary signals
Need a local analysis scheme for a time-frequency representation
(TFR) of nonstationary signals
Windowed F.T. or Short Time F.T. (STFT) : Segmenting the signal into
narrow time intervals, narrow enough to be considered stationary, and
then take the Fourier transform of each segment, Gabor 1946.
Followed by other TFRs, which differed from each other by the
selection of the windowing function
Short Time Fourier Transform
(STFT)
1.
2.
3.
4.
5.
6.
Choose a window function of finite length
Place the window on top of the signal at t=0
Truncate the signal using this window
Compute the FT of the truncated signal, save.
Incrementally slide the window to the right
Go to step 3, until window reaches the end of the signal
For each time location where the window is centered, we
obtain a different FT
Hence, each FT provides the spectral information of a
separate time-slice of the signal, providing simultaneous time
and frequency information
STFT
Time
parameter
Frequency
parameter
Signal to
be analyzed
FT Kernel
(basis function)
STFTx (t , ) x(t ) W (t t ) e jt dt
t
STFT of signal x(t):
Computed for each
window centered at t=t’
Windowing
function
Windowing function
centered at t=t’
STFT
t’=-8
t’=-2
t’=4
t’=8
STFT at Work
STFT At Work
STFT At Work
STFT
STFT provides the time information by computing a different FTs for
consecutive time intervals, and then putting them together
Time-Frequency Representation (TFR)
Maps 1-D time domain signals to 2-D time-frequency signals
Consecutive time intervals of the signal are obtained by truncating the
signal using a sliding windowing function
How to choose the windowing function?
What shape? Rectangular, Gaussian, Elliptic…?
How wide?
Wider window require less time steps low time resolution
Also, window should be narrow enough to make sure that the portion of
the signal falling within the window is stationary
Can we choose an arbitrarily narrow window…?
Selection of STFT Window
STFTx (t , ) x(t ) W (t t ) e jt dt
t
Two extreme cases:
W(t) infinitely long: W (t ) 1
STFT turns into FT, providing
excellent frequency information (good frequency resolution), but no time
information
W(t) infinitely short:
W (t ) (t )
jt
jt
STFTx (t , ) x(t ) (t t ) e
t
dt x(t ) e
STFT then gives the time signal back, with a phase factor. Excellent
time information (good time resolution), but no frequency information
Wide analysis window poor time resolution, good frequency resolution
Narrow analysis windowgood time resolution, poor frequency resolution
Once the window is chosen, the resolution is set for both time and frequency.
Heisenberg Principle
t f
Time resolution: How well
two spikes in time can be
separated from each other in
the transform domain
1
4
Frequency resolution: How
well two spectral components
can be separated from each
other in the transform domain
Both time and frequency resolutions cannot be arbitrarily high!!!
We cannot precisely know at what time instance a frequency component is
located. We can only know what interval of frequencies are present in which time
intervals
The Wavelet Transform
Overcomes the preset resolution problem of the STFT by using a
variable length window
Analysis windows of different lengths are used for different
frequencies:
Analysis of high frequencies Use narrower windows for
better time resolution
Analysis of low frequencies Use wider windows for better
frequency resolution
This works well, if the signal to be analyzed mainly consists of slowly
varying characteristics with occasional short high frequency bursts.
Heisenberg principle still holds!!!
The function used to window the signal is called the wavelet
The Wavelet Transform
A normalization
Translation parameter, Scale parameter,
Signal to be
measure of time
measure of frequency constant
analyzed
CWTx ( , s ) x ( , s )
Continuous wavelet transform
of the signal x(t) using the
analysis wavelet (.)
t
x
t
dt
s
s t
1
The mother wavelet. All kernels are
obtained by translating (shifting) and/or
scaling the mother wavelet
Scale = 1/frequency
CWTx ( , s )
x ( , s )
t
x
t
dt
s
s t
1
WT at Work
Low frequency
(large scale)
High frequency
(small scale)
WT at Work
WT at Work
WT at Work
Matlab Demos on CWT
Discrete Wavelet Transform
CWT computed by computers is really not CWT, it is a discretized
version of the CWT.
The resolution of the time-frequency grid can be controlled (within
Heisenberg’s inequality), can be controlled by time and scale step sizes.
Often this results in a very redundant representation
How to discretize the continuous time-frequency plane, so that the
representation is non-redundant?
Sample the time-frequency plane on a dyadic (octave) grid
CWTx ( , s )
x ( , s )
t
x
t
dt
s
s t
1
kn (t ) 2 k 2 k n
k, n Z
Discrete Wavelet Transform
Dyadic sampling of the time –frequency plane results in a
very efficient algorithm for computing DWT:
Subband coding using multiresolution analysis
Dyadic sampling and multiresolution is achieved through a
series of filtering and up/down sampling operations
x[n]
H
y[n]
y[n ] x[n ] * h[n ] h[n ] * x[n ]
N
x[k ] h[n k ]
k 1
N
h[k ] x[n k ]
k 1
Discrete Wavelet Transform
Implementation
~
yhigh[k ] g[n 2k ]
yhigh[k ] x[n] g[n 2k ]
x[n]
k
n
~
G
2
~
H
2
~
ylow[k ] x[n] h[ n 2k ]
~
G
2
2
G
~
H
2
2
H
+
2
G
2
H
+
yhigh[k ] g[n 2k ]
k
n
Decomposition
x[n]
Reconstruction
G
Half band high pass filter
2 Down-sampling
H
Half band low pass filter
2 Up-sampling
2-level DWT decomposition. The decomposition can be continues as long
as there are enough samples for down-sampling.
DWT - Demystified
|H(jw)|
Length: 512
B: 0 ~
g[n]
Length: 256
B: /2 ~ Hz
h[n]
2
d1: Level 1
DWT
Coeff.
Length: 128
B: /4 ~ /2 Hz
a1
g[n]
w
|G(jw)|
h[n]
2
2
Length: 128
B: 0 ~ /4 Hz
a2
d2: Level 2
DWT
Coeff.
Length: 64
B: /8 ~ /4 Hz
Length: 256
B: 0 ~ /2 Hz
2
/2
-/2
g[n]
2
d3: Level 3
DWT
Coeff.
-
-/2
h[n]
2
Length: 64
B: 0 ~ /8 Hz
…a3….Level 3 approximation
Coefficients
/2
w
Implementation of DWT on MATLAB
Choose wavelet
and number
of levels
Load signal
Hit Analyze
button
s=a5+d5+…+d1
Approx. coef.
at level 5
Level 1 coeff.
Highest freq.
(Wavedemo_signal1)
Applications of Wavelets
Compression
De-noising
Feature Extraction
Discontinuity Detection
Distribution Estimation
Data analysis
Biological data
NDE data
Financial data
Compression
DWT is commonly used for compression, since most DWT
are very small, can be zeroed-out!
Compression
Compression
ECG- Compression
Denoising Implementation
in Matlab
First, analyze
the signal with
appropriate
wavelets
Hit
Denoise
(Noisy Doppler)
Denoising Using Matlab
Choose thresholding
method
Choose noise type
Choose thrsholds
Hit
Denoise
Denosing Using Matlab
Discontinuity Detection
(microdisc.mat)
Discontinuity Detection
with CWT
(microdisc.mat)
Application Overview
Data Compression
Wavelet Shrinkage Denoising
Source and Channel Coding
Biomedical Engineering
EEG, ECG, EMG, etc analysis
MRI
Nondestructive Evaluation
Ultrasonic data analysis for nuclear power plant pipe inspections
Eddy current analysis for gas pipeline inspections
Numerical Solution of PDEs
Study of Distant Universes
Galaxies form hierarchical structures at different scales
Application Overview
Wavelet Networks
Real time learning of unknown functions
Learning from sparse data
Turbulence Analysis
Analysis of turbulent flow of low viscosity fluids flowing at high speeds
Topographic Data Analysis
Analysis of geo-topographic data for reconnaissance / object identification
Fractals
Daubechies wavelets: Perfect fit for analyzing fractals
Financial Analysis
Time series analysis for stock market predictions
History Repeats Itself…
1807, J.B. Fourier:
All periodic functions can be expressed as a weighted sum of
trigonometric function
Denied publication by Lagrange, Legendre and Laplace
1822: Fourier’s work is finally published
…
…
143 years
…
…
1965, Cooley & Tukey: Fast Fourier Transform
History Repeats Itself:
Morlet’s Story
1946, Gabor: STFT analysis:
high frequency components using a narrow window, or
low frequency components using a wide window, but not both
Late 1970s, Morlet’s (geophysical engineer) problem:
Time - frequency analysis of signals with high frequency components for short
time spans and low frequency components with long time spans
STFT can do one or the other, but not both Solution: Use different
windowing functions for sections of the signal with different frequency
content
Windows to be generated from dilation / compression of prototype small,
oscillatory signals wavelets
Criticism for lack of mathematical rigor !!!
Early 1980s, Grossman (theoretical physicist): Formalize the transform and
devise the inverse transformation First wavelet transform !
Rediscovery of Alberto Calderon’s 1964 work on harmonic analysis
1980s
1984, Yeves Meyer :
Similarity between Morlet’s and Colderon’s work, 1984
Redundancy in Morlet’s choice of basis functions
1985, Orthogonal wavelet basis functions with better time and
frequency localization
Rediscovery of J.O. Stromberg’s 1980 work the same basis
functions (also a harmonic analyst)
Yet re-rediscovery of Alfred Haar’s work on orthogonal
basis functions, 1909 (!).
Simplest known orthonormal wavelets
Transition to the
Discrete Signal Analysis
Ingrid Daubechies:
Discretization of time and scale parameters
of the wavelet transform
Wavelet frames, 1986
Orthonormal bases of compactly
supported wavelets (Daubechies wavelets),
1988
Liberty in the choice of basis functions at
the expense of redundancy
Stephane Mallat:
Multiresolution analysis w/ Meyer, 1986
Ph.D. dissertation, 1988
Discrete wavelet transform
Cascade algorithm for computing DWT
…However…
Decomposition of a discrete into dyadic frequencies (MRA) , known to
EEs under the name of “Quadrature Mirror Filters”, Croisier, Esteban and
Galand, 1976 (!)
Transition to the Discrete Signal
Analysis
Martin Vetterli & Jelena Kovacevic
Wavelets and filter banks, 1986
Perfect reconstruction of signals
using FIR filter banks, 1988
Subband coding
Multidimensional filter banks, 1992
1990s
Equivalence of QMF and MRA, Albert Cohen, 1990
Compactly supported biorthogonal wavelets, Cohen, Daubechies, J.
Feauveau, 1993
Wavelet packets, Coifman, Meyer, and Wickerhauser, 1996
Zero Tree Coding, Schapiro 1993 ~ 1999
Search for new wavelets with better time and frequency localization
properties.
Super-wavelets
Matching Pursuit, Mallat, 1993 ~ 1999
New & Noteworthy
Zero crossing representation
signal classification
computer vision
data compression
denoising
Super wavelet
Linear combination of known basic wavelets
Zero Tree Coding, Schapiro
Matching Pursuit , Mallat
Using a library of basis functions for decomposition
New MPEG standard