596 Chapter 13. Fourier and Spectral Applications ksp[j++]=C2*a[i]+C1*a[i+nh]+C0*a[i+1]+C3*a[i+nh1]; ksp[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); 2 For larger sets of wavelet coefficients.the wrap-around of the last rows or columns is a programming inconvenience.An efficient implementation would handle the wrap-arounds as special cases,outside of the main loop.Here,we will content ourselves with a more general scheme involving some extra arithmetic at run time.The following routine sets up any particular wavelet coefficients whose nted for values you happen to know. 1800 typedef struct int ncof,ioff,joff; float *cc,*cri to any /Cambridge from NUMERICAL RECIPES IN 18881992 wavefilt; wavefilt wfilt; Defining declaration of a structure. (Nort void pwtset(int n) server computer, e University Press. THE Initializing routine for pwt,here implementing the Daubechies wavelet filters with 4,12,and America make one paper 20 coefficients,as selected by the input value n.Further wavelet filters can be included in the ART obvious manner.This routine must be called (once)before the first use of pwt.(For the case n=4,the specific routine daub4 is considerably faster than pwt.) 9 Programs void nrerror(char error_text[]); for their int k: float sig =-1.0; static f1oatc4[5]=[0.0,0.4829629131445341,0.8365163037378079, 0.2241438680420134,-0.1294095225512604): to dir stat1cf1oatc12[13]=[0.0,0.111540743350,0.494623890398,0.751133908021, 0.315250351709,-0.226264693965,-0.129766867567， 0.097501605587,0.027522865530,-0.031582039318 OF SCIENTIFIC COMPUTING(ISBN 0.000553842201,0.004777257511,-0.001077301085] stat1cf1oatc20[21]=10.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, Further reproduction,or 1988-1992 by Numerical Recipes 10-621 43108 0.001992405295,-0.000685856695,-0.000116466855, 0.000093588670,-0.000013264203}: static float c4r[5],c12r[13],c20r[21]; wfilt.ncof=n Software. if(n==4)[ wfilt.cc=c4; wfilt.cr=c4r; (outside North America) e1se1f(n=12)[ wfilt.cc=c12; wfilt.cr=c12r; e1se1f(n==20){ wfilt.cc=c20; wfilt.cr=c20r else nrerror("unimplemented value n in pwtset"); for(k=1;k<=n;k++){
596 Chapter 13. Fourier and Spectral Applications Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). wksp[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); } For larger sets of wavelet coefficients, the wrap-around of the last rows or columns is a programming inconvenience. An efficient implementation would handle the wrap-arounds as special cases, outside of the main loop. Here, we will content ourselves with a more general scheme involving some extra arithmetic at run time. The following routine sets up any particular wavelet coefficients whose values you happen to know. 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, and 20 coefficients, as selected by the input value n. Further wavelet filters can be included in the obvious manner. This routine must be called (once) before the first use of pwt. (For the case n=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++) {
13.10 Wavelet Transforms 597 wfilt.cr[wfilt.ncof+1-k]=sig*wfilt.cc[k]; s1g■-s1g: 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 putset with n-4 use different default centerings. Once pwtset has been called,the following routine can be used as a specific instance of wtstep. #include "nrutil.h" typedef struct int ncof,ioff,joff; granted for 18881992 float *cc,*cr; wavefilt; 100 extern wavefilt wfilt; Defined in pwtset. to any Cambridge from NUMERICAL RECIPES IN 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).Used hierarchically by routines wt1 and wtn. (Nort server The actual filter is determined by a preceding (and required)call to pwtset,which initializes the structure wfilt. America computer, THE float ai,ai1,*wksp; ART unsigned long i,ii,j,jf,jr,k,n1,ni,nj,nh,nmod; if (n 4)return; Programs 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; to dir 1f(1s1gn>=0)[ Apply filter. for(1i=1,1=1;i<n;i+=2,ii++)[ ni=i+nmod+wfilt.ioff; Pointer to be incremented and wrapped-around. OF SCIENTIFIC COMPUTING(ISBN nj=i+nmod+wfilt.joff; 1788-1982 for (k=1;k<=wfilt.ncof;k++){ jf=n1&(ni+k); We use bitwise and to wrap-around the point- jr=n1&(nj+k); ers. 10-621 wksp[ii]+wfilt.cc[k]*a[if+1]; wksp[ii+nh]+wfilt.cr[k]*a[jr+1]; .Further reproduction, Numerical Recipes 43108 else Apply transpose filter. for(1i=1,1=1;1<n;1+=2,11++)[ (outside 膜 aisa[ii]; ailsa[iitnh]; ni=i+nmod+wfilt.ioff; See comments above. nj=i+nmod+wfilt.joff; North Amer Software. for (k=1;k<=wfilt.ncof;k++){ jf=(n12(ni+k))+1; jr=(n1&(nj+k))+1: wksp[jf]+wfilt.cc[k]*ai; wksp[jr]+wfilt.cr[k]wail; for (j=1;j<=n;j++)a[j]=wksp[j]; Copy the results back from workspace. free_vector(wksp,1,n);
13.10 Wavelet Transforms 597 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). wfilt.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 specific instance 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). Used hierarchically by routines wt1 and wtn. The actual filter is determined by a preceding (and required) call to pwtset, which initializes the 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); }