Wavelet analyses and applications

14 Pages·2009·0.5 MB·English

IOPPUBLISHING EUROPEANJOURNALOFPHYSICS
Eur.J.Phys.30(2009)1049–1062 doi:10.1088/0143-0807/30/5/013
Wavelet analyses and applications
Cristian C Bordeianu1, Rubin H Landau2 and
Manuel J Paez3
1FacultyofPhysics,UniversityofBucharest,Bucharest,RO077125,Romania
2DepartmentofPhysics,OregonStateUniversity,Corvallis,OR97331,USA
3DepartmentofPhysics,UniversityofAntioquia,Medellin,Colombia
E-mail:cristian.bordeianu@brahms.fizica.unibuc.ro,[email protected]
mpaez@fisica.udea.edu.co
Received27May2009,infinalform16June2009
Published20July2009
Onlineatstacks.iop.org/EJP/30/1049
Abstract
It is shown how a modern extension of Fourier analysis known as wavelet
analysis is applied to signals containing multiscale information. First, a
continuouswavelettransformisusedtoanalysethespectrumofanonstationary
signal(onewhoseformchangesintime). Thespectralanalysisofsuchasignal
givesthestrengthofthesignalineachfrequencyasafunctionoftime. Next,the
theoryisspecializedtodiscretevaluesoftimeandfrequency,andtheresulting
discrete wavelet transform is shown to be useful for data compression. This
paper is addressed to a broad community, from undergraduate to graduate
studentstogeneralphysicistsandtospecialistsinotherfieldsthanwavelets.
1. Introduction
Fourier decomposition is a standard tool in physics and is excellent for analysing periodic
signals whose functional forms do not change in time (stationary signals). However, if the
signalisnonstationary,suchasthatinfigure 1,
y(t)=sin2πt, for 0(cid:2)t (cid:2)2, (1)
=5sin2πt +10sin4πt, for 2(cid:2)t (cid:2)8, (2)
=2.5sin2πt +6sin4πt +10sin6πt, for 8(cid:2)t (cid:2)12, (3)
then the standard Fourier analysis may tell us which frequencies are present, but it will not
tellwhenthedifferentfrequencieswerepresent. Clearly,ifwewanttoobtaintimeresolution
aswellasfrequencyresolution,thenweneedtointroduceanothervariableintotheanalysis,
for otherwise the standard Fourier analysis has all of the√constituent frequencies occurring
simultaneously. In addition, because the basis set eiωt/ 2π extends over an infinite time
0143-0807/09/051049+14$30.00 (cid:2)c 2009IOPPublishingLtd PrintedintheUK 1049
1050 CCBordeianuetal
0 1 2 3 4 5 6 7 8 9
Figure1. Thetimesignal(1)containinganincreasingnumberoffrequenciesastimeincreases.
Theboxesarepossibleplacementsofwindowsforshort-timeFouriertransforms.
domainwithaconstantamplitude,thedeductionoftheFouriertransformY(ω),theamount
ofeiωt inasignal,requiresintegrationoverinfinitetime:
(cid:2)
+∞ e−iωt
Y(ω)= dt√ y(t)=(cid:4)ω|y(cid:5). (4)
−∞ 2π
We recognize the transform as the overlap of the basis function with signal y or the ω
representationofthesignal(thenegativeexponentialappearsbecause(cid:4)ω|occursandnot|ω(cid:5)).
Accordingly, the value of Y(ω) at any one ω is correlated to its value at other ω’s. This
latterfactoftenmakesitdifficulttoreconstructacomplicatedsignalwithoutincludingavery
largenumberofFouriercomponents,whichinturnrequiresthelengthycomputationofmany
componentsandtheefficientstorageofverylargedatafiles.
In this paper, we outline an extension of Fourier analysis known as wavelet analysis.
The purpose of the paper is not to present new research results, but rather to introduce the
basicideasofwaveletanalysistoabroadercommunityandindicatehowitovercomessome
of the shortcomings of Fourier analysis by adding another variable to the theory. Indeed,
wavelets have found value in areas as diverse as time-dependent quantum mechanics, brain
waves and √gravitational waves. Basically, rather than expanding a signal in terms of basis
waveseiωt/ 2π withinfiniteextent,weexpandintermsofbasiswavepacketsoffiniteextent,
theso-calledwavelets(examplesinfigure 2).
Bycentringdifferentwaveletsatdifferenttimes(figure 3(left)),weareabletodeducethe
timeresolutionofeachfrequencyintheinputsignal. Byadjustingthewaveletsateachtimeto
containdifferentfrequencies(figure 3(right)),weareabletodeducethefrequencydependence
ofthesignalateachtime. Andbecausethededucedwaveletcomponentsaresortedaccording
totimeandfrequencywithouthighcorrelation,foragivenlevelofresolutionrequiredinthe
reconstructedsignal,itispossibletoeliminatethecomponentsthatcontributeonlyforhigher
resolution;thispermitssignificantdatacompression(itisthemathematicsbehindJPEG-2000
standard). The subject of this paper is how this is done. In section 2, we provide some
backgroundinformationonshort-timeFouriertransforms,wavepacketsandfilters. Insection
3, we formulate the wavelet transform for continuous times and frequencies (CWT), and in
section 4 we formulate the wavelet transform for discrete times and frequencies (DWT). In
section5, weexplain thepyramidschemeimplementationandinsection6wecalculate the
Daub4 filter coefficients. While the CWT displays the basic principles of wavelet analysis
in a straightforward way, it is not optimized for computation. The DWT is optimized for
computation,butatthepriceoftechnicalcomplicationsthatsomereadermayfindchallenging
onfirstreading.
Waveletanalysesandapplications 1051
0.5
ψ
ψ0.0
0.0
–0.5
–1.0
–6 –4 –2 0 2 4 6 –4 0 4
t
t
Daub4 e6
0.1
ψ0.0 0
–1.0 –0.1
–4 0 4 0 1000
t
Figure 2. Four sample mother wavelets. Clockwise from top: Morlet (real part), Mexican
hat,Daub4e6(explainedlater)andHaar. Waveletbasisfunctionsaregeneratedbyscalingand
translatingthesemotherwavelets.
s = 1,τ= 0 1.0
0.5
0.0 0.0
–0.5
–1.0
–1.0
–6 –4 –2 0 2 4 6 –6 –4 –2 0 2 4 6
t t
s = 1,τ= 6 0.6 s = 2,τ= 0
0.0 0.0
–0.6
–1.0
–4 –2 0 2 4 6 8 10 –6 –4 –2 0 2 4 6
t t
Figure3. Fourwaveletbasisfunctionsgeneratedbyscaling(s)andtranslating(τ)theoscillating
Gaussianmotherwavelet. Notehows<1isawaveletwithhigherfrequency,whiles>1hasa
lowerfrequencythanthes=1mother.Likewise,theτ =6waveletisjustatranslatedversionof
theτ =0onedirectlyaboveit.
2. Short-timetransformsandwavelets
Aswehavesaid,becausetheFourier√transform(4)extendsoverinfinitetimewithaconstant
amplitudeforthebasisfunctioneiωt/ 2π,ifweapplyittoasignalsuchasthatinfigure 1,
inwhichtherearedifferentfrequenciespresentatdifferenttimes,weloseallinformationas
1052 CCBordeianuetal
towhenintimeaparticularfrequencyoccurred. Anobvioussolutiontothisproblem,known
as the short-time Fourier transform, was a precursor to wavelet analysis. An example of it
wouldbetoperformaseparatetransformforeachoftheboxesorwindowsshowninfigure 1,
andthusobtainadifferentspectrumforeachtimeinterval. Ratherthan‘chopup’asignalby
handforeachcase,weintroduceawindowfunctionw(t)thatdiffersfrom1foronlyashort
periodoftimeorlifetime(cid:5)t aroundt =0(imagineaboxoraGaussian). Thenthefunction
w(t −τ)wouldbeawindow centredattimeτ. Usingthiswindow function, wedefine the
short-timeFouriertransforma(cid:2)s
+∞ e−iωt
Y(ST)(ω,τ)= dt√ w(t −τ)y(t), (5)
−∞ 2π
wheretheexactformofthewindowfunctionisnotessentialaslongasithasashortlifetime
andsocanbethoughtofasatransparentboxonanopaquebackground. In(5),thevalueof
the translation time τ corresponds to the location of the window w over the signal, so that
anysignalwithinthewindowwidthistransformed,whilethesignallyingoutsidethewindow
does not enter. Evidently, the short-time transform is a function of two variables, and so a
surfaceor3Dplotisneededtoviewitasafunctionofbothωandτ.
Intheshort-termFouriertransform(5),wemultiplyawindowfunctionw(t−τ),whichis
finiteintimeabouttimeτ,bytheexponentiale−iωt,whichisoscillatoryintimewithfrequency
ω. Inwaveletanalysis,wecombinebothtypesofbehaviourintoawaveletthathasfinitewidth
intime(cid:5)t andalsooscillates, butmayhave avarietyofforms[1]. Forexample, awavelet
maybeanoscillatingGaussian(Morlet: topleftinfigure 3),
(cid:6)(t)=e2πite−t2 =(cos2πt +isin2πt)e−t2, (6)
thesecondderivativeofaGaussian(Mexicanhat,toprightinfigure 3),
d2
(cid:6)(t)=− e−t2 =(1−t2)e−t2, (7)
dt2
anup-and-downstepfunction(lowerleftinfigure 3)orafractal-likeshape(bottomright).
Ofcourse,oscillatoryfunctionsthatarenonzeroforfinitelengthsoftime(cid:5)tarejustwhat
arecommonlycalledwavepackets,andarefamiliarfromtime-dependentquantummechanics
and signal processing. It is well known [2] that the Fourier transform of a wave packet is a
pulse in the frequency domain of width (cid:5)ω, with the time and frequency widths related by
theso-calledHeisenberguncertaintyprinciple,
(cid:5)t(cid:5)ω(cid:3)2π. (8)
Accordingly, as a wavelet is made more localized in time (smaller (cid:5)t), it becomes less
localizedinfrequency(larger(cid:5)ω),andviceversa. Forexample,thefunctiony(t)=sinω t
0
iscompletelylocalizedinfrequency,(cid:5)ω=0,andsohasaninfiniteextentintime,(cid:5)t (cid:6)∞,
whilethewaveletsinfigure 2havefinite(cid:5)ωand(cid:5)t.
Let us now take this one step further. We want to use wavelet basis functions that
oscillate in time and also contain a range of frequencies. While the example functions in
(6) and (7) are evidently such functions, it is not clear how to use them as a set of basis
functions. More con√cretely, the basis functions for the short-term Fourier transform in (5)
exp(−iωt)w(t −τ)/ 2π containbothafrequencyωandtimevariableτ thatalsoappearon
the lhs, while the wavelets (cid:6)(t) we have given so far are just functions of time. Actually,
the wavelets (cid:6)(t) we have given as examples so far are what are called mother wavelets or
analysingfunctions,andeachisusedtogenerateanentiresetofbasisfunctionsordaughter
wavelets ψ (t) via a simple, yet clever, change of variable that scales and translates the
s,τ
motherwavelet: (cid:3) (cid:4)
1 t −τ
ψ (t)= √ (cid:6) . (9)
s,τ
s s
Waveletanalysesandapplications 1053
Forexample,
(cid:5) (cid:6)
1 8(t −τ)
(cid:6)(t)=sin(8t)e−t2/2 ⇒ ψ (t)= √ sin e−(t−τ)2/2s2. (10)
s,τ
s s
Translating and scalingone ofthese mother wavelets generate anentiresetofchild wavelet
basisfunctions,witheachindividualfunctioncoveringadifferentfrequencyrangeatadifferent
time. Weseethatlargerorsmallervaluesofsexpandorcontractthemotherwavelet,while
different values of τ shift the centre of the wavelet. Because the wavelets are inherently
oscillatory, the scaling leads to the same number of oscillations occurring in different time
spa√ns,whichisequivalenttohavingbasisstateswithdifferingfrequencies. Herethedivision
by s ismadetoensurethatthereisequal‘power’(orenergyorintensity)ineachregionof
s,althoughothernormalizationscanalsobefoundintheliterature.
Notethatwehavemadeachangeinnotationthatisstandardforwavelets,butsomewhat
confusingatfirst. Ratherthanusingthefrequencyvariableω wehaveswitchedtothescale
variable s, which has the dimension of time. You may think of the two as simply inversely
related:
2π 2π
ω= , s = (scale–frequencyrelation). (11)
s ω
If we are interested in the time details of a signal, we would say that we are interested in
whatishappeningatsmallvaluesofthescales;equation(11)indicatesthatsmallvaluesofs
correspondtohigh-frequencycomponentsofthesignal. Thatbeingthecase,thetimedetails
ofthesignalareinthehigh-frequency,orlow-scale,components.
3. Thecontinuouswavelettransform
ThewavelettransformY(s,t)ofatimesignaly(t)isdefinedmuchliketheFouriertransform
(4),butwithawaveletbasisratherthananexponential:
(cid:2) (cid:2) (cid:3) (cid:4)
+∞ +∞ 1 t −τ
Y(s,τ)= dtψ∗ (t)y(t)= dt√ (cid:6) y(t)=(cid:4)s,t|y(cid:5). (12)
−∞ s,τ −∞ s s
Because each wavelet is localized in time, each acts as its own window function. Because
eachwaveletisoscillatory,eachcontainsitsownsmallrangeoffrequenciesorscalesvalues.
Equation (12) says that the wavelet transform Y(s,τ) is a measure of the amount of basis
functionψ (t)presentinsignaly(t). Theτ variableindicatesthetimeportionofthesignal
s,τ
beingdecomposed,whilethescalevariablesisequivalenttothefrequencypresentduringthat
time. Althoughwedonotderiveithere,theinversetransformis
(cid:2) (cid:2)
1 +∞ +∞ ψ∗ (t)
y(t)= dτ ds s,τ Y(s,τ), (13)
C −∞ 0 s3/2
wherethenormalizationconstantCdepends onthespecificwaveletclassused. Infigure 4,
we show a 2D and a 3D representation of the results of evaluating numerically the wavelet
transform (12) for the sample signal (1). The Java code used for the transform is from our
textbook[4]. Recallthatthissignalcontainsonefrequencyduringthefirstwindowoftime,
two frequencies during the second time window and three frequencies for the last window.
Indeed,atshorttimesτ weseeasinglehumpathighscales (cid:6)0.8(lowfrequency);atlatter
timesthisisjoinedbyanotherhumpatsmallers(higherfrequency),withathirdpeakentering
atthelatesttimeswithayetsmallersvalue(higherfrequency). Althoughitishardtoevaluate
relativeheightsfromthisplot,therelativestrengthsofthethreefrequencycomponentsappear
tobeinagreementwith(1).
1054 CCBordeianuetal
1 Y
0
–1
0
4
8
0
12 1 s
2
Figure4. Left: thecontinuouswaveletspectrum3Dobtainedbyanalysingtheinputsignalwith
Morletwavelets. Observehowatsmallvaluesoftimeτ thereispredominantlyonefrequency
present;howasecond,higherfrequency(smaller-scale)componententersatintermediatetimes
andhowatlargertimes,astillhigherfrequencycomponententers. Right: thespectrogram2Dof
theinputsignalmakesitclearhowthereisone,twoandthenthreefrequenciespresent.
4. Thediscretewavelettransform
As is true for discrete Fourier transforms, if a time signal y is measured at only N discrete
times,
y(t )≡y , m=1,...,N, (14)
m m
then we can determine only N independent components of the transform Y. The discrete
wavelet transform (DWT) is a clever way to compute only the N independent components
requiredtoreproducetheinputsignal,consistentwiththeuncertaintyprinciple,byevaluating
the transforms at discrete values for the scaling parameter s and the time translation
parameterτ:
(cid:6)[(t −k2j)/2j] (cid:6)(t/2j −k)
ψ (t)= √ ≡ √ , (15)
j,k
2j 2j
k
s =2j, τ = , k,j =0,1,.... (16)
2j
Notethatincontrasttothecontinuouswavelettransform(12),theDWT(15)employsadiscrete
set of basis wavelets denoted by the integer subscripts j and k, whose maximum values are
yettobedetermined. Also,wehaveassumedthatthetotaltimeintervalT = 1,sothattime
isalwaysmeasuredinintegervalues. Thischoiceofsandτ basedonpowersof2iscalleda
dyadicgridarrangementandisseentoautomaticallyperformthescalingsandtranslationsat
thedifferenttimescalesthatareattheheartofwaveletanalysis(somereferencesscaledown
withincreasingj,incontrasttoourscalingup). Wenowusethetrapezoidintegrationruleto
expresstheDWTasthesimplesum
(cid:2)
+∞ (cid:7)
Y = dtψ (t)y(t)(cid:6) ψ (t )y(t )h (DWT). (17)
j,k j,k j,k m m
−∞
m
Foranorthonormalwaveletbasis,thentheinversediscretetransformis[1]
(cid:7)+∞
y(t)= Y ψ (t) (inverseDWT). (18)
j,k j,k
j,k=−∞
Waveletanalysesandapplications 1055
y
c
n
e
u
q
e
Fr
Time
Figure5. Timeandfrequencyresolutions. Eachboxrepresentsanequalportionofthetime–
frequencyplanebutwithdifferentproportionsoftimeandfrequency.
It is interesting to note that this inversion will exactly reproduce the input signal at the N
inputpointsifwesumoveraninfinitenumberofwavelets,butwillbelessexactinpractical
calculations.
Notein(15)and(17)thatupuntilthispointwehavekeptthetimevariabletinthewavelet
basisfunctionscontinuous,eventhoughsandτ arediscrete. Thisisusefulinestablishingthe
orthonormalityofthebasisfunctions,
(cid:2)
+∞
dtψj∗,k(t)ψj(cid:10),k(cid:10)(t)=δjj(cid:10)δkk(cid:10), (19)
−∞
whereδ istheKroneckerdeltafunction. Havingψ (t)’snormalizedto1meansthateach
m,n j,k
waveletbasishas‘unitenergy’;beingorthogonalmeansthateachbasisfunctionisindependent
oftheothers. Becausewaveletsarelocalizedintime,thedifferenttransformcomponentshave
lowlevelsofcorrelationwitheachother. Altogether, thisleadstoefficientandflexibledata
storage.
The use of a discrete wavelet basis makes it clear that we need only sample the input
signalatthediscretevaluesoftimedeterminedbytheintegersj andk. Ingeneral,youwant
time steps that sample the signal at enough times in each interval to obtain a reconstructed
signalwithapredeterminedlevelofprecision. Aruleofthumbistostartwith100stepsto
covereachmajorfeature. Ideally,theneededtimescorrespondtothetimesatwhichthesignal
wassampled,althoughthismayrequiresomeforethought.
Considerfigure 5. Wemeasureasignalatanumberofdiscretetimeswithintheintervals
(korτ values)correspondingtotheverticalcolumnsoffixedwidthalongthetimeaxis. For
eachtimeinterval,wewanttosamplethesignalatanumberofscales(frequenciesorjvalues).
However,thebasicmathematicsofFouriertransformsindicatesthatthewidth(cid:5)t ofawave
packet ψ(t) and the width (cid:5)ω of its Fourier transform Y(ω) are related by an uncertainty
principle
(cid:5)ω(cid:5)t (cid:3)2π. (20)
This relation constrains the number of times we can meaningfully sample a signal in order
to determine a number of Fourier components. So while we may want a high-resolution
1056 CCBordeianuetal
NSamples
Input
N/2 N/2
c(1) L H d(1)
Coefficients Coefficients
N/4 N/4
c(2) L H d(2)
Coefficients Coefficients
N/8 N/8
c(3) L H d(3)
Coefficients Coefficients
2 2
c(n) L H d(n)
Coefficients Coefficients
Figure 6. A signal is processed by high (H) and low (L) band filters, and the outputs are
downsampled(reducedbyafactorof2)bykeepingeveryotherpoint. Thefilteringcontinues
untilthereareonlytwoHpointsandtwoLpoints.Thetotalnumberofoutputdataequalsthetotal
numberofsignalpoints.
reproductionofoursignal,wedonotwanttostoremoredatathanareneededtoobtainthat
reproduction. Ifwesamplethesignalfortimescentredaboutsomeτ inanintervalofwidth
(cid:5)τ (figure 5) and compute the transform at a number of scales s or frequencies ω = 2π/s
covering a range of height (cid:5)ω, then the relation between the height and width is restricted
bytheuncertaintyrelation,whichmeansthateachoftherectanglesinfigure 5hasthesame
area(cid:5)ω(cid:5)t = 2π. Theincreasingheightsoftherectanglesathigherfrequenciesmeanthat
alargerrangeoffrequenciesshouldbesampledateachtimeasthefrequencyincreases. The
inherent assumption in wavelet analysis is that the low-frequency components provide the
grossorsmoothoutlineofthesignalwhich,beingsmooth,doesnotrequiremuchdetail,while
thehigh-frequencycomponentsgivethedetailsofthesignaloverashorttimeintervalandso
requiremanycomponentsinordertorecordthesedetailswithhighresolution.
Industrial-strength wavelet analyses do not compute explicit integrals but instead apply
a filtering technique known as multiresolution analysis (figure 6). It is a type of pyramid
algorithm that takes the N sample values of a signal and passes it successively through a
numberoffilters. SinceeachfilterisequivalenttoaDWT,thefinaloutputisthetransformat
allscalesandtimes. Herewedefineafilterasadevicethatconvertsaninputsignalf(t)toan
outputsignalg(t)viaanintegrationoverthesignal:
(cid:2)
+∞
g(t)= dτf(τ)h(t −τ)=f(t)∗h(t). (21)
−∞
Here,thefunctionh(t)iscalledtheresponseortransferfunctionofthefilterandtheoperation
in (21) iscalled aconvolution and isdenoted by an asterisk*. Equation (21)states that the
outputofafilterequalstheinputconvolutedwiththetransferfunction,whiletheconvolution
Waveletanalysesandapplications 1057
theoremstatesthattheFouriertransformoftheconvolutiong(t)isproportionaltotheproduct
ofthetransformsoff(t)andh(t):
√
G(ω)= 2πF(ω)H(ω). (22)
Acomparisonofthedefinitionofafiltertothedefinitionofawavelettransform(12)shows
that transforming a signal is equivalent to filtering it, with each filter outputting a different
scale. Theexplicitoutputisthesumoftheproductsofintegrationweightswiththevalueofthe
waveletsattheintegrationpoints,andeachproductiscalledafiltercoefficientc . Therefore,
i
ratherthantabulatingexplicitwaveletfunctions,asetoffiltercoefficientsisallthatisneeded
toperformdiscretewavelettransforms.
5. Pyramidschemeimplementation
Thepyramidalgorithmdecomposesasignalintolow-frequencyorsmoothcomponentss and
i
high-frequency or detailed components d . Because high-resolution reproduction requires
i
moreinformationaboutthedetailsinasignalthanaboutitsgrossshape,mostofthesmooth
components can be discarded and thus the data compressed. Furthermore, because the
components of different resolutions are independent of each other, it is possible to adjust
the resolution by changing the number of high resolution components stored. So when we
lookatfigure 6wecanthinkofeachLorHfilterasseparatingoutthesmoothanddetailed
componentsofthesignalrespectively, withrepeatedfilteringseffectivelyloweringthescale
togetmoredetailedinformation(thereasonthistechniqueiscalledmultiresolutionanalysis).
Each ↓ filters out half of the signal, a process known as factor-of-2 decimation, and this
producesthedatacompression. Theendresultisavectorwiththecomponentsofthewavelet
transform.
Inpractice,thefiltersLandHarefiltermatricesthatmultiplytheinputvector(wediscuss
someinthefollowingsection). Thisisfollowedbyareductionoftheoutputbyone-half,in
whichthelow-frequencypartsarediscarded,andthenareorderingoftheoutputsothatitmay
befilteredagain. Forexample,letusconsiderasignalvectoroflengthN = 8,whichwould
alsocorrespondtothelasttwostepsinthepyramidalgorithm:
⎛ ⎞ ⎛ ⎞ ⎛ ⎞ ⎛ ⎞
⎛y1⎞ ⎜s1(1)⎟ ⎜s1(1)⎟ ⎜s1(2)⎟ ⎜s1(2)⎟
⎜⎜⎜⎜⎜⎜yy23⎟⎟⎟⎟⎟⎟ ⎜⎜⎜⎜⎜⎜ds21((11))⎟⎟⎟⎟⎟⎟ ⎜⎜⎜⎜⎝ss23((11))⎟⎟⎟⎟⎠ fi−l→ter ⎜⎜⎜⎜⎝ds21((22))⎟⎟⎟⎟⎠ o−r→der ⎜⎜⎜⎜⎝ds21((22))⎟⎟⎟⎟⎠
⎜⎜⎜y4⎟⎟⎟fi−l→ter⎜⎜⎜d2(1)⎟⎟⎟o−r→der ⎛s4(1)⎞ ⎛d2(2)⎞ ⎛d2(2)⎞ . (23)
⎜⎜y5⎟⎟ ⎜⎜s3(1)⎟⎟ ⎜d1(1)⎟ ⎜d1(1)⎟ ⎜d1(1)⎟
⎜⎜⎜y6⎟⎟⎟ ⎜⎜⎜d3(1)⎟⎟⎟ ⎜⎜⎜d2(1)⎟⎟⎟ ⎜⎜⎜d2(1)⎟⎟⎟ ⎜⎜⎜d2(1)⎟⎟⎟
⎝y7⎠ ⎜⎝s4(1)⎟⎠ ⎜⎝d3(1)⎟⎠ ⎜⎝d3(1)⎟⎠ ⎜⎝d3(1)⎟⎠
y8 d(1) d(1) d(1) d(1)
4 4 4 4
Thefirstfilteringleadstoanoutputvectorwiththesmoothanddetailedcomponentsintermixed.
It gets reordered into the lower vector of just detailed components, which is saved, and a
smallervectoroffirst-ordersmoothedcomponents,whichthengetfilteredagain. Thissecond
filteringleadstointermixedsecond-orderfilteredcomponents,whichthengetreorderedinto
just smooth and detailed components. The process ends when there are just two smooth
componentsleft,aswehavehere.
1058 CCBordeianuetal
Asarealisticillustration,imaginethatwehavesampledthechirpsignaly(t)=sin(60t2)
at1024times. Thesuccessivefilteringofthissignalisillustratedschematicallyasapassage
fromthetoptothebottomoffigure 6,orwithactualdatainfigure 7. First,theoriginal1024
signal samples are fi(cid:14)ltered(cid:15)and ordered into 512 detailed and 512 smooth components. All
detailedcomponents d(1) aresaved,andtheremaining512once-filteredsmoothcomponents
(cid:14) (cid:15) i
s(1) areprocessedfurther. Thisstepintheprocessisknownasdownsampling. Weseein
i
figure 7 that at this stage the smooth components still contain many high-frequency parts,
whilethedetailedcomponentsaresmallerinmagnitude. Thesetwoset(cid:14)sare(cid:15)shownbelowthe
original signal. In the next level down, the 512 smooth components s(1) are filtered and
(cid:14) (cid:15) i
ordered. The resulting 256, twice-filtered detailed components d(2) are saved, while the
(cid:14) (cid:15) i
256smoothcomponents s(2) arefurtherprocessed,withthesecomponentsdisplayedlower
i
downinfigure 6. Theresultingoutputfromeachsuccessivestepissimilar,butwithcoarser
features for the smooth coefficients and larger values for the details. Note that in the upper
graphswehaveconnectedthepointstomaketheoutputlookcontinuous,whileinthelower
graphs,withfewerpoints,wehaveplottedtheoutputashistogramstomakethepointsmore
evident. Theprocesscontinuesuntilthereareonlytwosmoothandtwodetailedcoefficients
left. Since this last filtering is done with the broadest wavelet, it is of the lowest resolution
andthereforerequirestheleastinformation. Whenwearedone,therewillbeatotalof512+
256+128+64+32+16+4=1024transformvaluessaved,thesamenumberasthenumber
oftimesignalsmeasured.
Toreconstructtheoriginalsignalfromthe1024transformvalues(calledsignalsynthesis),
areversedprocessisfollowed. Weessentiallyfollowfigure 6;onlynowthedirectionsofall
the arrows are reversed and the inverses to the filter matrices are used. We begin with the
lastsequenceoffourcoefficients,upsamplethem,passthemthroughinversefilterstoobtain
new levels of coefficients and repeat until all of the N = 1024 values of the original signal
arerecovered. Ifthedatawerecompressedbyeliminatingsomeofthehigherorderdetailed
components,thenwewouldstarttheprocesspartiallyupthepyramid.
6. Daubechieswaveletsfiltering
We have spoken about how the discrete wavelet transformation can be implemented as
digital filtering, with the filtering being simply a matrix multiplication of the signal vector
by an appropriate filter matrix. We now indicate how one of the standard filter matrices,
correspondingtotheDaubechieswaveletbases,wasdeterminedbytheBelgianmathematician
Ingrid Daubechies in 1988 [3]. To keep things simple, we discuss just the Daub4 class
containing the four coefficients c ,c ,c and c . We deduce how these coefficients can be
0 1 2 3
used to filter just four input signal values {y ,y ,y ,y }. We have already seen that the
1 2 3 4
pyramidalgorithmfordeterminingaDWTisbasedonaseriesofhigh-andlow-passfilters,
so we start by trying to express these two operations in terms of the filter coefficients. The
basic idea is that replacing each signal value by the average of the values surrounding it
tends to smooth the signal, and thus acts as the low-pass filter L. Likewise, replacing each
signal value by the differences from surrounding values tends to emphasize the variation or
detailsinthesignal,andthusactsasthehigh-passfilterH.Forexample,ifthesignalwerea
sawtooth function with sharp points, the average would be rounded. Likewise, if the signal
wereconstant,thedifferenceswouldallvanish,astherearenodetailsinthesignal. Achoice
offiltercoefficientsthatdothisis
L=(+c +c +c +c ) (24)
0 1 2 3