f13-10

584 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Giunta,G.andMurli,A.1987, ACMTransactionsonMathematicalSoftware ,vol.13,pp.97–107. [4] Lyness, J.N. 1987, in Numerical Integration , P. Keast and G. Fairweather, eds. (Dordrecht: Reidel). [5] Pantis, G. 1975, Journal of Computational Physics , vol. 17, pp. 229–233. [6] Blakemore, M., Evans, G.A., and Hyslop, J. 1976, Journal of Computational Physics , vol. 22, pp. 352–376. [7] Lyness,J.N., andKaper,T.J.1987, SIAM JournalonScientific andStatistical Computing ,vol.8, pp. 1005–1011. [8] Thakkar,A.J., andSmith,V.H. 1975, ComputerPhysicsCommunications , vol.10,pp.73–79.[9] 13.10 Wavelet Transforms Like thefast Fouriertransform(FFT),thediscrete wavelettransform(DWT)is afast,linearoperationthatoperatesonadatavectorwhoselengthisanintegerpower of two, transforming it into a numerically different vector of the same length. AlsoliketheFFT,thewavelettransformisinvertibleandinfactorthogonal—theinverse transform, when viewed as a big matrix, is simply the transpose of the transform. Both FFT and DWT, therefore, can be viewed as a rotation in function space, from the input space (or time) domain, where the basis functions are the unit vectors e i, or Dirac delta functions in the continuumlimit, to a different domain. For the FFT,this new domain has basis functions that are the familiar sines and cosines. In the wavelet domain, the basis functions are somewhat more complicated and have the fanciful names “mother functions” and “wavelets.” Ofcoursethereareaninfinityofpossiblebasesforfunctionspace,almostallof themuninteresting! Whatmakesthewaveletbasisinterestingisthat, unlikesinesand cosines, individual wavelet functions are quite localized in space; simultaneously, likesines and cosines, individual wavelet functions are quite localized in frequency or (more precisely) characteristic scale. As we will see below, the particular kindof dual localization achieved by wavelets renders large classes of functions and operatorssparse,orsparsetosomehighaccuracy,whentransformedintothewavelet domain. Analogously with the Fourier domain, where a class of computations, likeconvolutions, become computationally fast, there is a large class of computations — those that can take advantage of sparsity — that become computationally fast in the wavelet domain [1]. Unlike sines and cosines, which define a unique Fourier transform, there is not one single unique set of wavelets; in fact, there are infinitely many possiblesets. Roughly, the different sets of wavelets make different trade-offs between how compactly they are localized in space and how smooth they are. (There are further fine distinctions.) Daubechies Wavelet FilterCoefficients A particular set of wavelets is specified by a particular set of numbers, called wavelet filter coefficients . Here, we will largely restrict ourselves to wavelet filters in a class discovered by Daubechies [2]. This class includes members ranging from highlylocalizedtohighlysmooth. Thesimplest(andmostlocalized)member,often calledDAUB4, has onlyfourcoefficients, c0,...,c 3. For the momentwe specialize to this case for ease of notation. 13.10WaveletTransforms 585Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Consider the following transformation matrix acting on a column vector of data to its right:  c 0c1c2c3 c3−c2c1−c0 c0c1c2c3 c3−c2c1−c0 ......... c0c1c2c3 c3−c2c1−c0 c2c3 c0c1 c1−c0 c3−c2 (13.10.1 ) Hereblankentriessignifyzeroes. Notethestructureofthismatrix. Thefirstrow generatesonecomponentofthedataconvolvedwiththefiltercoefficients c 0...,c 3. Likewise the third, fifth, and other odd rows. If the even rows followedthis pattern,offset by one, then the matrix would be a circulant, that is, an ordinary convolution that could be done by FFT methods. (Note how the last two rows wrap around like convolutions with periodic boundary conditions.) Instead of convolving with c 0,...,c 3,however,theevenrowsperformadifferentconvolution,withcoefficients c3,−c2,c1,−c0. The action of the matrix, overall, is thus to perform two related convolutions, then to decimate each of them by half (throw away half the values), and interleave the remaining halves. It is usefulto thinkof the filter c0,...,c 3as beinga smoothingfilter, call it H, something like a moving average of four points. Then, because of the minus signs, the filter c3,−c2,c1,−c0, call it G,i snota smoothing filter. (In signal processing contexts, HandGarecalled quadraturemirrorfilters [3].) Infact,the c’sarechosen so as to make Gyield, insofar as possible, a zeroresponse to a sufficiently smooth datavector. Thisisdonebyrequiringthesequence c3,−c2,c1,−c0tohaveacertain number of vanishing moments. When this is the case for pmoments (starting with the zeroth), a set of wavelets is said to satisfy an “approximationconditionof order p.” This results in the output of H, decimated by half, accurately representing the data’s “smooth” information. The output of G, also decimated, is referred to as the data’s “detail” information [4]. For such a characterization to be useful, it must be possible to reconstruct the original data vector of length Nfrom its N/2smooth or s-componentsand its N/2 detail or d-components. That is effected by requiring the matrix (13.10.1) to beorthogonal, so that its inverse is just the transposed matrix  c 0c3 ··· c2c1 c1−c2 ··· c3−c0 c2c1c0c3 c3−c0c1−c2 ... c2c1c0c3 c3−c0c1−c2 c2c1c0c3 c3−c0c1−c2 (13.10.2 ) 586 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Onesees immediatelythatmatrix(13.10.2)isinversetomatrix(13.10.1)ifand only if these two equations hold, c2 0+c2 1+c2 2+c2 3=1 c2c0+c3c1=0(13.10.3 ) If additionally we require the approximation condition of order p=2, then two additional relations are required, c3−c2+c1−c0=0 0c3−1c2+2c1−3c0=0(13.10.4 ) Equations (13.10.3) and (13.10.4) are 4 equations for the 4 unknowns c0,...,c 3, first recognized and solved by Daubechies. The unique solution (up to a left-right reversal) is c0=( 1+√ 3)/4√ 2 c1=( 3+√ 3)/4√ 2 c2=( 3−√ 3)/4√ 2 c3=( 1−√ 3)/4√ 2(13.10.5 ) In fact, DAUB4 is only the most compact of a sequence of wavelet sets: If we had six coefficients instead of four, there would be three orthogonality requirements in equation (13.10.3) (with offsets of zero, two, and four), and we could require the vanishingof p=3momentsinequation(13.10.4). Inthiscase,DAUB6,thesolution coefficients can also be expressed in closed form, c0=( 1+√ 10 +/radicalbig 5+2√ 10) /16√ 2 c1=( 5+√ 10 + 3/radicalbig 5+2√ 10) /16√ 2 c2=( 1 0−2√ 10 + 2/radicalbig 5+2√ 10) /16√ 2 c3=( 1 0−2√ 10−2/radicalbig 5+2√ 10) /16√ 2 c4=( 5+√ 10−3/radicalbig 5+2√ 10) /16√ 2 c5=( 1+√ 10−/radicalbig 5+2√ 10) /16√ 2 (13.10.6 ) Forhigher p,upto10,Daubechies [2]hastabulatedthecoefficientsnumerically. The number of coefficients increases by two each time pis increased by one. Discrete Wavelet Transform We have not yet defined the discrete wavelet transform (DWT), but we are almost there: The DWT consists of applying a wavelet coefficient matrix like (13.10.1) hierarchically ,firsttothefulldatavectoroflength N,thentothe“smooth” vector of length N/2, then to the “smooth-smooth” vector of length N/4, and so on until only a trivial number of “smooth- ...-smooth” components (usually 2) remain. The procedure is sometimes called a pyramidal algorithm [4], for obvious reasons. The output of the DWT consists of these remaining components and all the “detail” components that were accumulated along the way. A diagram should make the procedure clear: 13.10WaveletTransforms 587Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica). y1 y2 y3 y4 y5 y6 y7 y8 y9 y10 y11 y12 y13 y14 y15 y16 13.10.1−→ s1 d1 s2 d2 s3 d3 s4 d4 s5 d5 s6 d6 s7 d7 s8 d8 permute−→ s1 s2 s3 s4 s5 s6 s7 s8 d1 d2 d3 d4 d5 d6 d7 d8 13.10.1−→ S1 D1 S2 D2 S3 D3 S4 D4 d1 d2 d3 d4 d5 d6 d7 d8 permute−→ S1 S2 S3 S4 D1 D2 D3 D4 d1 d2 d3 d4 d5 d6 d7 d8 etc.−→ S1 S2 D1 D2 D1 D2 D3 D4 d1 d2 d3 d4 d5 d6 d7 d8  (13.10.7 ) If the length of the data vector were a higher power of two, there would be morestagesofapplying(13.10.1)(oranyotherwaveletcoefficients)andpermuting. The endpoint will always be a vector with two S’s and a hierarchy of D’s,D’s, d’s, etc. Notice that once d’s are generated, they simply propagate through to all subsequent stages. A value d iof any level is termed a “wavelet coefficient” of the original data vector;thefinalvalues S1,S2shouldstrictlybecalled“mother-functioncoefficients,” although the term “wavelet coefficients” is often used loosely for both d’s and final S’s. Since the full procedure is a composition of orthogonal linear operations, the whole DWT is itself an orthogonal linear operator. ToinverttheDWT,onesimplyreversestheprocedure,startingwiththesmallest level of the hierarchy and working (in equation 13.10.7) from right to left. The inverse matrix (13.10.2)is of course used instead of the matrix (13.10.1). Asalreadynoted,thematrices(13.10.1)and(13.10.2)embodyperiodic(“wrap- around”) boundary conditions on the data vector. One normally accepts this as a minorinconvenience: thelast few waveletcoefficientsat eachlevelof thehierarchyare affected by data from both ends of the data vector. By circularly shifting the matrix (13.10.1) N/2columns to the left, one can symmetrize the wrap-around; but this does not eliminate it. It is in fact possible to eliminate the wrap-aroundcompletely by altering the coefficients in the first and last Nrows of (13.10.1), giving an orthogonal matrix that is purely band-diagonal [5]. This variant, beyond our scope here, is useful when, e.g., the data varies by many orders of magnitude from one end of the data vector to the other. Here is a routine, wt1, that performs the pyramidal algorithm (or its inverse ifisignis negative) on some data vector a(1:n). Successive applications of the wavelet filter, and accompanying permutations, are done by an assumed routine wtstep, which must be provided. (We give examples of several different wtstep routines just below.) SUBROUTINE wt1(a,n,isign,wtstep) INTEGER isign,nREAL a(n)EXTERNAL wtstep C USES wtstep One-dimensional discrete wavelet transform. This routine implements the pyramid algo-rithm, replacing a(1:n)by its wavelet transform (for isign=1), or performing the inverse operation (for isign=-1 ). Note that nMUST be an integer power of 2. The subroutine 588 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).wtstep, whose actual name must be supplied in calling this routine, is the underlying wavelet filter. Examples of wtsteparedaub4and (preceded by pwtset)pwt. INTEGER nn if (n.lt.4) returnif (isign.ge.0) then Wavelet transform. nn=n Start at largest hierarchy, 1 if (nn.ge.4) then call wtstep(a,nn,isign) and work towards smallest. nn=nn/2 goto 1 endif else Inverse wavelet transform. nn=4 Start at smallest hierarchy, 2 if (nn.le.n) then call wtstep(a,nn,isign)nn=nn*2 and work towards largest. goto 2 endif endifreturn END Here, as a specific instanceof wtstep, is a routinefor theDAUB4 wavelets: SUBROUTINE daub4(a,n,isign) INTEGER n,isign,NMAX NMAX is the maximum allowed value of n. REAL a(n),C3,C2,C1,C0 PARAMETER (C0=0.4829629131445341,C1=0.8365163037378079, * C2=0.2241438680420134,C3=-0.1294095225512604,NMAX=1024) Applies the Daubechies 4-coefficient wavelet filter to data vector a(1:n)(forisign=1)o r applies its transpose (for isign=-1 ). Used hierarchically by routines wt1andwtn. REAL wksp(NMAX)INTEGER nh,nh1,i,j if(n.lt.4)return if(n.gt.NMAX) pause ’wksp too small in daub4’nh=n/2 nh1=nh+1 if (isign.ge.0) then Apply filter. i=1do 11j=1,n-3,2 wksp(i)=C0*a(j)+C1*a(j+1)+C2*a(j+2)+C3*a(j+3) wksp(i+nh)=C3*a(j)-C2*a(j+1)+C1*a(j+2)-C0*a(j+3)i=i+1 enddo 11 wksp(i)=C0*a(n-1)+C1*a(n)+C2*a(1)+C3*a(2) wksp(i+nh)=C3*a(n-1)-C2*a(n)+C1*a(1)-C0*a(2) else Apply transpose filter. wksp(1)=C2*a(nh)+C1*a(n)+C0*a(1)+C3*a(nh1) wksp(2)=C3*a(nh)-C0*a(n)+C1*a(1)-C2*a(nh1)j=3 do 12i=1,nh-1 wksp(j)=C2*a(i)+C1*a(i+nh)+C0*a(i+1)+C3*a(i+nh1)wksp(j+1)=C3*a(i)-C0*a(i+nh)+C1*a(i+1)-C2*a(i+nh1)j=j+2 enddo 12 endif do13i=1,n a(i)=wksp(i) enddo 13 return END 13.10WaveletTransforms 589Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).For larger sets of wavelet coefficients, the wrap-around of the last rows or columns is a programming inconvenience. An efficient implementation wouldhandle the wrap-arounds as special cases, outside of the main loop. Here, we will content ourselves with a more general scheme involving some extra arithmetic at run time. The following routine sets up any particular wavelet coefficients whosevalues you happen to know. SUBROUTINE pwtset(n) INTEGER n,NCMAX,ncof,ioff,joffPARAMETER (NCMAX=50) Maximum number of wavelet coefficients passed to pwt. REAL cc(NCMAX),cr(NCMAX) COMMON /pwtcom/ cc,cr,ncof,ioff,joff Initializing routine for pwt, here implementing the Daubechies wavelet filters with 4, 12, and 20 coefficients, as selected by the input value n. Further wavelet filters can be included in the obvious manner. This routine must be called (once) before the first use of pwt.( F o r the case n=4, the specific routine daub4is considerably faster than pwt.) INTEGER k REAL sig,c4(4),c12(12),c20(20) SAVE c4,c12,c20,/pwtcom/DATA c4/0.4829629131445341, 0.8365163037378079, * 0.2241438680420134,-0.1294095225512604/ DATA c12 /.111540743350, .494623890398, .751133908021, * .315250351709,-.226264693965,-.129766867567,* .097501605587, .027522865530,-.031582039318, * .000553842201, .004777257511,-.001077301085/ DATA c20 /.026670057901, .188176800078, .527201188932, * .688459039454, .281172343661,-.249846424327,* -.195946274377, .127369340336, .093057364604, * -.071394147166,-.029457536822, .033212674059, * .003606553567,-.010733175483, .001395351747,* .001992405295,-.000685856695,-.000116466855, * .000093588670,-.000013264203 / ncof=nsig=-1.do 11k=1,n if(n.eq.4)then cc(k)=c4(k) else if(n.eq.12)then cc(k)=c12(k) else if(n.eq.20)then cc(k)=c20(k) else pause ’unimplemented value n in pwtset’ endifcr(ncof+1-k)=sig*cc(k)sig=-sig enddo 11 ioff=-n/2 These values center the “support” of the wavelets at each level. Alternatively, the “peaks” of the wavelets can be approx- imately centered by the choices ioff=-2and joff=-n+2 . Note that daub4andpwtsetwith n=4 use different default centerings.joff=-n/2 return END Once pwtsethas been called, the following routine can be used as a specific instance of wtstep. SUBROUTINE pwt(a,n,isign) INTEGER isign,n,NMAX,NCMAX,ncof,ioff,joff PARAMETER (NMAX=2048,NCMAX=50)REAL a(n),wksp(NMAX),cc(NCMAX),cr(NCMAX) COMMON /pwtcom/ cc,cr,ncof,ioff,joff 590 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Partial wavelet transform: applies an arbitrary wavelet filter to data vector a(1:n)(for isign=1) or applies its transpose (for isign=-1 ). Used hierarchically by routines wt1 andwtn. The actual filter is determined by a preceding (and required) call to pwtset, which initializes the common block pwtcom. INTEGER i,ii,j,jf,jr,k,n1,ni,nj,nh,nmod REAL ai,ai1 SAVE /pwtcom/if (n.lt.4) returnnmod=ncof*n A positive constant equal to zero mod n. n1=n-1 Mask of all bits, since na power of 2. nh=n/2do 11j=1,n wksp(j)=0. enddo 11 if (isign.ge.0) then Apply filter. ii=1 do13i=1,n,2 ni=i+nmod+ioff Pointer to be incremented and wrapped-around. nj=i+nmod+joffdo 12k=1,ncof jf=iand(n1,ni+k) We use bitwise and to wrap-around the pointers. jr=iand(n1,nj+k)wksp(ii)=wksp(ii)+cc(k)*a(jf+1) wksp(ii+nh)=wksp(ii+nh)+cr(k)*a(jr+1) enddo 12 ii=ii+1 enddo 13 else Apply transpose filter. ii=1do 15i=1,n,2 ai=a(ii) ai1=a(ii+nh)ni=i+nmod+ioff See comments above. nj=i+nmod+joff do 14k=1,ncof jf=iand(n1,ni+k)+1jr=iand(n1,nj+k)+1wksp(jf)=wksp(jf)+cc(k)*ai wksp(jr)=wksp(jr)+cr(k)*ai1 enddo 14 ii=ii+1 enddo 15 endif do16j=1,n Copy the results back from workspace. a(j)=wksp(j) enddo 16 return END What Do Wavelets LookLike? We are now in a position actually to see some wavelets. To do so, we simply run unit vectors through any of the above discrete wavelet transforms, with isign negative so that the inverse transform is performed. Figure 13.10.1 shows the DAUB4 wavelet that is the inverse DWT of a unit vector in the 5th componentof avector of length 1024, and also the DAUB20 wavelet that is the inverse of the 22nd component. (One needs to go to a later hierarchical level for DAUB20, to avoid a wavelet with a wrapped-around tail.) Other unit vectors would give wavelets with the same shapes, but different positions and scales. 13.10WaveletTransforms 591Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).0 100 200 300 400 500 600 700 800 900 1000 0 100 200 300 400 500 600 700 800 900 1000−.1−.050.05.1 −.1−.050.05.1 DAUB20 e 22DAUB4 e 5 Figure 13.10.1. Wavelet functions, that is, single basis functions from the wavelet families DAUB4 and DAUB20. A complete, orthonormal wavelet basis consists of scalings and translations of either one of these functions. DAUB4 has an in finite number of cusps; DAUB20 would show similar behavior in a higher derivative. One sees that both DAUB4 and DAUB20 have wavelets that are continuous. DAUB20waveletsalsohavehighercontinuousderivatives. DAUB4hasthepeculiarpropertythatitsderivativeexistsonly almosteverywhere. Examplesofwhereitfails to exist are the points p/2 n, where pandnare integers; at such points, DAUB4 is left differentiable,but not right differentiable! This kind of discontinuity —at least in some derivative —is a necessary feature of wavelets with compact support, like the Daubechies series. For every increase in the number of wavelet coef ficients by two, the Daubechies wavelets gain about halfa derivative of continuity. (But not exactly half; the actual orders of regularity are irrational numbers!) Note that the fact that wavelets are not smooth does not prevent their having exactrepresentationsforsomesmoothfunctions,asdemandedbytheirapproximation order p. The continuity of a wavelet is not the same as the continuity of functions thataset ofwaveletscanrepresent. Forexample,DAUB4canrepresent(piecewise) linear functions of arbitrary slope: in the correct linear combinations, the cusps all cancel out. Every increase of two in the number of coef ficients allows one higher order of polynomial to be exactly represented. Figure 13.10.2 shows the result of performing the inverse DWT on the input vectore10+e58, again for the two different particularwavelets. Since 10 lies early in the hierarchical range of 9−16, that wavelet lies on the left side of the picture. Since58liesinalater(smaller-scale)hierarchy,itisanarrowerwavelet;intherange of33–64it is towards the end, so it lies on the right side of the picture. Note that smaller-scale wavelets are taller, so as to have the same squared integral. 592 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).0 100 200 300 400 500 600 700 800 900 1000 0 100 200 300 400 500 600 700 800 900 1000−.20.2 DAUB4 e 10 + e58 −.20.2 Lemarie e 10 + e58 Figure 13.10.2. More wavelets, here generated from the sum of two unit vectors, e10+e58, which are in different hierarchical levels of scale, and also at different spatial positions. DAUB4 wavelets (a) are defined by a filter in coordinate space (equation 13.10.5), while Lemarie wavelets (b) are de fined by afilter most easily written in Fourier space (equation 13.10.14). Wavelet Filtersinthe FourierDomain The Fourier transform of a set of filter coefficients cjis given by H(ω)=/summationdisplay jcjeijω(13.10.8 ) Here His a function periodic in 2π, and it has the same meaning as before: It is the wavelet filter, now written in the Fourier domain. A very useful fact is that the orthogonality conditions for the c’s (e.g., equation 13.10.3 above) collapse to two simple relations in the Fourier domain, 1 2|H(0)|2=1 ( 13.10.9 ) and 1 2/bracketleftbig |H(ω)|2+|H(ω+π)|2/bracketrightbig =1 ( 13.10.10 ) Likewise the approximation condition of order p(e.g., equation 13.10.4 above) has a simple formulation, requiring that H(ω)have a pth order zero at ω=π, or (equivalently) H(m)(π)=0 m=0,1,...,p −1( 13.10.11 ) 13.10WaveletTransforms 593Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).ItisthusrelativelystraightforwardtoinventwaveletsetsintheFourierdomain. You simply invent a function H(ω)satisfying equations (13.10.9) –(13.10.11). To find the actual cj’s applicable to a data (or s-component) vector of length N, and withperiodicwrap-aroundasinmatrices(13.10.1)and(13.10.2),youinvertequation (13.10.8) by the discrete Fourier transform cj=1 NN−1/summationdisplay k=0H(2πk N)e−2πijk/N(13.10.12 ) The quadrature mirror filterG(reversed cj’s with alternating signs), incidentally, has the Fourier representation G(ω)=e−iωH*(ω+π)( 13.10.13 ) where asterisk denotes complex conjugation. In general the above procedure will notproduce wavelet filters with compact support. In other words, all Nof the cj’s,j=0,...,N −1will in general be nonzero (though they may be rapidly decreasing in magnitude). The Daubechies wavelets,orotherwaveletswithcompactsupport,arespeciallychosensothat H(ω) is a trigonometric polynomial with only a small number of Fourier components,guaranteeing that there will be only a small number of nonzero c j’s. On the other hand,there is sometimes no particularreason to demandcompact support. Giving it up in fact allows the ready construction of relatively smoother wavelets (higher values of p). Even without compact support, the convolutions implicit in the matrix (13.10.1)can be done ef ficiently by FFT methods. Lemarie’s wavelet (see [4]) has p=4, does not have compact support, and is defined by the choice of H(ω), H(ω)=/bracketleftbigg 2(1−u)4315−420u+ 126 u2−4u3 315−420v+ 126 v2−4v3/bracketrightbigg1/2 (13.10.14 ) where u≡sin2ω 2v≡sin2ω (13.10.15 ) It is beyond our scope to explain where equation (13.10.14) comes from. An informaldescriptionis thatthequadraturemirror filterG(ω)derivingfromequation (13.10.14)hasthepropertythatitgivesidenticallyzerowhenappliedtoanyfunction whose odd-numbered samples are equal to the cubic spline interpolation of its even-numbered samples. Since this class of functions includes many very smoothmembers,itfollowsthat H(ω)doesagoodjoboftrulyselectingafunction ’ssmooth informationcontent. Sample Lemarie wavelets are shown in Figure 13.10.2. 594 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).wavelet amplitude.0001 10−710−610−5.001.01.1110 0 100 200 300 400 500 600 700 800 900 10000 100 200 300 400 500 600 700 800 90010000.511.5 wavelet number Figure 13.10.3. (a) Arbitrary test function, with cusp, sampled on a vector of length 1024. (b) Absolute value of the 1024 wavelet coef ficients produced by the discrete wavelet transform of (a). Note log scale. The dotted curve plots the same amplitudes when sorted by decreasing size. One sees that only 130 out of 1024 coef ficients are larger than 10−4(or larger than about 10−5times the largest coefficient, whose value is ∼10). Truncated Wavelet Approximations Most of the usefulness of wavelets rests on the fact that wavelet transforms can usefully be severely truncated, that is, turned into sparse expansions. The case of Fourier transforms is different: FFTs are ordinarily used without truncation, to compute fast convolutions, for example. This works because the convolutionoperator is particularly simple in the Fourier basis. There are not, however, any standardmathematicaloperationsthat are especially simple in the wavelet basis. To see how truncation works, consider the simple example shown in Figure 13.10.3. The upper panel shows an arbitrarily chosen test function, smooth except for a square-root cusp, sampled onto a vector of length 2 10. The bottom panel (solid curve) shows, on a log scale, the absolute value of the vector ’s components after it has been run through the DAUB4 discrete wavelet transform. One notes, from right to left, the different levels of hierarchy, 513 –1024, 257 –512, 129–256, etc. Withineachlevel,thewaveletcoef ficientsarenon-negligibleonlyverynearthe location of the cusp, or very near the left and right boundaries of the hierarchical range (edge effects). ThedottedcurveinthelowerpanelofFigure13.10.3plotsthesameamplitudes as the solid curve, but sorted into decreasing order of size. One can read off, forexample, that the 130th largest wavelet coef ficient has an amplitude less than 10 −5 of the largest coef ficient, whose magnitude is ∼10(power or square integral ratio less than 10−10). Thus, the example function can be represented quite accurately by only 130, rather than 1024, coef ficients—the remaining ones being set to zero. Note that this kind of truncation makes the vector sparse, but not shorter than 1024. It is veryimportant that vectors in wavelet space be truncatedaccording to theamplitude of the components, not their position in the vector. Keeping the 13.10WaveletTransforms 595Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).first 256 components of the vector (all levels of the hierarchy except the last two) would give an extremely poor, and jagged, approximation to the function. Whenyou compress a function with wavelets, you have to record both the values and the positions of the nonzero coef ficients. Generally, compact (and therefore unsmooth) wavelets are better for lower accuracy approximation and for functions with discontinuities (like edges), while smooth(andthereforenoncompact)waveletsarebetterforachievinghighnumerical accuracy. This makes compact wavelets a good choice for image compression, for example,whileitmakessmoothwaveletsbestforfastsolutionofintegralequations. Wavelet TransforminMultidimensions A wavelet transform of a d-dimensional array is most easily obtained by transformingthearraysequentiallyonits firstindex(forallvaluesofitsotherindices), then on its second, and so on. Each transformation corresponds to multiplication by an orthogonal matrix. By matrix associativity, the result is independent of theorder in which the indices were transformed. The situation is exactly like that for multidimensionalFFTs. AroutineforeffectingthemultidimensionalDWTcanthus be modeled on a multidimensional FFT routine like fourn: SUBROUTINE wtn(a,nn,ndim,isign,wtstep) INTEGER isign,ndim,nn(ndim),NMAX REAL a(*) EXTERNAL wtstepPARAMETER (NMAX=1024) C USES wtstep Replaces aby its ndim-dimensional discrete wavelet transform, if isignis input as 1. nn is aninteger array oflength ndim, containing the lengths of eachdimension (number ofreal values), which MUST all be powers of 2. ais a real array of length equal to the product of these lengths, in which the data are stored as in a multidimensional real FORTRAN array. Ifisignis input as −1,ais replaced by its inverse wavelet transform. The subroutine wtstep, whose actual name must be supplied in calling this routine, is the underlying wavelet filter. Examples of wtsteparedaub4and (preceded by pwtset)pwt. INTEGER i1,i2,i3,idim,k,n,nnew,nprev,nt,ntotREAL wksp(NMAX)ntot=1 do 11idim=1,ndim ntot=ntot*nn(idim) enddo 11 nprev=1 do16idim=1,ndim Main loop over the dimensions. n=nn(idim)nnew=n*nprev if (n.gt.4) then do 15i2=0,ntot-1,nnew do14i1=1,nprev i3=i1+i2 do12k=1,n Copy the relevant row or column or etc. into workspace. wksp(k)=a(i3) i3=i3+nprev enddo 12 if (isign.ge.0) then Do one-dimensional wavelet transform. nt=n 1 if (nt.ge.4) then call wtstep(wksp,nt,isign) nt=nt/2goto 1 endif 596 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).else Or inverse transform. nt=4 2 if (nt.le.n) then call wtstep(wksp,nt,isign)nt=nt*2 goto 2 endif endifi3=i1+i2 do 13k=1,n Copy back from workspace. a(i3)=wksp(k)i3=i3+nprev enddo 13 enddo 14 enddo 15 endif nprev=nnew enddo 16 returnEND Here, as before, wtstepis an individualwavelet step, either daub4orpwt. Compressionof Images An immediate application of the multidimensional transform wtnis to image compression. The overall procedure is to take the wavelet transform of a digitized image, and then to “allocate bits ”among the wavelet coef ficients in some highly nonuniform,optimized,manner. In general,largewaveletcoef ficientsget quantized accurately, while small coef ficients are quantized coarsely with only a bit or two —or else are truncated completely. If the resulting quantization levels are still statistically nonuniform, they may then be further compressed by a technique like Huffman coding ( §20.4). Whileamoredetaileddescriptionofthe “backend”ofthisprocess,namelythe quantizationandcodingoftheimage,isbeyondourscope,itisquitestraightforward to demonstratethe “front-end”wavelet encodingwith a simpletruncation: We keep (withfullaccuracy)allwaveletcoef ficientslargerthansomethreshold,andwedelete (set to zero) all smaller wavelet coef ficients. We can then adjust the threshold to vary the fraction of preserved coef ficients. Figure13.10.4showsasequenceofimagesthatdifferinthenumberofwavelet coefficients that have been kept. The original picture (a), which is an of ficial IEEE test image, has 256 by 256 pixels with an 8-bit grayscale. The two reproductions following are reconstructed with 23% (b) and 5.5% (c) of the 65536 wavelet coefficients. The latter image illustrates the kind of compromises made by the truncated wavelet representation. High-contrast edges (the model ’s right cheek and hairhighlights,e.g.) aremaintainedatarelativelyhighresolution,whilelow-contrastareas (the model ’s left eye and cheek, e.g.) are washed out into what amounts to large constant pixels. Figure 13.10.4 (d) is the result of performing the identical procedure with Fourier, instead of wavelet, transforms: The figure is reconstructed from the 5.5% of 65536 real Fourier components having the largest magnitudes. One sees that, since sines and cosines are nonlocal,the resolutionis uniformlypoor acrossthepicture;also,thedeletionofanycomponentsproducesamottled “ringing” everywhere. (Practical Fourier image compression schemes therefore break up an 13.10WaveletTransforms 597Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Figure13.10.4. (a)IEEEtestimage, 256×256pixelswith8-bitgrayscale. (b)Theimageistransformed into thewavelet basis; 77% of its wavelet components are set to zero (those of smallest magnitude); it is then reconstructed from the remaining 23%. (c) Same as (b), but 94.5% of the wavelet components are deleted. (d) Same as (c), but the Fourier transform is used instead of the wavelet transform. Wavelet coefficients are better than the Fourier coef ficients at preserving relevant details. image into small blocks of pixels, 16×16, say, and do rather elaborate smoothing across block boundaries when the image is reconstructed.) Fast SolutionofLinear Systems One of the most interesting, and promising, wavelet applications is linear algebra. Thebasicidea [1]istothinkofanintegraloperator(thatis,alargematrix)as adigitalimage. Supposethattheoperatorcompresseswellunderatwo-dimensional wavelet transform, i.e., that a large fraction of its wavelet coef ficients are so small as to be negligible. Thenanylinearsystem involvingtheoperatorbecomesa sparse 598 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).system in the wavelet basis. In other words, to solve A·x=b (13.10.16 ) wefirst wavelet-transform the operator Aand the right-hand side bby /tildewideA≡W·A·WT,/tildewideb≡W·b (13.10.17 ) whereWrepresents the one-dimensional wavelet transform, then solve /tildewideA·/tildewidex=/tildewideb (13.10.18 ) andfinally transform to the answer by the inverse wavelet transform x=WT·/tildewidex (13.10.19 ) (Note that the routine wtndoes the complete transformation of Ainto/tildewideA.) A typical integral operator that compresses well into wavelets has arbitrary (or even nearlysingular) elements near to its main diagonal,but becomes smoothaway from the diagonal. An example might be Aij=/braceleftbigg−1 ifi=j |i−j|−1/2otherwise(13.10.20 ) Figure13.10.5showsagraphicalrepresentationofthewavelettransformofthis matrix, where iandjrange over 1...256, using the DAUB12 wavelets. Elements larger in magnitude than 10−3times the maximum element are shown as black pixels, while elements between 10−3and 10−6are shown in gray. White pixels are <10−6. The indices iandjeach number from the lower left. In thefigure, one sees the hierarchical decomposition into power-of-two sized blocks. At the edges or corners of the various blocks, one sees edge effects caused by the wrap-around wavelet boundary conditions. Apart from edge effects, within each block, the nonnegligible elements are concentrated along the block diagonals.This is a statement that, for this type of linear operator,a wavelet is coupledmainly to near neighbors in its own hierarchy (square blocks along the main diagonal) and near neighbors in other hierarchies (rectangular blocks off the diagonal). The number of nonnegligible elements in a matrix like that in Figure 13.10.5 scales only as N, the linear size of the matrix; as a rough rule of thumb it is about 10Nlog 10(1//epsilon1), where /epsilon1is the truncation level, e.g., 10−6. For a 2000 by 2000 matrix, then, the matrix is sparse by a factor on the order of 30. Various numerical schemes can be used to solve sparse linear systems of this “hierarchically band diagonal ”form. Beylkin, Coifman, and Rokhlin [1]make the interesting observations that (1) the product of two such matrices is itself hierarchically band diagonal (truncating, of course, newly generated elements thatare smaller than the predetermined threshold /epsilon1); and moreover that (2) the product can be formed in order Noperations. Fast matrix multiplication makes it possible to find the matrix inverse by Schultz’s (or Hotelling ’s) method, see §2.5. 13.10WaveletTransforms 599Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Figure 13.10.5. Wavelet transform of a 256×256matrix, represented graphically. The original matrix has a discontinuous cusp along its diagonal, decaying smoothly away on both sides of the diagonal. In wavelet basis,thematrix becomessparse: Components largerthan 10−3areshownasblack, components larger than 10−6as gray, and smaller-magnitude components are white. The matrix indices iand j number from the lower left. Otherschemesarealsopossibleforfastsolutionofhierarchicallybanddiagonal forms. For example, one can use the conjugate gradient method, implemented in§2.7 as linbcg. CITED REFERENCES AND FURTHER READING: Daubechies, I. 1992, Wavelets (Philadelphia: S.I.A.M.). Strang, G. 1989, SIAM Review , vol. 31, pp. 614–627. Beylkin, G., Coifman, R., and Rokhlin, V. 1991, Communications on Pure and Applied Mathe- matics, vol. 44, pp. 141–183. [1] Daubechies, I. 1988, Communications on Pureand AppliedMathematics , vol. 41, pp. 909–996. [2] Vaidyanathan, P.P. 1990, Proceedings of the IEEE , vol. 78, pp. 56–93. [3] Mallat, S.G. 1989, IEEE Transactions on Pattern Analysis and Machine Intelligence , vol. 11, pp. 674–693. [4] Freedman, M.H., and Press, W.H. 1992, preprint. [5] 600 Chapter13. FourierandSpectralApplicationsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).13.11 Numerical Use of the SamplingTheorem In§6.10 we implemented an approximating formula for Dawson ’s integral due to Rybicki. Now that we have become Fourier sophisticates, we can learn that the formuladerives from numerical application of the sampling theorem ( §12.1), normally considered to be a purely analytic tool. Our discussion is identical to Rybicki [1]. For present purposes, the sampling theorem is most conveniently stated as follows: Consider an arbitrary function g(t)and the grid of sampling points tn=α+nh, where n ranges over the integers and αis a constant that allows an arbitrary shift of the sampling grid. We then write g(t)=∞/summationdisplay n=−∞g(tn)s i n cπ h(t−tn)+e(t)( 13.11.1 ) where sincx≡sinx/x. The summation over the sampling points is called the sampling representation ofg(t), and e(t)is its error term. The sampling theorem asserts that the sampling representation is exact, that is, e(t)≡0, if the Fourier transform of g(t), G(ω)=/integraldisplay∞ −∞g(t)eiωtdt (13.11.2 ) vanishes identically for |ω|≥π/h. Whencan sampling representations be used toadvantage fortheapproximate numerical computation of functions? In order that the error term be small, the Fourier transform G(ω) must be suf ficiently small for |ω|≥π/h. On the other hand, in order for the summation in (13.11.1) to be approximated by a reasonably small number of terms, the function g(t) itself should be very small outside of a fairly limited range of values of t. Thus we are led to two conditions to be satis fied in order that (13.11.1) be useful numerically: Both the function g(t)and its Fourier transform G(ω)must rapidly approach zero for large values of their respective arguments. Unfortunately, thesetwoconditions aremutuallyantagonistic —theUncertaintyPrinci- pleinquantum mechanics. Thereexiststrictlimitsonhow rapidlythesimultaneous approach to zero can be in both arguments. According to a theorem of Hardy [2],i fg(t)=O(e−t2) as|t|→∞andG(ω)=O(e−ω2/4)as|ω|→∞, then g(t)≡Ce−t2, where Cis a constant. This can be interpreted as saying that of all functions the Gaussian is the mostrapidly decaying in both tandω, and in this sense is the “best”function to be expressed numerically as a sampling representation. Let us then write for the Gaussian g(t)=e −t2, e−t2=∞/summationdisplay n=−∞e−t2 nsincπ h(t−tn)+e(t)( 13.11.3 ) The error e(t)depends on the parameters handαas well as on t, but it is suf ficient for the present purposes to state the bound, |e(t)|<e−(π/2h)2(13.11.4 ) which can be understood simply as the order of magnitude of the Fourier transform of the Gaussian at the point where it “spills over ”into the region |ω|>π / h. When the summation in (13.11.3) is approximated by one with finite limits, say from N0−NtoN0+N, where N0is the integer nearest to −α/h, there is a further truncation error. However, if Nis chosen so that N>π / (2h2), the truncation error in the summation is less than the bound given by (13.11.4), and, since this bound is an overestimate, weshall continue to use it for (13.11.3) as well. The truncated summation gives a remarkablyaccurate representation for the Gaussian even for moderate values of N. For example, |e(t)|<5×10 −5forh=1/2andN=7;|e(t)|<2×10−10forh=1/3andN=1 5; and|e(t)|<7×10−18forh=1/4andN=2 5.