Fourier analysis and wavelet

13 Pages·1997·0.289 MB·English

by  Walker J.C.

walker.qxp 3/24/98 2:16 PM Page 658
Fourier Analysis
and Wavelet
Analysis
James S. Walker
In this article we will compare the classical Frequency Information, Denoising
methods of Fourier analysis with the newer
As an example of the importance of frequency in-
methods of wavelet analysis. Given a sig-
formation, we will examine how Fourier analysis can
nal, say a sound or an image, Fourier analy-
be used for removing noise from signals. Consider
sis easily calculates the frequencies and
a signal f(x)defined over the unit interval (where here
the amplitudes of those frequencies which make
x stands for time). The periodP 1 Fourier series ex-
utphpoer t thcahen ats rifgaoncrt atelhr. eiTsohtriiecsts ip coraofl vctiohdnees ssi daig ebnrraaotlai,o dwn oshv.i Hcehrov wiiesew vim eorf-, pcnan=sRi0o1nf( xo)fe −fi2i(cid:25)sn xddexfi. nEeadch bFyo urienr2 cZocenfefiic2i(cid:25)ennxt,, cwni,t ihs
an amplitude associated with the frequency nof the
although Fourier inversion is possible under
exponential ei2(cid:25)nx. Although each of these expo-
certain circumstances, Fourier methods are not
nentials has a precise frequency, they all suffer from
always a good tool to recapture the signal, par-
a complete absence of time localization in that their
ticularly if it is highly nonsmooth: too much
magnitudes, jei2(cid:25)nxj, equal 1 for all time x.
Fourier information is needed to reconstruct
To see the importance of frequency information,
the signal locally. In these cases, wavelet analy-
let us examine a problem in noise removal. In
sis is often very effective because it provides a
Figure 1(a)[top] we show the graph of the signal
simple approach for dealing with local aspects
of a signal. Wavelet analysis also provides us with
(1) f(x)=
new methods for removing noise from signals
that complement the classical methods of (5cos2(cid:25)(cid:23)x)[e−640(cid:25)(x−1=8)2
Fourier analysis. These two methodologies are +e−640(cid:25)(x−3=8)2+e−640(cid:25)(x−4=8)2
major elements in a powerful set of tools for the-
oretical and applied analysis. +e−640(cid:25)(x−6=8)2+e−640(cid:25)(x−7=8)2 ]
This article contains many graphs of discrete
where the frequency, (cid:23), of the cosine factor is 280.
signals. These graphs were created by the com-
puter program FAWAV, A Fourier–Wavelet An- Such a signal might be used by a modem for trans-
alyzer, being developed by the author. mitting the bit sequence 1011011. The Fourier co-
efficients for this signal are shown in Figure 1(b)[top].
The highest magnitude coefficients are concentrated
around the frequencies (cid:6)280. Suppose that when this
James S. Walker is professor of mathematics at the
signal is received, it is severely distorted by added
University of Wisconsin-Eau Claire. His e-mail address
noise; see Figure 1(a)[middle]. Using Fourier analy-
is [email protected].
sis, we can remove most of this noise. Computing the
The author would like to thank Hugo Rossi, Steven
noisy signal’s Fourier coefficients, we obtain the
Krantz, his colleague Marc Goulet, and two anony-
mous reviewers for their helpful comments during the graph shown in Figure 1(b)[middle]. The original sig-
writing of this article. nal’s largest magnitude Fourier coefficients are clus-
658 NOTICESOFTHEAMS VOLUME44, NUMBER6
walker.qxp 3/24/98 2:16 PM Page 659
Figure 1. (a)[top] Signal. (b)[top] Fourier coefficients of signal. (a)[middle] Signal after adding noise.
(b)[middle] Fourier coefficients of noisy signal and filter function. (b)[bottom] Fourier coefficients
after multiplication by filter function. (a)[bottom] Recovered signal.
tered around the frequency positions (cid:6)280. The for n=−1M+1;:::;0;:::;1M. The discrete
2 2
Fourier coefficients of the added noise are lo- Fourier coefficients ˆfncan be calculated by a fast
calized around the origin, and they decrease in Fourier transform (FFT) algorithm and are the
magnitude until they are essentially zero near discrete analog of the Fourier coefficients cn
the frequencies (cid:6)280. Thus, the original sig- for F, when fj =F(j=M). Moreover, ˆfn is just a
nal’s coefficients and the noise’s coefficients Riemann sum approximation of the integral that
are well separated. To remove the noise from the defines cn.1 The magnitudes of the discrete
signal, we multiply the noisy signal’s coefficients Fourier coefficients for this transient damp down
by a filter function, which is 1where the signal’s to zero very slowly (their graph is a very wide
coefficients are concentrated and 0 where the bell-shaped curve with maximum at the origin).
noise’s coefficients are concentrated. We then re- Consequently, to represent the transient well, one
cover essentially all of the signal’s coefficients; must retain most if not all of these Fourier co-
see Figure 1(b)[bottom]. Performing a Fourier efficients. In Figure 2(a)[bottom] we show the re-
series partial sum with these recovered coeffi- sults obtained from trying to compress the tran-
cients, we obtain the denoised signal, which is sPient by computing a discrete partial sum
shown in Figure 1(a)[bottom]. Clearly, the bit se- 1n0=4−104ˆfnei2(cid:25)nj using only one-fifth of the
quence 1011011can now be determined from Fourier coefficients.2 Clearly, even a moderate
the denoised signal, and the denoised signal is compression ratio of 5:1 is not effective.
a close match of the original signal. In the sec- Wavelets, however, are often very effective at
tion “Signal Denoising” we shall look at another representing transients. This is because they are
example of this method and also discuss how designed to capture information over a large
wavelets can be used for noise removal. range of scales. A wavelet series expansion of a
function f is defined by
Signal Compression X
As the example above shows, Fourier analysis is (cid:12)nk2n=2ˆ(2nx−k)
very effective in problems dealing with frequency n;k2Z
location. However, it is often very ineffective at Z
with 1
representing functions. In particular, there are (cid:12)n = f(x)2n=2ˆ(2nx−k) dx:
severe problems with trying to analyze tran- k −1
sient signals using classical Fourier methods.
The function ˆ(x)is called the wavelet, and the
For example, in Figure 2(a)[top] we show a dis-
coefficients (cid:12)n are called the wavelet coeffi-
crete signal obtained from M =1024 values k
ffj =F(j=M)gMj=−01 of the function F(x)= cients. The function 2n=2ˆ(2nx−k) is the
e−105(cid:25)(x−:6)2. For this example, we compute the 1For further discussion of discrete Fourier coefficients,
discrete Fourier series coefficients ˆfndefined by see [2] or [P11].
ˆfn =M−1PMj=−01fje−i2(cid:25)nj=M c2rTehtee sFuomur ier5n c1=o2−e5f1fi1ciˆfennetsi2, (cid:25)eqnuj,a wlsh ficjh. uses all of the dis-
JUNE/JULY1997 NOTICESOFTHEAMS 659
walker.qxp 3/24/98 2:17 PM Page 660
Figure 2. (a)[top] Signal. (a)[bottom] Fourier series using 205 coefficients, 5:1 compression. (b)[top]
Wavelet coefficients for signal. (b)[bottom] Wavelet series using only largest 4% in magnitude of
wavelet coefficients, 25:1 compression.
wavelet shrunk by a factor of 2n if nis positive The Haar Wavelet
(magnified by a factor of 2−nif nis negative) and
In order to understand how wavelet analysis
shifted by k2−n units. The factor 2n=2in the ex-
works, it is best to begin with the simplest
pression 2n=2ˆ(2nx−k)preserves the L2-norm.
wavelet, the Haar wavelet. Let 1A(x)denote the
Since the wavelet series depends on two pa-
indicator function of the set A, defined by
rameters of scale and translation, it can often be
very effective in analyzing signals. These para- 1A(x)=1 if x2A and 1A(x)=0 if x2= A. The
Haar wavelet ˆ is defined by
meters make it possible to analyze a signal’s
behavior at a dense set of time locations and with ˆ(x)=1[0;1)(x)−1[1;1)(x). It is 0 outside of
2 2
respect to a vast range of scales, thus providing [0;1), so it is well localized in time, and it sat-
the ability to zoom inon the transient behavior isfies
of the signal. For example, let us examine the ear- Z1 Z1
lier transient using a discretized version of a ˆ(x)dx=0; jˆ(x)j2dx=1:
wavelet series. We shall use a Daubechies order −1 −1
4 wavelet (Daub4 for short; see the section The Haar wavelet ˆ(x) is closely related to the
“Daubechies Wavelets”). In Figure 2(b)[top] we
function (cid:30)(x) defined by (cid:30)(x)=1 (x). This
show all of the 1024wavelet coefficients of this [0;1)
function (cid:30)(x)is called the Haar scaling function.
transient and observe that most of these coef-
Clearly, the Haar wavelet and scaling function
ficients are close to 0 in magnitude. Conse-
satisfy the identities
quently, by retaining only the largest magnitude
coefficients for use in a wavelet series, we ob-
ˆ(x)=(cid:30)(2x)−(cid:30)(2x−1);
tain significant compression. In Figure 2(b)[bot- (2)
tom] we show the reconstruction of the transient (cid:30)(x)=(cid:30)(2x)+(cid:30)(2x−1);
using only the top 4% in magnitude of the wavelet
and the scaling function satisfies
coefficients, a 25:1 compression ratio. Notice
how accurately the transient is represented. In Z1 Z1
fact, the maximum error at all computed points (cid:30)(x)dx=1; j(cid:30)(x)j2dx=1:
is less than 9:95(cid:2)10−14. There is an important −1 −1
application here to the field of signal transmis- The Haar wavelet ˆ(x) generates the system
sion. By transmitting only these 4% of the wavelet
of functions f2n=2ˆ(2nx−k)g. It is possible to
coefficients, the information in the signal can be
show directly that f2n=2ˆ(2nx−k)g is an or-
transmitted 25 times faster than if we transmit-
thonormal basis for L2(R), but it is more illu-
ted all of the original signal.This provides a con-
minating to put the discussion on an axiomatic
siderable boost in efficiency of transmission.
We shall look at more examples of com- level. This axiomatic approach leads to the
pression in the section “Compression of Sig- Daubechies wavelets and many other wavelets
nals”, but first we shall describe how wavelet as well. We begin by defining the subspaces
analysis works. fVngn2Z of L2(R)in the following way:
660 NOTICESOFTHEAMS VOLUME44, NUMBER6
walker.qxp 3/24/98 2:17 PM Page 661
n
Vn = stepfunctionsinL2(R); constant o tNhoeter mthoarte t,h esˆe˜ nw;ka(cid:17)ve0lets fhoarv e pner<iod0 , 1. aFunrd-
Tlohwisin sge tf iovfe soaunxbitoshmpeasci ne[6tse] :rfvVanlgsn[22kZns;akt2i+snf1ie);sk th2eZ fo:l- itˆ˜esrnfvy;ka +l2 ˆ˜[n0n=;;k1ˆ(˜)x, n)t;h=kef2 opnre= ar2ilˆol dk(2i2cn xHZa−aankr d)w anfvo(cid:21)erl e0tn.s O (cid:21)ˆ˜n0n t;hkaes naintd--
k=0;1;:::;2n−1. So we have the following the-
Axioms for a Multi-Resolution Analysis
orem as a consequence of Theorem 1.
(MRA)
Scaling: f(x)2Vnif and only if f(2x)2Vn+1. Theorem 2. The functions 1and ˆ˜n;kfor n(cid:21)0
and k=0;1;:::;2n−1are an orthonormal basis
Inclusion: Vn (cid:26)V(cid:26)n+1S, for (cid:27)each n. for L2[0;1).
Density:closure Vn =L2(R).
T n2Z
Maximality: Vn =f0g. Remark. In the section “Daubechies Wavelets”
n2Z we will make use of periodized scaling func-
Basis: 9(cid:30)(x)such that f(cid:30)(x−k)gk2Zis an tions, (cid:30)˜n;k, defined by
orthonormal basis for V0. X (cid:0) (cid:1)
To satisfy the basis axiom, we shall use the (4) (cid:30)˜ (x)= 2n=2(cid:30) 2n(x+j)−k
n;k
Haar scaling function (cid:30)defined above. Then, by
j2Z
combining the scaling axiom with the basis
axiom, we find that f2n=2(cid:30)(2nx−k)gk2Z is an forn(cid:21)0andk=0;1;:::;2n−1:
orthonormal basis for Vn. But the totality of all
these orthonormal bases, consisting of the set
Fast Haar Transform
f2n=2(cid:30)(2nx−k)gk;n2Z, is not an orthonormal
basis for L2(R) because the spaces Vn are not The relation between the Haar scaling function
mutually orthogonal. To remedy this difficulty, (cid:30)and wavelet ˆleads to a beautiful set of re-
we need what are called wavelet subspaces.De- lations between their coefficients as bases. Let
fine the wavelet subspace Wnto be the orthog- f(cid:11)ng and f(cid:12)ng be defined by
onal complement of Vnin Vn+1. That is, Wnsat- k kZ
isfies the equation Vn+1=Vn(cid:8)Wn where (cid:8) (cid:11)n = 1 f(x)2n=2(cid:30)(2nx−k)dx;
dspeancoetse.s F rtohme sthuem d eonf smityu tauxaiollmy aonrtdh roegpoenaatel ds uabp-- (5) k Z−11
(cid:12)n = f(x)2n=2ˆ(2nx−k)dx:
pLl2i(cRa)t=ioVn0 (cid:8)ofL t1nh=e0 lWans:t DeqecuoamtiLpoons, inwge Vo0btinai na k −1
similar way, we obtain L2(R)= n2Z Wn:Thus, Substituting 2nxin place of x in the identities
L2(R)is an orthogonal sum of the wavelet sub- in (2), we obtain
spaces Wn.
Using (2) and the MRA axioms, it is easy to 2n=2(cid:30)(2nx)= p1 [2n2+1(cid:30)(2n+1x)]
prove the following lemma. 2
Lemma 1. The functions fˆ(x−k)gk2Z are an + p12[2n2+1(cid:30)(2n+1x−1)]
orthonormal basis for the subspace W0:
f2nIt=2 fˆo(l2lnoxw−s kf)rgokm2Z tihs ea ns coarltihnogn aoxrmioaml btahsaist 2n=2ˆ(2nx)= p12[2n2+1(cid:30)(2n+1x)]
sfourm W onf. Talhl etrheef owrea,v seilnect es uLb2(sRp)aicse tsh eW onr, twhoeg hoanvael − p12[2n2+1(cid:30)(2n+1x−1)]:
obtained the following result.
Theorem 1. The functions It then follows that
f2n=2ˆ(2nx−k)gk;n2Z (cid:11)n = p1 (cid:11)n+1+ p1 (cid:11)n+1 ;
k 2 2k 2 2k+1
are an orthonormal basis for L2(R): (6)
1 1
This orthonormal basis is the Haar basisfor (cid:12)n = p (cid:11)n+1− p (cid:11)n+1 :
k 2 2k 2 2k+1
L2(R). There is also a Haar basis for L2[0;1):To
obtain it, we first define periodic wavelets ˆ˜
n;k This result shows that the nth level coefficients
by
(cid:11)n and (cid:12)n are obtained from the (n+1)stlevel
X (cid:0) (cid:1) k k
(3) ˆ˜n;k(x)= 2n=2ˆ 2n(x+j)−k : coefficients (cid:11)kn+1 through multiplication by the
j2Z following orthogonal matrix:
JUNE/JULY1997 NOTICESOFTHEAMS 661
walker.qxp 3/24/98 2:17 PM Page 662
0p1 p1 1 cient (cid:11)00 and a single wavelet coefficient (cid:12)00,
(7) O=B@ 2 2CA: and at this last step the permutation P2 is un-
p1 p−1 necessary. The complete transformation, de-
2 2 noted by H, satisfies
Successively applying this orthogonal matrix O, H =(H2(cid:8)IM−2)
we obtain f(cid:11)ng and f(cid:12)ng starting from some (cid:1)(cid:1)(cid:1)(PM=2(cid:8)IM=2)(HM=2(cid:8)IM=2)PMHM
k k
highest level coefficients f(cid:11)kNg for some large where IN is the N(cid:2)N identity matrix.
N. Because of the density axiom, Ncan be cho-
These matrix multiplications can be per-
sPen large enough to approximate f by formed rapidly on a computer. Multiplying by
k2Z(cid:11)kN2N=2(cid:30)(2Nx−k)in L2-norm as closely HM requires only O(M) operations, since HM
as desired.
consists mostly of zeroes. Similarly, the per-
Let us now discretize these results. Suppose
that we are working with data ffjgMj=−01 associ- mfourtea,t itohne P wMhroeqleu irtersa nOs(fMo)rmopaetriaotnio nrse.q Tuhierrees-
ated with the time values fxj =j=MgjM=−01on the O(M)+O(M=2)+(cid:1)(cid:1)(cid:1)+O(2)=O(M) operations.
unit interval. If we shrink the Haar scaling func- The transformation H is called a fast Haar
tion (cid:30)(x)=1 (x) enough, it covers only the transform. It should be noted that FFTs, which
[0;1)
first point, x0=0. Consequently, by choosing a have revolutionized scientific practice during
large enough N, we may assume that our scal- the last thirty years, are O(MlogM)algorithms.
ikn=g 0c;o1e;:f:f:i;cMien−t1s . fA(cid:11)sskNugmsinagti tshfayt M(cid:11)kNis= af pk,o wfoerr so Sisi necvee reya pche rmHkutiast iaonn oPrkth, iotg foonllaolw ms athtraitx ,H anids
of 2, say M =2R, it follows that N =R. invertible. Its inverse is
cieTnth ere nlaetxito nstse ipn i(n6v) oinlv aes m eaxtprrixe sfsoirnmg .t Lheet c Aoe(cid:8)ffBi- H−1=HMT PMT (HMT=2(cid:8)IM=2)(PMT=2(cid:8)IM=2)
stand for the orthogon(cid:16)al su(cid:17)m of the matrices A (cid:1)(cid:1)(cid:1)(H2T(cid:8)IM−2):
tahnde BM, th(cid:2)aMt is,o Art(cid:8)hoBg=onaA0lB0 m. aNtorwix l edt eHfMindeedn obtye Therefore, the inverse operation is also an O(M)
HM =O(cid:8)O(cid:8)(cid:1)(cid:1)(cid:1)(cid:8)O, where the orthogonal operation.
matrix sums are applied M=2times and Ois the
Discrete Haar Series
matrix defined in (7). Then, by applying the co-
The fast Haar transform can be used for com-
efficient relations in (6) and using the fact that
f =(cid:11)R, we obtain puting partial sums of the discretized version
k k of the following Haar wavelet series in L2[0;1):
hHM[f0;f1;:::;fM−1]T X1 2Xn−1
= (cid:11)R−1;(cid:12)R−1;(cid:11)R−1;(cid:12)R−1; (8) (cid:11)0+ (cid:12)n2n=2ˆ(2nx−k):
0 0 1 1 0 k
(cid:21)
T n=0 k=0
:::;(cid:11)R−1 ;(cid:12)R−1 :
1M−1 1M−1 Let us assume, as in the previous section, that
2 2
we have a discrete signal ffjgMj=−01 associated
Twoe saoprpt ltyh ea nco eMff(cid:2)iciMentpse prrmoupteartlyio inn tmo tawtroi xg rPoMupass, iwnittehr vthael. tSiumbes vtiatluuteins gf xthje=sej =tiMmgeMj =v−0a1luoens tihnet ou n(8it)
follows:
(cid:20) (cid:21) and restricting the upper limit of n, we obtain
T
PM (cid:11)0R−1;(cid:12)R0−1;:::;(cid:11)R1−M1−1;(cid:12)R1−M1−1 RX−12Xn−1
=(cid:20)(cid:11)0R−1;:::;(cid:11)R1−M1−1;(cid:12)0R−21;:::;(cid:12)2R1−M1−1(cid:21)T: fj =(cid:11)00CM +n=0 k=0 (cid:12)nkCM2n=2ˆ(2nxj −k):
2 2
The right side of this equation is just the trans-
If we go to the next lower level, the transfor- formation H−1H applied to ffjg. The first
mations just described are repeated, only now part of this transformation, Hffjg, produces the
the matrices used are HM=2and PM=2, and they coefficients (cid:11)0;(cid:12)0;(cid:12)1;(cid:12)1;:::;(cid:12)R−1 , and the
operate only on the M=2-length vector 0 0 0 1 1M−1
wf(cid:11)a0Rv−el1e;t: :c:o;e(cid:11)fRf12i−Mc1−ie1ngtst of (cid:12)oR0b−ta2;in: :t:h;(cid:12)e R1n−Me2−x1tg leavnedl sidseu acc oesncsa dtl hep efa aorctrt,io gtrih nweah l aidcphap telain csfauftjirogen.s T tohhfae t2 H ctohne− s1ct,o anrnestpt aCrnoMt-
scaling coefficients f(cid:11)0R−2;:::;(cid:11)R1−M24−1g. These vaerec tuonr iCt vMeacntodr tsh ien vRecMt,o urss ifnCgM thp2en =s2tˆan(d2anrxdj i−nnke)gr
4
operations continue until we can no longer product. Consequently, CM = 1=M.
divide the number of components by 2. At the There are many ways of forming partial sums
R=log M step, we obtain a single scale coeffi- of discretized Haar series. The simplest ones con-
2
662 NOTICESOFTHEAMS VOLUME44, NUMBER6
walker.qxp 3/24/98 2:17 PM Page 663
Figure 3. (a)[top] Haar series partial sum, 229 terms. (b)[top] Fourier series partial sum, 229 terms.
(a)[bottom] Haar series partial sum, 92 terms. (b)[bottom] Daub4 series partial sum, 22 terms.
sist of multiplying the data by H, then setting created by the energy method, which contains
some of the resulting coefficients equal to 0, and 99.5% of the energy of this signal. This partial
then multiplying by H−1. A widely used method sum, which used 229 coefficients out of a pos-
involves specifying a threshold. All coefficients sible 8192, provides an acceptable visual repre-
whose magnitudes lie below this threshold are sentation of the signal. In fact, the sum of the
set equal to 0. This method is frequently used squares of the errors is 2:0(cid:2)10−3. By compar-
for noise removal, where coefficients whose ison a 229 coefficient Fourier series partial sum
magnitudes are significant only because of the suffers from serious drawbacks (see Figure
added noise will often lie below a well-chosen 3(b)[top]). The sum of the squares of the errors
threshold. We shall give an example of this in the is 4:5(cid:2)10−1, and there is severe oscillation and
section “Signal Denoising”. A second method a Gibbs’ effect near x=0:5. Although these lat-
keeps only the largest magnitude coefficients, ter two defects could be ameliorated using other
while setting the rest equal to 0. This method summation methods ([11], Ch. 4), there would
is convenient for making comparisons when it still be a significant deviation from the original
is known in advance how many terms are needed. signal (especially near x=0:5).
We used it in the compression example in the sec- This example illustrates how wavelet analy-
tion “Signal Compression”. A third method, sis homes right in on regions of high variability
which we shall call the energy method, involves of signals and that Fourier methods try to
specifying a fraction of the signal’s energy, where smooth them out. The size of a function’s Fourier
the energy is the square root of the sum of the coefficients is related to the frequency content
squares of the coefficients.3 We then retain the of the function, which is measured by integra-
least number of the largest magnitude coeffi- tion of the function against completely unlo-
cients whose energy exceeds this fraction of the calized basis functions. For a function having a
signal’s energy and set all other coefficients discontinuity, or some type of transient behav-
equal to 0. The energy method is useful for the- ior, this produces Fourier coefficients that de-
oretical purposes: it is clearly helpful to be able crease in magnitude at a very slow rate. Conse-
to specify in advance what fraction of the sig- quently, a large number of Fourier coefficients
nal’s energy is contained in a partial sum. We are needed to accurately represent such sig-
shall use the energy method frequently below. nals.4 Wavelet series, however, use compactly
Let us look at an example. Suppose our sig- supported basis functions which, at increasing
nal is levels of resolution, have rapidly decreasing sup-
ffj =F(j=8192)g8j=1091 ports and can zoom in on transient behavior. The
transient behaviors contribute to the magnitude
where F(x)=x1[0;:5)(x)+(x−1)1[:5;1)(x). In Fig- of only a small portion of the wavelet coefficients.
ure 3(a)[top] we show a Haar series partial sum,
4In recent years, though, significant improvements
3The Haar transform is orthogonal, so it makes sense have been achieved using local cosine bases [3, 7,
to specify energy in this way. 9, 5, 1].
JUNE/JULY1997 NOTICESOFTHEAMS 663
walker.qxp 3/24/98 2:17 PM Page 664
Figure 4. (a) Magnitudes of highest level coefficients for a function F: [top] Haar coefficients, [middle]
Daub4 coefficients; [bottom] Graph of F. (b) Sums of squares of all coefficients: [top] Haar coefficients,
[middle] Daub4 coefficients; [bottom] Graph of F. Vertical scales for highest level coefficients and sums
of squares are logarithmic.
Consequently, a small number of wavelet coeffi- DaubR4 wavelet is supportRed on [0;3)and satis-
cients are needed to accurately represent such sig- fies −11ˆ(x)dx=0 and −11xˆ(x)dx=0. Con-
nals. sequently, if the signal is constant or linear over
The Haar system performs well when the signal an interval [a;b) which contains
is constant over long stretches. This is because the [k2−n;3(k+1)2−n), then the wavelet coefficient
Haar wavelet is supported onR [0;1) and satisfies a (cid:12)n will equal 0. In Figure 4(a) we show graphs
0thorder moment condition, −11ˆ(x)dx=0. There- ofk the magnitudes of the highest level Haar co-
fore, if the signal ffjg is constant over an interval efficients and Daub4 coefficients. Each magni-
a(cid:20)xj <bsuch that [k2−n;(k+1)2−n)(cid:26)[a;b), then tude j(cid:12)1k2jis plotted at the x-coordinate k2−12
the wavelet coefficient (cid:12)nequals 0. For example, sup- for k=0;1;:::;212−1. These graphs show that
k
the highest level Haar coefficients are near 0over
pose ffj =F(j=8192)g8j=1091 where the constant parts of F, while the highest level
F(x)=(8x−1)1 (x) Daub4 coefficients are near 0over the constant
[:125;:25)
and linear parts of F. In Figure 4(b) we show
+1[:25;:75)(x)+(7−8x)1[:75;:875)(x): graphs of the sums of the squares of all the co-
efficients, which show that almost all the Daub4
In Figure 3(a)[bottom] we show a Haar series partial coefficients are near 0over the constant and lin-
sum which contains 99.5% of the energy of this sig- ear parts of F, while the Haar coefficients are
nal and uses only 92 coefficients out of a possible near 0 only over the constant parts of F. Fur-
8192. The fact that the signal is constant over three thermore, the largest magnitude Daub4 coeffi-
large subintervals of [0;1)accounts for the excellent cients are concentrated around the locations of
compression in this example. In order to obtain the points of nondifferentiability of F. This kind
wavelet bases that provide considerably more com- of local analysis illustrates one of the powerful
pression, we need a compactly supported wavelet features of wavelet analysis.
ˆ(x) which has more moments equal to zero. That Looking again at Figure 3(b)[bottom], we see
is, we want that the most serious defects of the Daub4 com-
Z pressed signal are near the points where Fis non-
1
(9) xjˆ(x)dx=0; forj =0;1;:::;L−1 differentiable. If, however, we consider the in-
−1 terval [0:4;0:6]where F is constant, the Daub4
compressed signal values differ from the values
for an integer L(cid:21)2. We say that such a wavelet has of F by no more than 1:2(cid:2)10−15 at all of the
its first L moments equal to zero. For example, a 1641 discrete values of x in this interval. In
Daub4 wavelet has its first 2 moments equal to zero. contrast, a Fourier series partial sum using 23
Using a Daub4 wavelet series for the signal above, it coefficients differs by more than 10−3 at 1441
is possible to capture 99.5% of the energy using only of these 1641 values of x. The Fourier series par-
22 coefficients! See Figure 3(b)[bottom]. This im- tial sum exhibits oscillations of amplitude
provement in compression is due to the fact that the 6:5(cid:2)10−3 around the value 1 over this subin-
664 NOTICESOFTHEAMS VOLUME44, NUMBER6
walker.qxp 3/24/98 2:18 PM Page 665
Figure 5. (a)[top] Signal. (b)[top] 37-term Daub4 approximation. (a)[bottom] 257-term Fourier cosine
series approximation. (b)[bottom] Highest level Daub4 wavelet coefficients of signal.
terval. Because the Fourier coefficients for Fare The MRA axioms tell us that (cid:30)(x) must gener-
only O(n−2), using just the first 23 coefficients ate a subspace V0 and that V0(cid:26)V1. Therefore,
produces an oscillatory approximation to Fover
X p
all of [0;1), including the subinterval [0:4;0:6]. (11) (cid:30)(x)= c 2(cid:30)(2x−k)
The highest magnitude wavelet coefficients, k
k2Z
however, are concentrated at the corner points
for F, and their terms affect only a small por- for some constants fc g. A wavelet ˆ(x), for
k
tion of the partial sum (since their basis func- which fˆ(x−k)g spans the wavelet subspace
tions are compactly supported). Consequently, W0, can be defined by5
the wavelet series provides an extremely close X p
approximation of F over the subinterval (12) ˆ(x)= (−1)kc1−k 2(cid:30)(2x−k):
[0:4;0:6]. k2Z
A major defect of the Haar wavelet is its dis-
continuity. For one thing, it is unsatisfying to use
Equations (11) and (12) generalize the equations
discontinuous functions to approximate con-
in (2) for the Haar case. The orthogonality of (cid:30)
tinuous ones. Even with discrete signals there can and ˆleads to the following equation
be undesirable jumps in Haar series partial sum
X
values (see Figure 3(a)[bottom]). Therefore, we
want to have a wavelet that is continuous. In the (13) (−1)kc1−kck=0:
next section we will describe the Daubechies k2Z
wavelets, which have their first L(cid:21)2 moments This equation, and the second equation in (14)
equal to zero and are continuous. below, imply the orthogonality of the matrices,
WN, used in the fast wavelet transform which we
Daubechies Wavelets shall discuss later in this section.
It is possible to generalize the construction of Combining (11) with the first two integrals in
the Haar wavelet so as to obtain a continuous (10), it follows that
scaling function (cid:30)(x)and a continuous wavelet
ˆ(x). Moreover, Daubechies has shown how to X p X
(14) c = 2; jc j2=1:
make them compactly supported. We will briefly k k
sketch the main ideas; more details can be found k2Z k2Z
in [4, 5, 8, 10]. Similarly, assuming that L=2, the equations in
Generalizing from the case of the Haar (9) combined with (12) imply
wavelets, we require that (cid:30)(x)and ˆ(x) satisfy
X X
Z1 Z1 (15) (−1)kck=0; k(−1)kck=0:
(cid:30)(x)dx=1; j(cid:30)(x)j2dx=1; k2Z k2Z
(10) −1 Z1 −1 5A simple proof, based on the MRA axioms, that
jˆ(x)j2dx=1: fˆ(x−k)g spans W0 can be found in [10].
−1
JUNE/JULY1997 NOTICESOFTHEAMS 665
walker.qxp 3/24/98 2:18 PM Page 666
R
And, for L>2, equations (9) and (12) yield ad- (cid:11)˜n by (cid:11)˜n = 1f(x)(cid:30)˜ (x)dx. And, we periodi-
k k 0 n;k
ditional equations similar to the ones in (15).
cally extend f with period 1, also denoting this
There is a finite set of coefficients that solves
periodic extension by f. Then, for n(cid:21)0 and
the equations in (14) and (15), namely,
k=0;1;:::;2n−1, we have
p p
Z
(16) cc02== 1344+−ppp2233;; cc13==314+4−pp2p233; (20) (cid:11)˜kn == X01fZ(xj)+21nf=(2xjX−2Zj(cid:30))2(cid:0)n2=2n(cid:30)(x(2+njx)−−kk)(cid:1)ddxx
j2Z j
with all other c =0. Using these values of c ,
k k =(cid:11)n:
the following iterative solution of (11) k
(cid:30)0(x)=1X[0;1)(xp); Similar arguments show that (cid:12)˜nk =(cid:12)nk and
(17) (cid:30)n(x)= ck 2(cid:30)n−1(2x−k); (cid:11)˜kn+2n =(cid:11)˜kn and (cid:12)˜nk+2n =(cid:12)˜nk for n(cid:21)0 and
k2Z k=0;1;:::;2n−1. After periodizing (11) and
forn(cid:21)1; (12), it follows that
X
pcoonrtveedr goens t[0o; a3 ]c. oItn tthineuno fuosl lfouwnsc tfiroonm (cid:30) (1(x2)) stuhpat- (cid:11)˜kn = cm(cid:11)˜mn++12k;
the wavelet ˆ(x) is also continuous and com- (21) mX2Z
pactly supported on [0;3]. This wavelet ˆ we (cid:12)˜nk = (−1)mc1−m(cid:11)˜mn++12k:
have been referring to as the Daub4 wavelet. m2Z
The set of coefficients fc g in (16) is the small-
k
est set of coefficients that produce a continuous
Remark. In the section “The Haar Wavelet” we
compactly supported scaling function. Other saw, for the Haar wavelet ˆ, that
sets of coefficients, related to higher values of ˆ˜ (x)=2n=2ˆ(2nx−k) for n(cid:21)0 and
L, are given in [4] and [12]. n;k
k=0;1;:::;2n−1. Almost exactly the same re-
Once the scaling function (cid:30)(x) and the
wavelet function ˆ(x) have been found, then sult holds for the Daubechies wavelets. For in-
stance, if ˆis the Daub4 wavelet, then ˆis sup-
we proceed as we did above in the Haar case. We
define the coefficients f(cid:11)ng and f(cid:12)ng by the ported on [0;3]. It follows, for n(cid:21)2, that on the
equations in (5), where nokw (cid:30) and ˆk are the unit interval ˆ˜n;0is supported on [0;3(cid:1)2−n]. On
Daubechies scaling function and wavelet, re- the unit interval, we then have
spectively. The scaling identity (11) and the ˆ˜n;k(x)=2n=2ˆ(2nx−k) for k=0;1;:::;2n−3.
wavelet definition (12) yield the following coef- Hence, for n(cid:21)2, the periodized Daub4 wavelets
ficient relations: ˆ˜n;k are identical in L2[0;1) with the wavelet
X functions 2n=2ˆ(2nx−k), except when
(cid:11)kn = cm(cid:11)mn++12k; k=2n−2and 2n−1. Similar results hold for all
(18) mX2Z the Daubechies wavelets.
(cid:12)nk = (−1)mc1−m(cid:11)mn++12k: We can discretize the series in (19) by analogy
m2Z
with the Haar series. The coefficient relations in
In order to perform calculations in L2[0;1), we (21) yield a fast wavelet transform, W, an or-
define the periodized wavelet ˆ˜ and the pe- thogonal matrix defined by
n;k
riodized scaling function (cid:30)˜ by equations (3)
and (4), only now using the nD;kaubechies wavelet W =(W2(cid:8)IM−2)
ˆ and scaling function (cid:30) in place of the Haar (cid:1)(cid:1)(cid:1)(PM=2(cid:8)IM=2)(WM=2(cid:8)IM=2)PMWM
wavelet and scaling function. Theorem 2 re-
mains valid using these periodic wavelets, but where each matrix WN is an N(cid:2)N orthogonal
the proof is more involved (see section 4.5 of [5]
matrix (as follows from (13) and the second equa-
or section 3.11 of [8]). Therefore, for each
f 2L2[0;1) we can write ttihoen i(nN (−141)))s. tTlheev mela tcroixe fWfiNciiesn utsse df (cid:11)˜toN −p1rgodauncde
f(cid:12)˜N−1g from the Nth level scaling cokefficients
X1 2Xn−1 k
(19) f(x)=(cid:11)˜00+ (cid:12)˜nkˆ˜n;k(x); f(cid:11)˜kNg as follows:
n=0 k=0 h i
T
R WN (cid:11)˜0N;:::;(cid:11)˜2NN−1
where (cid:11)˜0= 1f(x)dx and (cid:12)˜n = h i
R01f(x)ˆ˜n;k(x)d0x. W0e also define the coefficieknts = (cid:11)˜0N−1;(cid:12)˜N0−1;:::;(cid:11)˜2NN−−11−1;(cid:12)˜N2N−−11−1 T:
666 NOTICESOFTHEAMS VOLUME44, NUMBER6
walker.qxp 3/24/98 2:18 PM Page 667
Z
1
If we use the coefficients c0, c1, c2, c3 defined F(x)2R(cid:30)(2R(x−x ))dx(cid:25)F(x ):
in (16), then for N >2, WN has the following −1 k k
structure :
This approximation will hold for all period 1 con-
tinuous functions F and will be more accurate the
larger the value of 2R =M.
The higher the order of a Daubechies wavelet,
the more of its moments are zero. A Daubechies
wavelet of order 2L is defined by 2L nonzero coef-
ficients fc g, has its first Lmoments equal to zero,
k
and is supported on the interval [0;2L−1]. Gener-
ally speaking, the more moments that are zero, the
more wavelet coefficients that are nearly vanishing
for smooth functions F. This follows from consid-
ering Taylor expansions. Suppose F(x) has an L-
(cid:16) (cid:17) term Taylor expansion about the point xk=k2−n.
For N =2;W2= cc21++cc03−cc30+−cc12 . For other That is,
Daubechies wavelets, there are other finite co-
efficient sequences fckg, and the matrices WN F(x)=LX−1 1F(j)(x )(x−x )j
are defined similarly. The permutation matrix PN j! k k
j=0
is the same one that we defined for the Haar
transform and is used to sort the (N−1)stlevel + 1F(L)(tx)(x−xk)L
coefficients so that WN−1 can be applied to the L!
scaling coefficients f(cid:11)˜N−1g.
k where tx lies between x and xk. Suppose also that
As initial data for the wavelet transform we ˆis supported on [−a;a]and that ˆhas its first L
can, as we did for the Haar transform, use dis- moments equal to 0 and that jF(L)(x)j is bounded
crete data of the form ffjgMj=−01. The equations by a constant B on [(k−a)2−n;(k+a)2−n]. It then
in (15) then provide a discrete analog of the zero
follows that
moment conditions in (9) [for L = 2]; hence the
(cid:18) (cid:19)
wavelet coefficients will be 0 where the data is B a L+1=2
linear. In the last section, we saw how this can (22) j(cid:12)˜nj(cid:20) p :
k L+1=2L! 2n
produce effective compression of signals when
just the 0th and 1stmoments of ˆare 0.
This inequality shows why ˆ(x) having zero mo-
It is often the case that the initial data are val-
ues of a measured signal, i.e. ff g=fF(x )g, for ments produces a large number of small wavelet co-
k k
x =k2−n, where F is a signal obtained from a efficients. If F has some smoothness on an interval
k (c;d), then wavelet coefficients (cid:12)˜ncorresponding to
measurement process. As shown in the previous k
the basis functions ˆ(2n(x−x )) whose supports
section, we can interpret the behavior of the dis- k
are contained in (c;d)will approach 0 rapidly as n
crete case based on properties of the function
increases to 1.
F. A measured signaRl F is often described by a
convolution: F(x)= −11g(t)(cid:22)(x−t)dt, where g In addition to the Daubechies wavelets, there is
another class of compactly supported wavelets called
is the signal being measured and (cid:22)is called the
coiflets. These wavelets are also constructed using
instrument function. Such convolutions generally
the method outlined above. A coiflet of order 3L is
have greater regularity than a typical function
in L2(R). For instance, if g2L1(R) and is sup- defined by 3Lnonzero coefficients fckgand has its
first Lmoments equal to zero and is supported on
ported on a finite interval and (cid:22)=1[−r;r] for the interval [−L;2L−1]. A coiflet of order 3Lis dis-
some positive r, then F is continuous and sup-
tinguished from a Daubechies wavelet of order 2L
ported on a finite interval. By a linear change of
in that, in addition to ˆhaving its first Lmoments
variables, we may then assume that F is sup-
ported on [0;1]. The data fF(x )g then provide equal to zero, the scaling function (cid:30)for the coiflet
afipcpiernotxsi mf(cid:11)˜atRigo.n Isf fwoer atshseu hmigeh tehsatkt leFvheal ss cpaelrei ocdo e1f-, aR−l1s1o xhja(cid:30)s( xL)d−x1=m0, ofomre jnt=s 1v;a::n:i;sLhi−ng1.. FIno rp aa crotiifcluelta orf,
k order 3L, supported on [−a;a], an argument simi-
then
lar to the one that proves (22) shows that
(cid:11)˜kR =2−R=2Z−11F(x)2R(cid:30)(cid:16)2R(x−xk)(cid:17)dx j(cid:11)˜kR−CMF(xk)j(cid:20) pL+B1=2L!(cid:18)2aR(cid:19)L+1=2:
(cid:25)CMF(xk):
p
Here we have replaced 1= 2Rby the scale fac- This inequality provides a stronger theoretical jus-
tor CMand used tification for using the data fCMF(xk)g in place of
JUNE/JULY1997 NOTICESOFTHEAMS 667