wtlib.c

su 的源代码库

			if (level==(StL1-1)) {

			  /* fprintf(stderr,"level()=%d\n",level); */

				for (i=0;i<ntot+1;i++) workcd[i]=0;
				for (iz=0;iz<sizez1/2;iz++)
			                for (jx=0;jx<sizex2;jx++) 
					        workcd[iz+jx*n1]=cd[iz+jx*n1];
				wtn(workcd,ToC,wfilt,wtsizes,justDC);
				for (i=0;i<ntot+1;i++) workf[i]=workcd[i];
				for (iz=0;iz<sizez1/2;iz++){
				        for (jx=0;jx<sizex2;jx++){
					        sum=0;
						for (i=0;i<leap;i++)
						        sum+=workf[iz*leap+i+jx*n1];
						f[iz+jx*n1]=sum/leap;
					}
				}    
			}
		}
	}

	if (StL2>StL1) {
	        for (level=StLevel;level<StL2;level++){
    
		  /* fprintf(stderr,"level=%d\n",level); */

			sizez1=n1>>StL1;
			sizex2=n2>>level;
			leap=2<<level;
			leap2=leap*leap;
              
			for (i=0;i<ntot+1;i++) workcd[i]=0;
			for (iz=0;iz<sizez1;iz++)
			        for (jx=sizex2/2;jx<sizex2;jx++)
				        workcd[iz+jx*n1]=cd[iz+jx*n1]; 
			wtn(workcd,ToC,wfilt,wtsizes,justDC);
			for (i=0;i<ntot+1;i++) workf[i]=workcd[i];
			for (iz=0;iz<sizez1;iz++){
			        for (jx=sizex2/2;jx<sizex2;jx++){
				        sum=0;
					for (i=0;i<leap;i++)
					        sum+=workf[iz+(jx*leap+i-n2)*n1];
					f[iz+jx*n1]=sum/leap;
				}
			}

			if (level==StL2-1) {

			        /* fprintf(stderr,"level()=%d\n",level); */

				for (i=0;i<ntot+1;i++) workcd[i]=0;
				for (iz=0;iz<sizez1;iz++)
				        for (jx=0;jx<sizex2/2;jx++) 
					        workcd[iz+jx*n1]=cd[iz+jx*n1];
				wtn(workcd,ToC,wfilt,wtsizes,justDC);
				for (i=0;i<ntot+1;i++) workf[i]=workcd[i];
				for (iz=0;iz<sizez1;iz++){
				        for (jx=0;jx<sizex2/2;jx++){
					        sum=0;
						for (i=0;i<leap;i++)
						        sum+=workf[iz+(jx*leap+i)*n1];
						f[iz+jx*n1]=sum/leap;
					}
				}
			}
		}
	}

	free1float(workf);
	free1float(workcd);
} 
		
void wto1d(float *cd,int npoints,enum ToCorD tocord,
            WtFilter *wfilt)
