c13-10 (Numerical Recipes in C)

13.10 Wavelet Transforms591The splitting point b must be chosen large enough that the remaining integral over (b, ∞) issmall. Successive terms in its asymptotic expansion are found by integrating by parts. Theintegral over (a, b) can be done using dftint.

You keep as many terms in the asymptoticexpansion as you can easily compute. See [6] for some examples of this idea. Morepowerful methods, which work well for long-tailed functions but which do not use the FFT,are described in [7-9].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 TransformsLike the fast Fourier transform (FFT), the discrete wavelet transform (DWT) isa fast, linear operation that operates on a data vector whose length is an integer powerof two, transforming it into a numerically different vector of the same length.

Alsolike the FFT, the wavelet transform is invertible and in fact orthogonal — the inversetransform, 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, fromthe 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 thewavelet domain, the basis functions are somewhat more complicated and have thefanciful names “mother functions” and “wavelets.”Of course there are an infinity of possible bases for function space, almost all ofthem uninteresting! What makes the wavelet basis interesting is that, unlike sines andcosines, individual wavelet functions are quite localized in space; simultaneously,like sines and cosines, individual wavelet functions are quite localized in frequencyor (more precisely) characteristic scale. As we will see below, the particular kindof dual localization achieved by wavelets renders large classes of functions andoperators sparse, or sparse to some high accuracy, when transformed into the waveletdomain. Analogously with the Fourier domain, where a class of computations, likeconvolutions, become computationally fast, there is a large class of computationsSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).CITED REFERENCES AND FURTHER READING:Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag),p. 88.

[1]592Chapter 13.Fourier and Spectral ApplicationsDaubechies Wavelet Filter CoefficientsA particular set of wavelets is specified by a particular set of numbers, calledwavelet filter coefficients. Here, we will largely restrict ourselves to wavelet filtersin a class discovered by Daubechies [2]. This class includes members ranging fromhighly localized to highly smooth. The simplest (and most localized) member, oftencalled DAUB4, has only four coefficients, c0 , . .

. , c3 . For the moment we specializeto this case for ease of notation.Consider the following transformation matrix acting on a column vector ofdata to its right:c0 c3 . . .c2c1c1−c2...c2c1c0c3c3−c0c1−c2c2c1c3−c0...c0c3c3−c0c1−c2c2c1c0c3c3 −c0 c1−c2(13.10.1)Here blank entries signify zeroes. Note the structure of this matrix.

The first rowgenerates one component of the data convolved with the filter coefficients c0 . . . , 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 convolutionthat could be done by FFT methods. (Note how the last two rows wrap aroundlike convolutions with periodic boundary conditions.) Instead of convolving withc0 , . . .

, c3 , however, the even rows perform a different convolution, with coefficientsc3 , −c2 , c1 , −c0 . The action of the matrix, overall, is thus to perform two relatedconvolutions, 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 processingcontexts, H and G are called quadrature mirror filters [3].) In fact, the c’s are chosenso as to make G yield, insofar as possible, a zero response to a sufficiently smoothdata vector.

This is done by requiring the sequence c3 , −c2 , c1 , −c0 to have a certainnumber of vanishing moments. When this is the case for p moments (starting withthe zeroth), a set of wavelets is said to satisfy an “approximation condition of orderSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited.

To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).— those that can take advantage of sparsity — that become computationally fastin the wavelet domain [1].Unlike sines and cosines, which define a unique Fourier transform, there isnot one single unique set of wavelets; in fact, there are infinitely many possiblesets. Roughly, the different sets of wavelets make different trade-offs betweenhow compactly they are localized in space and how smooth they are.

(There arefurther fine distinctions.)59313.10 Wavelet Transformsc0 c1 c2 c3c3−c2c1−c0······c0c1c2c3c3−c2...c2c3c1−c0c0c1c2c3c3−c2c1−c0c0c1c1−c0 c (13.10.2)3−c2One sees immediately that matrix (13.10.2) is inverse to matrix (13.10.1) if andonly if these two equations hold,c20 + c21 + c22 + c23 = 1c 2 c0 + c3 c1 = 0(13.10.3)If additionally we require the approximation condition of order p = 2, then twoadditional relations are required,c 3 − c2 + c1 − c0 = 00c3 − 1c2 + 2c1 − 3c0 = 0(13.10.4)Equations (13.10.3) and (13.10.4) are 4 equations for the 4 unknowns c0 , .

. . , c3 ,first recognized and solved by Daubechies. The unique solution (up to a left-rightreversal) is√√√√c1 = (3 + 3)/4 2c0 = (1 + 3)/4 2(13.10.5)√√√√c3 = (1 − 3)/4 2c2 = (3 − 3)/4 2In fact, DAUB4 is only the most compact of a sequence of wavelet sets: If wehad six coefficients instead of four, there would be three orthogonality requirementsin equation (13.10.3) (with offsets of zero, two, and four), and we could requirethe vanishing of p = 3 moments in equation (13.10.4). In this case, DAUB6, thesolution coefficients can also be expressed in closed form,p√√√c0 = (1 + 10 + 5 p+ 2 10)/16 2√√√c2 = (10 − 2 10 +p2 5 + 2 10)/16 2√√√c4 = (5 + 10 − 3 5 + 2 10)/16 2p√√√c1 = (5 + 10 + 3 5p+ 2 10)/16 2√√√c3 = (10 − 2 10 −2 10)/16 2p2 5 +√√√c5 = (1 + 10 − 5 + 2 10)/16 2(13.10.6)Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).p.” This results in the output of H, decimated by half, accurately representing thedata’s “smooth” information. The output of G, also decimated, is referred to asthe data’s “detail” information [4].For such a characterization to be useful, it must be possible to reconstruct theoriginal data vector of length N from its N/2 smooth or s-components and its N/2detail or d-components. That is effected by requiring the matrix (13.10.1) to beorthogonal, so that its inverse is just the transposed matrix594Chapter 13.Fourier and Spectral ApplicationsFor higher p, up to 10, Daubechies [2] has tabulated the coefficients numerically.

Thenumber of coefficients increases by two each time p is increased by one.Discrete Wavelet Transformy1 y y23  y4  y5  y6  y7  y8  13.10.1 y9  −→ y10  y11  y12  y13  y14 y15y16 s1  s1  S1 S d s21   d2  s3   d3  s4   d4  permute s5  −→  d5  s6   d6  s7   d7 s8d8s s23   s4  s  13.10.1 5  −→ s6 s  7 s8  d1   d2  d3   d4  d5   d6 d7d8D S21   D2  permute S  −→ 3 D3 S  4 D4  d1   d2  d3   d4  d5   d6 d7d8Setc. S2  −→ 3 S4  D1   D2  D3   D4  d1   d2  d3   d4  d5   d6 d7d81S 1S2D 1 D2  D1   D2  D3   D4  d1   d2  d3   d4  d5   d6 d7d8(13.10.7)If the length of the data vector were a higher power of two, there would bemore 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.