c13-10 (Numerical Recipes in C), страница 4

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).1|H(ω)|2 + |H(ω + π)|2 = 1(13.10.10)2Likewise 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)60113.10 Wavelet Transforms0101.1.01.001.000110−510−610−71000100200200300300400500600400500600wavelet number70070080090010008009001000Figure 13.10.3.(a) Arbitrary test function, with cusp, sampled on a vector of length 1024.

(b)Absolute value of the 1024 wavelet coefficients produced by the discrete wavelet transform of (a). Notelog scale. The dotted curve plots the same amplitudes when sorted by decreasing size. One sees thatonly 130 out of 1024 coefficients are larger than 10−4 (or larger than about 10−5 times the largestcoefficient, whose value is ∼ 10).Truncated Wavelet ApproximationsMost of the usefulness of wavelets rests on the fact that wavelet transformscan usefully be severely truncated, that is, turned into sparse expansions. Thecase 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, anystandard mathematical operations that are especially simple in the wavelet basis.To see how truncation works, consider the simple example shown in Figure13.10.3. The upper panel shows an arbitrarily chosen test function, smooth exceptfor a square-root cusp, sampled onto a vector of length 210 . The bottom panel(solid curve) shows, on a log scale, the absolute value of the vector’s componentsafter 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.

Within each level, the wavelet coefficients are non-negligible only very near thelocation of the cusp, or very near the left and right boundaries of the hierarchicalrange (edge effects).The dotted curve in the lower panel of Figure 13.10.3 plots the same amplitudesas the solid curve, but sorted into decreasing order of size. One can read off, forexample, that the 130th largest wavelet coefficient has an amplitude less than 10−5of the largest coefficient, whose magnitude is ∼ 10 (power or square integral ratioless than 10−10 ).

Thus, the example function can be represented quite accuratelyby only 130, rather than 1024, coefficients — the remaining ones being set to zero.Note that this kind of truncation makes the vector sparse, but not shorter than 1024.It is very important that vectors in wavelet space be truncated according to theamplitude of the components, not their position in the vector. Keeping the first 256Sample 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).wavelet amplitude1.51.50602Chapter 13.Fourier and Spectral ApplicationsWavelet Transform in MultidimensionsA wavelet transform of a d-dimensional array is most easily obtained bytransforming the array sequentially on its first index (for all values of its other indices),then on its second, and so on.

Each transformation corresponds to multiplicationby an orthogonal matrix. By matrix associativity, the result is independent of theorder in which the indices were transformed. The situation is exactly like that formultidimensional FFTs. A routine for effecting the multidimensional DWT can thusbe modeled on a multidimensional FFT routine like fourn:#include "nrutil.h"void wtn(float a[], unsigned long nn[], int ndim, int isign,void (*wtstep)(float [], unsigned long, int))Replaces a by its ndim-dimensional discrete wavelet transform, if isign is input as 1. Herenn[1..ndim] is an integer array containing the lengths of each dimension (number of realvalues), which MUST all be powers of 2.

a is a real array of length equal to the product ofthese lengths, in which the data are stored as in a multidimensional real array. If isign is inputas −1, a is replaced by its inverse wavelet transform. The routine wtstep, whose actual namemust be supplied in calling this routine, is the underlying wavelet filter.

Examples of wtstepare daub4 and (preceded by pwtset) pwt.{unsigned long i1,i2,i3,k,n,nnew,nprev=1,nt,ntot=1;int idim;float *wksp;for (idim=1;idim<=ndim;idim++) ntot *= nn[idim];wksp=vector(1,ntot);for (idim=1;idim<=ndim;idim++) {Main loop over the dimensions.n=nn[idim];nnew=n*nprev;if (n > 4) {for (i2=0;i2<ntot;i2+=nnew) {for (i1=1;i1<=nprev;i1++) {for (i3=i1+i2,k=1;k<=n;k++,i3+=nprev) wksp[k]=a[i3];Copy the relevant row or column or etc.

into workspace.if (isign >= 0) {Do one-dimensional wavelet transform.for(nt=n;nt>=4;nt >>= 1)(*wtstep)(wksp,nt,isign);} else {Or inverse transform.for(nt=4;nt<=n;nt <<= 1)(*wtstep)(wksp,nt,isign);}for (i3=i1+i2,k=1;k<=n;k++,i3+=nprev) a[i3]=wksp[k];Copy back from workspace.}}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).components of the vector (all levels of the hierarchy except the last two) would givean extremely poor, and jagged, approximation to the function. When you compressa function with wavelets, you have to record both the values and the positions ofthe nonzero coefficients.Generally, compact (and therefore unsmooth) wavelets are better for loweraccuracy approximation and for functions with discontinuities (like edges), whilesmooth (and therefore noncompact) wavelets are better for achieving high numericalaccuracy. This makes compact wavelets a good choice for image compression, forexample, while it makes smooth wavelets best for fast solution of integral equations.13.10 Wavelet Transforms603}nprev=nnew;}free_vector(wksp,1,ntot);}Compression of ImagesAn immediate application of the multidimensional transform wtn is to imagecompression.

The overall procedure is to take the wavelet transform of a digitizedimage, and then to “allocate bits” among the wavelet coefficients in some highlynonuniform, optimized, manner. In general, large wavelet coefficients get quantizedaccurately, while small coefficients are quantized coarsely with only a bit or two— or else are truncated completely. If the resulting quantization levels are stillstatistically nonuniform, they may then be further compressed by a technique likeHuffman coding (§20.4).While a more detailed description of the “back end” of this process, namely thequantization and coding of the image, is beyond our scope, it is quite straightforwardto demonstrate the “front-end” wavelet encoding with a simple truncation: We keep(with full accuracy) all wavelet coefficients larger than some threshold, and we delete(set to zero) all smaller wavelet coefficients.

We can then adjust the threshold tovary the fraction of preserved coefficients.Figure 13.10.4 shows a sequence of images that differ in the number of waveletcoefficients that have been kept. The original picture (a), which is an official IEEEtest image, has 256 by 256 pixels with an 8-bit grayscale.