par_wtlib.c

来自「该程序是用vc开发的对动态数组进行管理的DLL」· C语言 代码 · 共 828 行 · 第 1/2 页

			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=1.0f*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.4829629131445341f,0.8365163037378079f,
			    0.2241438680420314f,-0.1294095225512604f};
	static float c6[7]={0.0,0.332670552950f,0.806891509311f,
			    0.459877502118f,-0.135011020010f,-0.085441273882f,
			    0.035226291882f};
	static float c8[9]={0.0f,0.230377813309f,0.714846570553f,0.630880767930f,
			    -0.027983769417f,-0.187034811719f,0.030841381836f,
			    0.032883011667f,-0.010597401785f};
	static float c10[11]={0.0f,0.160102397974f,0.603829269797f,0.724308528438f,
			      0.138428145901f,-0.242294887066f,-0.032244869585f,
			      0.077571493840f,-0.006241490213f,-0.012580751999f,
			      0.003335725285f};
	static float c12[13]={0.0f,0.111540743350f,0.494623890398f,0.751133908021f,
			      0.315250351709f,-0.226264693965f,-0.129766867567f,
			      0.097501605587f,0.027522865530f,-0.031582039318f,
			      0.000553842201f,0.004777257511f,-0.001077301085f};
	static float c14[15]={0.0f,0.077852054085f,0.396539319482f,0.729132090846f,
			      0.469782287405f,-0.143906003929f,-0.224036184994f,
			      0.071309219276f,0.080612609151f,-0.038029936935f,
			      -0.016574541631f,0.012550998556f,0.000429577973f,
			      -0.001801640704f,0.000353713800f};
	static float c16[17]={0.0f,0.054415842243f,0.312871590914f,0.675930736297f,
			      0.585354683654f,-0.015829105056f,-0.284015542962f,
			      0.000472484574f,0.128747426620f,-0.017369301002f,
			      -0.044088253931f,0.013981027917f,0.008746094047f,
			      -0.004870352993f,-0.000391740373f,0.000675449406f,
			      -0.000117476784f};
	static float c18[19]={0.0f,0.038077947364f,0.243834674613f,0.604823123690f,
			      0.657288078051f,0.133197385825f,-0.293273783279f,
			      -0.096840783223f,0.148540749338f,0.030725681479f,
			      -0.067632829061f,0.000250947115f,0.022361662124f,
			      -0.004723204758f,-0.004281503682f,0.001847646883f,
			      0.000230385764f,-0.000251963189f,0.000039347320f};
	static float c20[21]={0.0f,0.026670057901f,0.188176800078f,0.527201188932f,
			      0.688459039454f,0.281172343661f,-0.249846424327f,
			      -0.195946274377f,0.127369340336f,0.093057364604f,
			      -0.071394147166f,-0.029457536822f,0.033212674059f,
			      0.003606553567f,-0.010733175483f,0.001395351747f,
			      0.001992405295f,-0.000685856695f,-0.000116466855f,
			      0.000093588670f,-0.000013264203f};
	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);
}