/*********************************************************************
wto1d - wavelet transform operator, 1D
*********************************************************************
Function Prototype:
void wto1d(float *cd,int npoints,enum ToCorD tocord,
            WtFilter *wfilt);
*********************************************************************
Input:
float *cd		pointer to wavelet coefficients, or f(x)
int npoints		size of the input signal
enum ToCorD tocord	=ToC inverse  or =ToD  forward transform
WtFilter *wfilt		pointer to wavelet operator (filters)

Returns:
float *cd		pointer to f(x) or to wavelet coefficients
*********************************************************************
Notes:
1D wavelet transform operator, used to apply multi-dimensional  
wavelet transform, in which data of each dimension is extracted 
and 1D wavelet transform (forward: tocord==ToD; backward: tocord 
==ToC) is applied to this 1-D data. Before using wto1d, we must 
call wto1dset to set up the wavelet filters.			   
							   
*********************************************************************
Reference:                                                    
W. Press et al., Numberical Recipes in C (second edition),     
Cambridge University Press, 1988				   
*********************************************************************
 Author:  Zhaobo Meng, 11/25/95,  Colorado School of Mines 
*********************************************************************/
{
        float ai,ai1,*temp;
	int   i,icd,j,jf,jr,k,n1,ni,nj,nhalf,nmod;

	temp=ealloc1float(npoints+1);
	nmod=wfilt->order*npoints;
	n1=npoints-1;
	nhalf=npoints>>1;
	for (j=1;j<=npoints;j++) temp[j]=0.0;
	if (tocord==ToD){
	        for (icd=1,i=1;i<=npoints;i+=2,icd++){
		        ni=i+nmod+wfilt->ioff;
			nj=i+nmod+wfilt->joff;
			for (k=1;k<=wfilt->order;k++){
			        jf=n1&(ni+k);
				jr=n1&(nj+k);
				temp[icd]+=wfilt->cc[k]*cd[jf+1];
				temp[icd+nhalf]+=wfilt->cr[k]*cd[jr+1];
			}
		}
	}
	else {
	        for (icd=1,i=1;i<=npoints;i+=2,icd++){
		        ai=cd[icd];
			ai1=cd[icd+nhalf];
			ni=i+nmod+wfilt->ioff;
			nj=i+nmod+wfilt->joff;
			for (k=1;k<=wfilt->order;k++){
		                jf=(n1&(ni+k))+1;
				jr=(n1&(nj+k))+1;
				temp[jf]+=wfilt->cc[k]*ai;
				temp[jr]+=wfilt->cr[k]*ai1;
			}
		}
	}
	for (j=1;j<=npoints;j++) cd[j]=temp[j]; 
	free1float(temp);
}

void wto1dset(WtFilter *wfilt,WtSizes *wtsizes)
/*********************************************************************
wto1dset - setup wavelet operators
*********************************************************************
Function Prototype:
void wto1dset(WtFilter *wfilt,WtSizes *wtsizes);
*********************************************************************
Input:
WtFilter *wfilt		pointer to wavelet filter
WtSizes *wtsizes	pointer to wavelet filter sizes

Return:
wfilt->order		order of Daubechies wavelet
wtsizes->sizes		sizes of the wavelet operator
*********************************************************************
Notes:
This subroutine initialize the wavelet filters for wavelet
transform operator wto1d. According to the order you choose
in your main (wtsizes->order=4 up to 20), wto1dset set up
the filter to be the Daubechies wavelet of this order/length. 
This subroutine will be need before calling wt1 or wtn. 

*********************************************************************
Reference: 
W. Press et al., Numberical Recipes in C (second edition), 
Cambridge University Press, 1988
*********************************************************************
Author:  Zhaobo Meng, 11/25/95,  Colorado School of Mines
*********************************************************************/
{
        int   k;
	int   order; /*number of coefficient=order of wavelet*/
	float sig=-1.0;
	static float c4[5]={0.0,0.4829629131445341,0.8365163037378079,
			    0.2241438680420314,-0.1294095225512604};
	static float c6[7]={0.0,0.332670552950,0.806891509311,
			    0.459877502118,-0.135011020010,-0.085441273882,
			    0.035226291882};
	static float c8[9]={0.0,0.230377813309,0.714846570553,0.630880767930,
			    -0.027983769417,-0.187034811719,0.030841381836,
			    0.032883011667,-0.010597401785};
	static float c10[11]={0.0,0.160102397974,0.603829269797,0.724308528438,
			      0.138428145901,-0.242294887066,-0.032244869585,
			      0.077571493840,-0.006241490213,-0.012580751999,
			      0.003335725285};
	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 c14[15]={0.0,0.077852054085,0.396539319482,0.729132090846,
			      0.469782287405,-0.143906003929,-0.224036184994,
			      0.071309219276,0.080612609151,-0.038029936935,
			      -0.016574541631,0.012550998556,0.000429577973,
			      -0.001801640704,0.000353713800};
	static float c16[17]={0.0,0.054415842243,0.312871590914,0.675930736297,
			      0.585354683654,-0.015829105056,-0.284015542962,
			      0.000472484574,0.128747426620,-0.017369301002,
			      -0.044088253931,0.013981027917,0.008746094047,
			      -0.004870352993,-0.000391740373,0.000675449406,
			      -0.000117476784};
	static float c18[19]={0.0,0.038077947364,0.243834674613,0.604823123690,
			      0.657288078051,0.133197385825,-0.293273783279,
			      -0.096840783223,0.148540749338,0.030725681479,
			      -0.067632829061,0.000250947115,0.022361662124,
			      -0.004723204758,-0.004281503682,0.001847646883,
			      0.000230385764,-0.000251963189,0.000039347320};
	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],c6r[7],c8r[9],c10r[11],c12r[13],c14r[15],
	        c16r[17],c18r[19],c20r[21];

	order=wtsizes->order;
	wfilt->order=order;
	switch(order) {
	case 4:
	        wfilt->cc=c4;
		wfilt->cr=c4r;
		break;
	case 6: 
	        wfilt->cc=c6;
		wfilt->cr=c6r;
		break;
	case 8:
	        wfilt->cc=c8;
		wfilt->cr=c8r;
		break;
	case 10:
	        wfilt->cc=c10;
		wfilt->cr=c10r;
		break;
	case 12:
	        wfilt->cc=c12;
		wfilt->cr=c12r;
		break;
	case 14:
	        wfilt->cc=c14;
		wfilt->cr=c14r;
		break;
	case 16:
	        wfilt->cc=c16;
		wfilt->cr=c16r;
		break;
	case 18:
	        wfilt->cc=c18;
		wfilt->cr=c18r;
		break;
	case 20:
	        wfilt->cc=c20;
		wfilt->cr=c20r;
		break;
	default:
	        fprintf(stderr,"unimplemented order in wto1dset");
		exit(0);
	}
	for (k=1;k<=order;k++){
	        wfilt->cr[wfilt->order+1-k]=sig*wfilt->cc[k];
		sig=-sig;
	}

	wfilt->ioff=wfilt->joff=-(order>>1);
  
	wt_cascade(wfilt,wtsizes);
}

