13.10 Wavelet Transforms 591 The splitting point b must be chosen large enough that the remaining integral over(b,)is small.Successive terms in its asymptotic expansion are found by integrating by parts.The integral over(a,b)can be done using dftint.You keep as many terms in the asymptotic expansion as you can easily compute.See[6]for some examples of this idea.More powerful methods,which work well for long-tailed functions but which do not use the FFT, are described in [7-9]. CITED REFERENCES AND FURTHER READING: Stoer,J.,and Bulirsch,R.1980,Introduction to Numerical Analysis(New York:Springer-Verlag), p.88.[1] Narasimhan,M.S.and Karthikeyan,M.1984,IEEE Transactions on Antennas Propagation, vol.32,pp.404-408.[2] Filon,L.N.G.1928.Proceedings of the Royal Society of Edinburgh,vol.49,pp.38-47.[3] Giunta,G.and Murli,A.1987,ACM Transactions on Mathematical Software,vol.13,pp.97-107. 4] ICAL 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] RECIPES Blakemore,M.,Evans,G.A.,and Hyslop,J.1976,Journal of Computational Physics,vol.22. pp.352-376.7] 9 Lyness,J.N.,and Kaper,T.J.1987,SIAM Journal on Scientific and Statistical Computing,vol.8. pp.1005-1011.[8] Thakkar,A.J.,and Smith,V.H.1975,Computer Physics Communications,vol.10,pp.73-79.[9] 邑白 9 13.10 Wavelet Transforms Like the fast Fourier transform(FFT),the discrete wavelet transform(DWT)is a fast,linear operation that operates on a data vector whose length is an integer power of two,transforming it into a numerically different vector of the same length.Also like the FFT,the wavelet transform is invertible and in fact orthogonal-the inverse 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 Numerical Recipes 10621 43106 the input space(or time)domain,where the basis functions are the unit vectors ei, or Dirac delta functions in the continuum limit,to a different domain.For the FFT. this new domain has basis functions that are the familiar sines and cosines.In the (outside wavelet domain,the basis functions are somewhat more complicated and have the North fanciful names“mother functions'”and“wavelets. Of course there are an infinity of possible bases for function space,almost all of them uninteresting!What makes the wavelet basis interesting is that,unlike sines and cosines,individual wavelet functions are quite localized in space;simultaneously, like sines and cosines,individual wavelet functions are quite localized in frequency or (more precisely)characteristic scale.As we will see below,the particular kind of dual localization achieved by wavelets renders large classes of functions and operators sparse,or sparse to some high accuracy,when transformed into the wavelet domain.Analogously with the Fourier domain,where a class of computations,like convolutions,become computationally fast,there is a large class of computations
13.10 Wavelet Transforms 591 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). The splitting point b must be chosen large enough that the remaining integral over (b,∞) is small. Successive terms in its asymptotic expansion are found by integrating by parts. The integral over (a, b) can be done using dftint. You keep as many terms in the asymptotic expansion as you can easily compute. See [6] for some examples of this idea. More powerful methods, which work well for long-tailed functions but which do not use the FFT, are described in [7-9]. CITED REFERENCES AND FURTHER READING: Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag), p. 88. [1] Narasimhan, M.S. and Karthikeyan, M. 1984, IEEE Transactions on Antennas & Propagation, vol. 32, pp. 404–408. [2] Filon, L.N.G. 1928, Proceedings of the Royal Society of Edinburgh, vol. 49, pp. 38–47. [3] Giunta, G. and Murli, A. 1987, ACM Transactions on Mathematical Software, 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., and Kaper, T.J. 1987, SIAM Journal on Scientific and Statistical Computing, vol. 8, pp. 1005–1011. [8] Thakkar, A.J., and Smith, V.H. 1975, Computer Physics Communications, vol. 10, pp. 73–79. [9] 13.10 Wavelet Transforms Like the fast Fourier transform (FFT), the discrete wavelet transform (DWT) is a fast, linear operation that operates on a data vector whose length is an integer power of two, transforming it into a numerically different vector of the same length. Also like the FFT, the wavelet transform is invertible and in fact orthogonal — the inverse 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 continuum limit, 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.” Of course there are an infinity of possible bases for function space, almost all of them uninteresting! What makes the wavelet basis interesting is that, unlike sines and cosines, individual wavelet functions are quite localized in space; simultaneously, like sines and cosines, individual wavelet functions are quite localized in frequency or (more precisely) characteristic scale. As we will see below, the particular kind of dual localization achieved by wavelets renders large classes of functions and operators sparse, or sparse to some high accuracy, when transformed into the wavelet domain. Analogously with the Fourier domain, where a class of computations, like convolutions, become computationally fast, there is a large class of computations
592 Chapter 13.Fourier and Spectral Applications -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 possible sets.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 Filter Coefficients 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 [21.This class includes members ranging from highly localized to highly smooth.The simplest (and most localized)member,often 茶 called DAUB4,has only four coefficients,co,...,c3.For the moment we specialize to this case for ease of notation. Consider the following transformation matrix acting on a column vector of RECIPES I data to its right: 。●0 令 「co C1 C2 C3 Press. C3 -C2C1-C0 Co C1 C2 C3 C3 -C2C1 -C0 (13.10.1) SCIENTIFIC Co C1 C2 C3 C3 -C2 C1 -C0 C3 Co C1 LC1 -C0 C3 -c2」 Here blank entries signify zeroes.Note the structure of this matrix.The first row generates one component of the data convolved with the filter coefficients co...,c3. 10.621 Likewise the third,fifth,and other odd rows.If the even rows followed this pattern, Numerica offset by one,then the matrix would be a circulant,that is,an ordinary convolution 43106 that could be done by FFT methods.(Note how the last two rows wrap around Recipes like convolutions with periodic boundary conditions.Instead of convolving with (outside co,...,c3,however,the even rows perform a different convolution,with coefficients c3,-c2,c1,-co.The action of the matrix,overall,is thus to perform two related North convolutions,then to decimate each of them by half(throw away half the values), and interleave the remaining halves. It is useful to think of the filter co,...,ca as being a smoothing filter,call it H, something like a moving average of four points.Then,because of the minus signs, the filter c3,-c2,c1,-co,call it G,is not a smoothing filter.(In signal processing contexts,H and G are called quadrature mirror filters [3].)In fact,the c's are chosen so as to make G yield,insofar as possible,a zero response to a sufficiently smooth data vector.This is done by requiring the sequence c3,-c2,c1,-co to have a certain number of vanishing moments.When this is the case for p moments(starting with the zeroth),a set of wavelets is said to satisfy an"approximation condition of order
592 Chapter 13. Fourier and Spectral Applications Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). — 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 possible sets. 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 Filter Coefficients 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 highly localized to highly smooth. The simplest (and most localized) member, often called DAUB4, has only four coefficients, c0,...,c3. For the moment we specialize to this case for ease of notation. Consider the following transformation matrix acting on a column vector of data to its right: c0 c1 c2 c3 c3 −c2 c1 −c0 c0 c1 c2 c3 c3 −c2 c1 −c0 . . . . . . ... c0 c1 c2 c3 c3 −c2 c1 −c0 c2 c3 c0 c1 c1 −c0 c3 −c2 (13.10.1) Here blank entries signify zeroes. Note the structure of this matrix. The first row generates one component of the data convolved with the filter coefficients c 0 ...,c3. Likewise the third, fifth, and other odd rows. If the even rows followed this 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 c0,...,c3, however, the even rows perform a different convolution, with coefficients 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 useful to think of the filter c0,...,c3 as being a smoothing filter, 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, is not a smoothing filter. (In signal processing contexts, H and G are called quadrature mirror filters [3].) In fact, the c’s are chosen so as to make G yield, insofar as possible, a zero response to a sufficiently smooth data vector. This is done by requiring the sequence c 3, −c2, c1, −c0 to have a certain number of vanishing moments. When this is the case for p moments (starting with the zeroth), a set of wavelets is said to satisfy an “approximation condition of order
13.10 Wavelet Transforms 593 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[4l. For such a characterization to be useful,it must be possible to reconstruct the original data vector of length N from its N/2 smooth or s-components and its N/2 detail or d-components.That is effected by requiring the matrix(13.10.1)to be orthogonal,so that its inverse is just the transposed matrix Co C3 C2 C1 C1 -C2 C3 -C0 C2 C1 Co C3 C3 -C0C1 -C2 (13.10.2) C2 C1 Co C3 C3 -C0C1-C2 C2 C1 C0 C3 RECIPES C3 -c0 C1 -C2 9 One sees immediately that matrix (13.10.2)is inverse to matrix (13.10.1)if and only if these two equations hold, 哈+c+c吃+c=1 (13.10.3) C2c0+c3C1=0 If additionally we require the approximation condition of order p =2,then two IENTIFIC additional relations are required, C3-C2+C1-c0=0 (13.10.4) 0c3-1c2+2c1-3c0=0 Equations (13.10.3)and (13.10.4)are 4 equations for the 4 unknowns co,...,c3, first recognized and solved by Daubechies.The unique solution(up to a left-right Sg② Numerica 10.621 reversal)is 0=(1+V3)/4v2 c1=(3+V3)/4v② 4310 (13.10.5) Recipes c2=(3-V/42 c3=(1-V3)/4V2 (outside 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 vanishing ofp =3 moments in equation(13.10.4).In this case,DAUB6,the solution coefficients can also be expressed in closed form, c0=(1+√o+V5+2Vo/16v2c1=(5+V0+3V5+2W1ō)/16v2 c2=(10-2W0+2V5+2yo)/16√2c3=(10-2面-2V5+2Wō/16√2 c4=(6+V10-3V5+2W1o)/16v2c=(1+√10-V5+2W1o)/16v2 (13.10.6)
13.10 Wavelet Transforms 593 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 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 N from its N/2 smooth or s-components and its N/2 detail or d-components. That is effected by requiring the matrix (13.10.1) to be orthogonal, so that its inverse is just the transposed matrix c0 c3 ··· c2 c1 c1 −c2 ··· c3 −c0 c2 c1 c0 c3 c3 −c0 c1 −c2 ... c2 c1 c0 c3 c3 −c0 c1 −c2 c2 c1 c0 c3 c3 −c0 c1 −c2 (13.10.2) One sees immediately that matrix (13.10.2) is inverse to matrix (13.10.1) if and 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 c 0,...,c3, 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 vanishing of p = 3 moments in equation (13.10.4). In this case, DAUB6, the solution coefficients can also be expressed in closed form, c0 = (1 + √10 + 5+2√10)/16√2 c1 = (5 + √10 + 35+2√10)/16√2 c2 = (10 − 2 √10 + 2 5+2√10)/16√2 c3 = (10 − 2 √10 − 2 5+2√10)/16√2 c4 = (5 + √10 − 3 5+2√10)/16√2 c5 = (1 + √10 − 5+2√10)/16√2 (13.10.6)
594 Chapter 13.Fourier and Spectral Applications For higher p,up to 10,Daubechies [2]has tabulated the coefficients numerically.The number of coefficients increases by two each time p is 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,first to the full data vector of length N,then to the "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: 100 17 S1 S17 S17 d 旺山 S2 S etc. 4 4 13.10.1 permute D2 (North America users to make one paper UnN电.t THE ye ds D 2 7 8 ART 8 13.10.1 da permute D 9 S5 器 only),or y10 ds 出 Programs 11 s 12 de d 91 ST s y14 dr 2 Copyright (C) y15 S8 dz d to dir y16 ds ds ds ds (13.10.7 ectcustser If the length of the data vector were a higher power of two.there would be OF SCIENTIFIC COMPUTING(ISBN 0-521- more stages ofapplying(13.10.1)(or any other wavelet coefficients)and permuting. The endpoint will always be a vector with two S's and a hierarchy of D's,D's, v@cam d's,etc.Notice that once d's are generated,they simply propagate through to all subsequent stages. 1988-1992 by Numerical Recipes A value di of any level is termed a "wavelet coefficient"of the original data 1-43108-.5 vector;the final values S1,S2 should strictly be called"mother-function coefficients," although the term"wavelet coefficients"is often used loosely for both d's and final (outside 膜 S's.Since the full procedure is a composition of orthogonal linear operations,the North Software. whole DWT is itself an orthogonal linear operator. To invert the DWT,one simply reverses the procedure,starting with the smallest level of the hierarchy and working(in equation 13.10.7)from right to left.The visit website inverse matrix(13.10.2)is of course used instead of the matrix(13.10.1). machine As already noted,the matrices(13.10.1)and(13.10.2)embody periodic("wrap- around")boundary conditions on the data vector.One normally accepts this as a minor inconvenience:the last few wavelet coefficients at each level of the hierarchy are affected by data from both ends of the data vector.By circularly shifting the matrix (13.10.1)N/2 columns to the left,one can symmetrize the wrap-around; but this does not eliminate it.It is in fact possible to eliminate the wrap-around
594 Chapter 13. Fourier and Spectral Applications Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). For higher p, up to 10, Daubechies [2] has tabulated the coefficients numerically. The number of coefficients increases by two each time p is 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, first to the full data vector of length N, then to the “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: 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 more stages of applying (13.10.1) (or any other wavelet coefficients) and permuting. 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 di of any level is termed a “wavelet coefficient” of the original data vector; the final values S1, S2 should strictly be called “mother-function coefficients,” 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. To invert the DWT, one simply reverses the procedure, starting with the smallest 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). As already noted, the matrices (13.10.1) and (13.10.2) embody periodic (“wraparound”) boundary conditions on the data vector. One normally accepts this as a minor inconvenience: the last few wavelet coefficients at each level of the hierarchy are affected by data from both ends of the data vector. By circularly shifting the matrix (13.10.1) N/2 columns to the left, one can symmetrize the wrap-around; but this does not eliminate it. It is in fact possible to eliminate the wrap-around
13.10 Wavelet Transforms 595 completely by altering the coefficients in the first and last N rows 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 if isign is 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. void wt1(float a[],unsigned long n,int isign, void(*wtstep)(f1oat☐，unsigned long,int)) One-dimensional discrete wavelet transform.This routine implements the pyramid algorithm, replacing a[1..n]by its wavelet transform (for isign=1),or performing the inverse operation (for isign=-1).Note that n MUST be an integer power of 2.The routine wtstep,whose actual name must be supplied in calling this routine,is the underlying wavelet filter.Examples of wtstep are daub4 and (preceded by pwtset)pwt. 3 unsigned long nn; RECIPES if (n 4)return; (Nort if(1s1gn>=0)[ Wavelet transform. for (nn=n;nn>=4;nn>>=1)(*wtstep)(a,nn,isign); server computer, Start at largest hierarchy,and work towards smallest. America one paper University Press. THE else Inverse wavelet transform. 9 ART for (nn=4;nn<=n;nn<<=1)(*wtstep)(a,nn,isign); Start at smallest hierarchy,and work towards largest. ictly proh Programs Here,as a specific instance of wtstep,is a routine for the DAUB4 wavelets: #include "nrutil.h" #define C00.4829629131445341 ectcustser 1881892 SCIENTIFIC COMPUTING(ISBN #def1neC10.8365163037378079 #def1neC20.2241438680420134 #define C3 -0.1294095225512604 10-621 void daub4(float a[],unsigned long n,int isign) Applies the Daubechies 4-coefficient wavelet filter to data vector a[1..n](for isign=1)or applies its transpose(for isign=-1).Used hierarchically by routines wt1 and wtn. Numerical Recipes 43106 float *wksp; unsigned long nh,nhi,i,j; (outside Software. if (n 4)return; wksp=vector(1,n); nh1=(nh=n>1)+1; 1f(1s1gm>=0){ Apply filter for(1=1,j=1;j<=n-3;j+=2,1++){ ksp[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] 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]; for (i=1,j=3;i<nh;i++){
13.10 Wavelet Transforms 595 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). completely by altering the coefficients in the first and last N rows 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 if isign is 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.) void wt1(float a[], unsigned long n, int isign, void (*wtstep)(float [], unsigned long, int)) One-dimensional discrete wavelet transform. This routine implements the pyramid algorithm, replacing a[1..n] by its wavelet transform (for isign=1), or performing the inverse operation (for isign=-1). Note that n MUST be an integer power of 2. The routine wtstep, whose actual name must be supplied in calling this routine, is the underlying wavelet filter. Examples of wtstep are daub4 and (preceded by pwtset) pwt. { unsigned long nn; if (n < 4) return; if (isign >= 0) { Wavelet transform. for (nn=n;nn>=4;nn>>=1) (*wtstep)(a,nn,isign); Start at largest hierarchy, and work towards smallest. } else { Inverse wavelet transform. for (nn=4;nn<=n;nn<<=1) (*wtstep)(a,nn,isign); Start at smallest hierarchy, and work towards largest. } } Here, as a specific instance of wtstep, is a routine for the DAUB4 wavelets: #include "nrutil.h" #define C0 0.4829629131445341 #define C1 0.8365163037378079 #define C2 0.2241438680420134 #define C3 -0.1294095225512604 void daub4(float a[], unsigned long n, int isign) Applies the Daubechies 4-coefficient wavelet filter to data vector a[1..n] (for isign=1) or applies its transpose (for isign=-1). Used hierarchically by routines wt1 and wtn. { float *wksp; unsigned long nh,nh1,i,j; if (n < 4) return; wksp=vector(1,n); nh1=(nh=n >> 1)+1; if (isign >= 0) { Apply filter. for (i=1,j=1;j<=n-3;j+=2,i++) { 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]; } 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]; for (i=1,j=3;i<nh;i++) {