c13-10 (779573), страница 3
Текст из файла (страница 3)
Used hierarchically by routines wt1 and wtn.The actual filter is determined by a preceding (and required) call to pwtset, which initializesthe structure wfilt.{float ai,ai1,*wksp;unsigned long i,ii,j,jf,jr,k,n1,ni,nj,nh,nmod;if (n < 4) return;wksp=vector(1,n);nmod=wfilt.ncof*n;A positive constant equal to zero mod n.n1=n-1;Mask of all bits, since n a power of 2.nh=n >> 1;for (j=1;j<=n;j++) wksp[j]=0.0;if (isign >= 0) {Apply filter.for (ii=1,i=1;i<=n;i+=2,ii++) {ni=i+nmod+wfilt.ioff;Pointer to be incremented and wrapped-around.nj=i+nmod+wfilt.joff;for (k=1;k<=wfilt.ncof;k++) {jf=n1 & (ni+k);We use bitwise and to wrap-around the pointjr=n1 & (nj+k);ers.wksp[ii] += wfilt.cc[k]*a[jf+1];wksp[ii+nh] += wfilt.cr[k]*a[jr+1];}}} else {Apply transpose filter.for (ii=1,i=1;i<=n;i+=2,ii++) {ai=a[ii];ai1=a[ii+nh];ni=i+nmod+wfilt.ioff;See comments above.nj=i+nmod+wfilt.joff;for (k=1;k<=wfilt.ncof;k++) {jf=(n1 & (ni+k))+1;jr=(n1 & (nj+k))+1;wksp[jf] += wfilt.cc[k]*ai;wksp[jr] += wfilt.cr[k]*ai1;}}}for (j=1;j<=n;j++) a[j]=wksp[j];Copy the results back from workspace.free_vector(wksp,1,n);}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).}598Chapter 13.Fourier and Spectral Applications.1.050DAUB4 e5−.101002003004005006007008009001000.1.050−.05DAUB20 e22−.101002003004005006007008009001000Figure 13.10.1. Wavelet functions, that is, single basis functions from the wavelet families DAUB4and DAUB20.
A complete, orthonormal wavelet basis consists of scalings and translations of either oneof these functions. DAUB4 has an infinite number of cusps; DAUB20 would show similar behaviorin a higher derivative.What Do Wavelets Look Like?We are now in a position actually to see some wavelets. To do so, we simplyrun unit vectors through any of the above discrete wavelet transforms, with isignnegative so that the inverse transform is performed.
Figure 13.10.1 shows theDAUB4 wavelet that is the inverse DWT of a unit vector in the 5th component of avector of length 1024, and also the DAUB20 wavelet that is the inverse of the 22ndcomponent. (One needs to go to a later hierarchical level for DAUB20, to avoid awavelet with a wrapped-around tail.) Other unit vectors would give wavelets withthe same shapes, but different positions and scales.One sees that both DAUB4 and DAUB20 have wavelets that are continuous.DAUB20 wavelets also have higher continuous derivatives.
DAUB4 has the peculiarproperty that its derivative exists only almost everywhere. Examples of where itfails to exist are the points p/2n , where p and n are integers; at such points, DAUB4is left differentiable, but not right differentiable! This kind of discontinuity — atleast in some derivative — is a necessary feature of wavelets with compact support,like the Daubechies series. For every increase in the number of wavelet coefficientsby two, the Daubechies wavelets gain about half a derivative of continuity. (But notexactly half; the actual orders of regularity are irrational numbers!)Note that the fact that wavelets are not smooth does not prevent their havingexact representations for some smooth functions, as demanded by their approximationorder p. The continuity of a wavelet is not the same as the continuity of functionsSample 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).−.0559913.10 Wavelet Transforms.20DAUB4 e10 + e58010020030040050060070080090010004005006007008009001000.20−.2Lemarie e10 + e580100200300Figure 13.10.2. More wavelets, here generated from the sum of two unit vectors, e10 + e58 , whichare 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 defined bya filter most easily written in Fourier space (equation 13.10.14).that a set of wavelets can represent. For example, DAUB4 can represent (piecewise)linear functions of arbitrary slope: in the correct linear combinations, the cusps allcancel out.
Every increase of two in the number of coefficients allows one higherorder of polynomial to be exactly represented.Figure 13.10.2 shows the result of performing the inverse DWT on the inputvector e10 + e58 , again for the two different particular wavelets. Since 10 lies earlyin the hierarchical range of 9 − 16, that wavelet lies on the left side of the picture.Since 58 lies in a later (smaller-scale) hierarchy, it is a narrower wavelet; in the rangeof 33–64 it is towards the end, so it lies on the right side of the picture. Note thatsmaller-scale wavelets are taller, so as to have the same squared integral.Wavelet Filters in the Fourier DomainThe Fourier transform of a set of filter coefficients cj is given byH(ω) =Xcj eijω(13.10.8)jHere H is a function periodic in 2π, and it has the same meaning as before: It isthe wavelet filter, now written in the Fourier domain.
A very useful fact is that theorthogonality conditions for the c’s (e.g., equation 13.10.3 above) collapse to twosimple relations in the Fourier domain,1|H(0)|2 = 12(13.10.9)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).−.2600Chapter 13.Fourier and Spectral ApplicationsandH (m) (π) = 0m = 0, 1, .
. . , p − 1(13.10.11)It is thus relatively straightforward to invent wavelet sets in the Fourier domain.You simply invent a function H(ω) satisfying equations (13.10.9)–(13.10.11). Tofind the actual cj ’s applicable to a data (or s-component) vector of length N , andwith periodic wrap-around as in matrices (13.10.1) and (13.10.2), you invert equation(13.10.8) by the discrete Fourier transformcj =N−11 X2πk −2πijk/N)eH(NN(13.10.12)k=0The quadrature mirror filter G (reversed cj ’s with alternating signs), incidentally,has the Fourier representationG(ω) = e−iω H * (ω + π)(13.10.13)where asterisk denotes complex conjugation.In general the above procedure will not produce wavelet filters with compactsupport. In other words, all N of the cj ’s, j = 0, . .
. , N − 1 will in general benonzero (though they may be rapidly decreasing in magnitude). The Daubechieswavelets, or other wavelets with compact support, are specially chosen so that H(ω)is a trigonometric polynomial with only a small number of Fourier components,guaranteeing that there will be only a small number of nonzero cj ’s.On the other hand, there is sometimes no particular reason to demand compactsupport.
Giving it up in fact allows the ready construction of relatively smootherwavelets (higher values of p). Even without compact support, the convolutionsimplicit in the matrix (13.10.1) can be done efficiently by FFT methods.Lemarie’s wavelet (see [4]) has p = 4, does not have compact support, and isdefined by the choice of H(ω),1/2315 − 420u + 126u2 − 4u3H(ω) = 2(1 − u)4(13.10.14)315 − 420v + 126v2 − 4v3whereωv ≡ sin2 ω(13.10.15)u ≡ sin22It is beyond our scope to explain where equation (13.10.14) comes from. Aninformal description is that the quadrature mirror filter G(ω) deriving from equation(13.10.14) has the property that it gives identically zero when applied to any functionwhose odd-numbered samples are equal to the cubic spline interpolation of itseven-numbered samples.
Since this class of functions includes many very smoothmembers, it follows that H(ω) does a good job of truly selecting a function’s smoothinformation content. Sample Lemarie wavelets are shown in Figure 13.10.2.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.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
