c13-10 (779573), страница 2
Текст из файла (страница 2)
Notice that once d’s are generated, they simply propagate through to allsubsequent stages.A value di of any level is termed a “wavelet coefficient” of the original datavector; 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 finalS’s. Since the full procedure is a composition of orthogonal linear operations, thewhole DWT is itself an orthogonal linear operator.To invert the DWT, one simply reverses the procedure, starting with the smallestlevel of the hierarchy and working (in equation 13.10.7) from right to left.
Theinverse 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 aminor inconvenience: the last few wavelet coefficients at each level of the hierarchyare affected by data from both ends of the data vector. By circularly shifting thematrix (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-aroundSample 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).We have not yet defined the discrete wavelet transform (DWT), but we arealmost 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, andso on until only a trivial number of “smooth-.
. .-smooth” components (usually 2)remain. The procedure is sometimes called a pyramidal algorithm [4], for obviousreasons. The output of the DWT consists of these remaining components and allthe “detail” components that were accumulated along the way. A diagram shouldmake the procedure clear:13.10 Wavelet Transforms595void 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, whoseactual name must be supplied in calling this routine, is the underlying wavelet filter. Examplesof 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.1294095225512604void daub4(float a[], unsigned long n, int isign)Applies the Daubechies 4-coefficient wavelet filter to data vector a[1..n] (for isign=1) orapplies 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++) {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).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, beyondour scope here, is useful when, e.g., the data varies by many orders of magnitudefrom one end of the data vector to the other.Here is a routine, wt1, that performs the pyramidal algorithm (or its inverseif isign is negative) on some data vector a[1..n].
Successive applications ofthe wavelet filter, and accompanying permutations, are done by an assumed routinewtstep, which must be provided. (We give examples of several different wtsteproutines just below.)596Chapter 13.Fourier and Spectral Applicationswksp[j++]=C2*a[i]+C1*a[i+nh]+C0*a[i+1]+C3*a[i+nh1];wksp[j++] = C3*a[i]-C0*a[i+nh]+C1*a[i+1]-C2*a[i+nh1];}}for (i=1;i<=n;i++) a[i]=wksp[i];free_vector(wksp,1,n);}typedef struct {int ncof,ioff,joff;float *cc,*cr;} wavefilt;wavefilt wfilt;Defining declaration of a structure.void pwtset(int n)Initializing routine for pwt, here implementing the Daubechies wavelet filters with 4, 12, and20 coefficients, as selected by the input value n.
Further wavelet filters can be included in theobvious manner. This routine must be called (once) before the first use of pwt. (For the casen=4, the specific routine daub4 is considerably faster than pwt.){void nrerror(char error_text[]);int k;float sig = -1.0;static float c4[5]={0.0,0.4829629131445341,0.8365163037378079,0.2241438680420134,-0.1294095225512604};static float c12[13]={0.0,0.111540743350, 0.494623890398, 0.751133908021,0.315250351709,-0.226264693965,-0.129766867567,0.097501605587, 0.027522865530,-0.031582039318,0.000553842201, 0.004777257511,-0.001077301085};static float c20[21]={0.0,0.026670057901, 0.188176800078, 0.527201188932,0.688459039454, 0.281172343661,-0.249846424327,-0.195946274377, 0.127369340336, 0.093057364604,-0.071394147166,-0.029457536822, 0.033212674059,0.003606553567,-0.010733175483, 0.001395351747,0.001992405295,-0.000685856695,-0.000116466855,0.000093588670,-0.000013264203};static float c4r[5],c12r[13],c20r[21];wfilt.ncof=n;if (n == 4) {wfilt.cc=c4;wfilt.cr=c4r;}else if (n == 12) {wfilt.cc=c12;wfilt.cr=c12r;}else if (n == 20) {wfilt.cc=c20;wfilt.cr=c20r;}else nrerror("unimplemented value n in pwtset");for (k=1;k<=n;k++) {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).For larger sets of wavelet coefficients, the wrap-around of the last rows orcolumns is a programming inconvenience. An efficient implementation wouldhandle the wrap-arounds as special cases, outside of the main loop. Here, we willcontent ourselves with a more general scheme involving some extra arithmetic atrun time.
The following routine sets up any particular wavelet coefficients whosevalues you happen to know.13.10 Wavelet Transforms597wfilt.cr[wfilt.ncof+1-k]=sig*wfilt.cc[k];sig = -sig;}wfilt.ioff = wfilt.joff = -(n >> 1);These values center the “support” of the wavelets at each level. Alternatively, the “peaks”of the wavelets can be approximately centered by the choices ioff=-2 and joff=-n+2.Note that daub4 and pwtset with n=4 use different default centerings.Once pwtset has been called, the following routine can be used as a specificinstance of wtstep.#include "nrutil.h"typedef struct {int ncof,ioff,joff;float *cc,*cr;} wavefilt;extern wavefilt wfilt;Defined in pwtset.void pwt(float a[], unsigned long n, int isign)Partial wavelet transform: applies an arbitrary wavelet filter to data vector a[1..n] (for isign =1) or applies its transpose (for isign = −1).
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