void wt1(float *cd,enum ToCorD tocord,int npoints,
     WtFilter *wfilt,WtSizes *wtsizes,int idim)
/***********************************************************************
wt1 - 1D wavelet transform
***********************************************************************
Function Prototype:
void wt1(float *cd,enum ToCorD tocord,int npoints,
     WtFilter *wfilt,WtSizes *wtsizes,int idim);
***********************************************************************
Input:
float *cd		pointer to wavelet coeff. or input f(x)
enum ToCorD tocord	=ToC (invers) =ToD (forward) wavelet transform
int npoints		size of signal f(x) or cd
WtFilter *wfilt		pointer to wavelet operator (filters)
WtSizes *wtsizes	pointer to wavelet sizes
int idim		index for dimension

Returns:
float *cd		pointer to input f(x) or to wavelet coefficients
***********************************************************************
Notes:
1D wavelet transform. cd[1..npoints] as the discrete values of a 1D
function f(x) will be replaced by its wavelet transform if tocord
==ToD, or the wavelet coefficients cd[1..npoints] will be replaced
by the inverse wavelet transform, thus to obtain the reconstructed
f. The size of f npoints MUST be an integer power of 2.

***********************************************************************
Reference:
W. Press et al., Numberical Recipes in C (second edition),
Cambridge University Press, 1988
***********************************************************************
Author:  Zhaobo Meng, 11/25/95,  Colorado School of Mines
***********************************************************************/
{
        int   nthis;    /*length of this version*/
	int   lenlastc; /*length of the last c version*/
  
	lenlastc=npoints>>wtsizes->Mraleveln[idim];
	if (lenlastc<4) {
	        fprintf(stderr,"lenlastc=%d\n",lenlastc);
		fprintf(stderr,"StLevel too large or npoints too small.\n");
		exit(0);
	}

	if (tocord==ToD)
	        for (nthis=npoints;nthis>=2*lenlastc;nthis>>=1)wto1d(cd,nthis,ToD,wfilt);
	else
	        for (nthis=2*lenlastc;nthis<=npoints;nthis<<=1)wto1d(cd,nthis,ToC,wfilt);
}

