c13-10 (779573), страница 5
Текст из файла (страница 5)
The two reproductionsfollowing are reconstructed with 23% (b) and 5.5% (c) of the 65536 waveletcoefficients. The latter image illustrates the kind of compromises made by thetruncated wavelet representation. High-contrast edges (the model’s right cheek andhair highlights, e.g.) are maintained at a relatively high resolution, while low-contrastareas (the model’s left eye and cheek, e.g.) are washed out into what amounts tolarge constant pixels.
Figure 13.10.4 (d) is the result of performing the identicalprocedure with Fourier, instead of wavelet, transforms: The figure is reconstructedfrom the 5.5% of 65536 real Fourier components having the largest magnitudes.One sees that, since sines and cosines are nonlocal, the resolution is uniformly pooracross the picture; also, the deletion of any components produces a mottled “ringing”everywhere. (Practical Fourier image compression schemes therefore break up animage into small blocks of pixels, 16 × 16, say, and do rather elaborate smoothingacross block boundaries when the image is reconstructed.)Fast Solution of Linear SystemsOne of the most interesting, and promising, wavelet applications is linearalgebra.
The basic idea [1] is to think of an integral operator (that is, a large matrix) asa digital image. Suppose that the operator compresses well under a two-dimensionalwavelet transform, i.e., that a large fraction of its wavelet coefficients are so smallas to be negligible. Then any linear system involving the operator becomes a sparseSample 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).Here, as before, wtstep is an individual wavelet step, either daub4 or pwt.604Chapter 13.Fourier and Spectral Applicationssystem in the wavelet basis. In other words, to solveA·x=b(13.10.16)we first wavelet-transform the operator A and the right-hand side b bye ≡ W · A · WT ,Aeb≡W·b(13.10.17)where W represents the one-dimensional wavelet transform, then solvee ·eAx=eb(13.10.18)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).Figure 13.10.4. (a) IEEE test image, 256 × 256 pixels with 8-bit grayscale. (b) The image is transformedinto the wavelet basis; 77% of its wavelet components are set to zero (those of smallest magnitude); itis then reconstructed from the remaining 23%. (c) Same as (b), but 94.5% of the wavelet componentsare deleted.
(d) Same as (c), but the Fourier transform is used instead of the wavelet transform. Waveletcoefficients are better than the Fourier coefficients at preserving relevant details.13.10 Wavelet Transforms605and finally transform to the answer by the inverse wavelet transformxx = WT · e(13.10.19)e(Note that the routine wtn does the complete transformation of A into A.)A typical integral operator that compresses well into wavelets has arbitrary (oreven nearly singular) elements near to its main diagonal, but becomes smooth awayfrom the diagonal. An example might be−1if i = j(13.10.20)Aij =|i − j|−1/2 otherwiseFigure 13.10.5 shows a graphical representation of the wavelet transform of thismatrix, where i and j range over 1 . .
. 256, using the DAUB12 wavelets. Elementslarger in magnitude than 10−3 times the maximum element are shown as blackpixels, while elements between 10−3 and 10−6 are shown in gray. White pixels are< 10−6 . The indices i and j each number from the lower left.In the figure, one sees the hierarchical decomposition into power-of-two sizedblocks. At the edges or corners of the various blocks, one sees edge effects causedby the wrap-around wavelet boundary conditions. Apart from edge effects, withineach block, the nonnegligible elements are concentrated along the block diagonals.This is a statement that, for this type of linear operator, a wavelet is coupled mainlySample 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).Figure 13.10.5. Wavelet transform of a 256 × 256 matrix, represented graphically. The original matrixhas a discontinuous cusp along its diagonal, decaying smoothly away on both sides of the diagonal. Inwavelet basis, the matrix becomes sparse: Components larger than 10−3 are shown as black, componentslarger than 10−6 as gray, and smaller-magnitude components are white.
The matrix indices i and jnumber from the lower left.606Chapter 13.Fourier and Spectral ApplicationsCITED 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 Mathematics, vol. 44, pp. 141–183. [1]Daubechies, I. 1988, Communications on Pure and Applied Mathematics, 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]13.11 Numerical Use of the Sampling TheoremIn §6.10 we implemented an approximating formula for Dawson’s integral due toRybicki. Now that we have become Fourier sophisticates, we can learn that the formuladerives from numerical application of the sampling theorem (§12.1), normally considered tobe 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 nranges over the integers and α is a constant that allows an arbitrary shift of the samplinggrid. We then writeg(t) =∞Xg(tn ) sincn=−∞π(t − tn ) + e(t)h(13.11.1)where sinc x ≡ sin x/x. The summation over the sampling points is called the samplingrepresentation of g(t), and e(t) is its error term. The sampling theorem asserts that thesampling representation is exact, that is, e(t) ≡ 0, if the Fourier transform of g(t),Z ∞G(ω) =g(t)eiωt dt(13.11.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. 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).to near neighbors in its own hierarchy (square blocks along the main diagonal) andnear neighbors in other hierarchies (rectangular blocks off the diagonal).The number of nonnegligible elements in a matrix like that in Figure 13.10.5scales only as N , the linear size of the matrix; as a rough rule of thumb it is about10N log10 (1/), where is the truncation level, e.g., 10−6 . For a 2000 by 2000matrix, 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] makethe interesting observations that (1) the product of two such matrices is itselfhierarchically band diagonal (truncating, of course, newly generated elements thatare smaller than the predetermined threshold ); and moreover that (2) the productcan be formed in order N operations.Fast matrix multiplication makes it possible to find the matrix inverse bySchultz’s (or Hotelling’s) method, see §2.5.Other schemes are also possible for fast solution of hierarchically band diagonalforms.
For example, one can use the conjugate gradient method, implemented in§2.7 as linbcg..
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