void wtn(float *cd,enum ToCorD tocord,
     WtFilter *wfilt,WtSizes *wtsizes,int dconly)
/*************************************************************************
wtn - n-D wavelet transform 
**************************************************************************
Function Prototype:
void wtn(float *cd,enum ToCorD tocord,
     WtFilter *wfilt,WtSizes *wtsizes, int dconly) 
*************************************************************************
Input:
float *cd		pointer to n-D wavelet coeff, or f(x_1,..,x_n)
enum ToCorD tocord	=ToC (inverse) =ToD (forward) wavelet transform
WtFilter *wfilt		pointer to wavelet operator (filters)
WtSizes *wtsizes	pointer to wavelet sizes 
int dconly              keep and transform back from mra dc component 
                        only

Returns:
float *cd		pointer to f(x_1,...,x_n) or n-D wavelet coeffs.
*************************************************************************
Notes:
Given a n-D discrete values cd[...] of a n-D function f(x_1,...,x_n)
n-D wavelet transform will be applied if tocord==ToD; or given
a n-D wavelet coefficients of an n-D function, inverse wavelet
transform will be applied to reconstruct the n-D function f. Must
call wto1dset before calling wtn.
*************************************************************************
Reference: 
W. Press et al., Numberical Recipes in C (second edition), 
Cambridge University Press, 1988
*************************************************************************
Author:  Zhaobo Meng, 11/25/95,  Colorado School of Mines
Modified: Carlos E. Theodoro, 06/25/97, Colorado School of Mines
          - modified to transform back only dc component, if required 
*************************************************************************/
{
        long int   i1,i2,i3;
	int   k,n,nt,ij,ji,iz,jx;
	long int   nprev=1;	/* size of the previous version */
	long int   nnew;	/* size of the current version */
	long int   ntot=1;    /* size of the signal length */
	int  lenlastcinidim;	/* length of the last c version in idim-D */
	long int   idim;      /* dimension index */
	float *temp;		/* work space */
	float *temp1;		/* work space */

	for (idim=1;idim<=wtsizes->NDim;idim++) 
	        ntot*=wtsizes->NPointsn[idim];
  
	temp=ealloc1float(ntot+1);
	temp1=ealloc1float(ntot+1);

	for (idim=1;idim<=wtsizes->NDim;idim++) {
	        lenlastcinidim=(wtsizes->NPointsn[idim])>>(wtsizes->Mraleveln[idim]);
		n=wtsizes->NPointsn[idim];
		nnew=n*nprev;
		if (n>lenlastcinidim) {
		        for (i2=0;i2<ntot;i2+=nnew){
			        for (i1=1;i1<=nprev;i1++){
				        for (i3=i1+i2,k=1;k<=n;k++,i3+=nprev) 
					        temp[k]=cd[i3];
					if (tocord==ToC){
					        for (nt=lenlastcinidim;nt<=n;nt<<=1) 
						        wto1d(temp,nt,tocord,wfilt);
					}	
					else{
					        for (nt=n;nt>=lenlastcinidim;nt>>=1) 
						        wto1d(temp,nt,tocord,wfilt);
					}
					for (i3=i1+i2,k=1;k<=n;k++,i3+=nprev) 
					        cd[i3]=temp[k];
				}
			}
		}
		nprev=nnew;
	}
	if (dconly==1 && tocord==ToD) {
	        for (ij=1;ij<=ntot;ij++) temp1[ij]=cd[ij];
		for (ji=1;ji<=ntot;ji++) cd[ji]=0.;
    
		for (iz=1;iz<=((wtsizes->NPointsn[1])>>(wtsizes->Mraleveln[1]));iz++) {
		        for (jx=0;jx<((wtsizes->NPointsn[2])>>(wtsizes->Mraleveln[2]));jx++)
			        cd[iz+(wtsizes->NPointsn[1])*jx]=temp1[iz+(wtsizes->NPointsn[1])*jx];
		}
	}

	free1float(temp);
	free1float(temp1);
}