daubechi.c

Dwt离散db小波 daubechi.c daubtst.c

#ifndef MATH_H
#include <math.h>
#endif

#include <malloc.h>

/************************************************************************/
/*	Daubechies smoothing filter H					*/
/*		   diffentation filter G				*/
/************************************************************************/



void H(double *sig,double *Hsig,double *h,int length_sig,int length_h)
{int k,l;
 for(k=0;k<length_sig/2;k++)
     {Hsig[k]=0.0;
      for(l=0;l<length_h;l++)
          Hsig[k]+=h[l]*sig[(l+2*k)%length_sig];
      }
 }

void G(double *sig,double *Gsig,double *g,int length_sig,int length_h)
{int k,l;
 for(k=0;k<length_sig/2;k++)
     {Gsig[k]=0.0;
      for(l=0;l<length_h;l++)
          Gsig[k]+=g[l]*sig[(l+2*k)%length_sig];
      }
 }


/****************************************************************************/
/*	Daubechies adjoint filters H*, G*				    */
/****************************************************************************/

void Had(double *sig,double *Hsig,double *h,int length_Hsig,int length_h)
{int k,l;
 
for(k=0;k<2*length_Hsig;k+=2)
   {sig[k]=0.0;
    for(l=0;l<length_h;l+=2)
      sig[k]+=h[l]*Hsig[((k-l)/2+length_Hsig)%length_Hsig];
    }
for(k=1;k<2*length_Hsig;k+=2)
   {sig[k]=0.0;
    for(l=1;l<length_h;l+=2)
       sig[k]+=h[l]*Hsig[((k-l)/2+length_Hsig)%length_Hsig];
   }             
}   

void Gad(double *sig,double *Gsig,double *g,int length_Gsig,int length_g)
{int k,l;
 
for(k=0;k<2*length_Gsig;k+=2)
   {sig[k]=0.0;
    for(l=0;l<length_g;l+=2)
      sig[k]+=g[l]*Gsig[((k-l)/2+length_Gsig)%length_Gsig];
    }
for(k=1;k<2*length_Gsig;k+=2)
   {sig[k]=0.0;
    for(l=1;l<length_g;l+=2)
       sig[k]+=g[l]*Gsig[((k-l)/2+length_Gsig)%length_Gsig];
   }             
}   
 
/****************************************************************************/
/*	Initialization of Daubechies filters h,g			    */
/*									    */
/*	order		order of the wavelet to be initialized		    */
/*	length_filt	length of the corresponding filters		    */
/*	h,g		corresponding filter coefficients		    */
/*									    */
/*	return:	0	O.K.						    */
/*		-1	fiters of the desired order not implemented	    */	
/****************************************************************************/

int daub_init(int order,int *length_filt,double *h,double *g)
{int i,sig;
 switch(order)
  {case 1 : (*length_filt)=2;
            h[0]=1.0/sqrt(2.0); h[1]=h[0];
	    break;
  case 2: (*length_filt)=4;
          h[3]=(1.0-sqrt(3.0))/(4.0*sqrt(2.0));
	  h[2]=(3.0-sqrt(3.0))/(4.0*sqrt(2.0));
	  h[1]=(3.0+sqrt(3.0))/(4.0*sqrt(2.0));
	  h[0]=(1.0+sqrt(3.0))/(4.0*sqrt(2.0));
	  break;
 case 3: (*length_filt)=6;
         h[0]=0.3326705529500825;
	 h[1]=0.8068915093110924;
	 h[2]=0.4598775021184914;
	 h[3]=-0.1350110200102546;
	 h[4]=-0.0854412738820267;
	 h[5]=0.0352262918857095;
	 break;
 case 4: (*length_filt)=8;
         h[0]=0.2303778133088964;
	 h[1]=0.7148465705529154;
	 h[2]=0.6308807679398587;
	 h[3]=-0.0279837694168599;
	 h[4]=-0.1870348117190931;
	 h[5]=0.0308413818355607;
         h[6]=0.0328830116668852;
	 h[7]=-0.0105974017850690;
	 break;
 default: return(1);
  }
sig=-1;  
for(i=0;i<*length_filt;i++)  	 
 {sig*=-1;g[i]=sig*h[*length_filt-1-i];}
return(0);
} 
	 
	 
/****************************************************************************/
/*	Wavelt decomposition						    */
/*									    */
/*	sig		signal to be decomposed				    */
/*	length_sig	signal length					    */
/*	length_filt	filter_length (set by daub_init)		    */
/*	h		smoothing filter (set by daub_init)		    */
/*	g		differencing filter (set by daub_init)		    */
/*									    */
/*	output:		sig overwritten by the wavelet spectrum		    */
/*			decomp : 0 	OK				    */
/*				 1	not enough memory		    */
/****************************************************************************/

int decomp(double *sig,int length_sig,int length_filt,double *h,double *g)
{
 double *Hsig, *Gsig;
 int i,j;
 
 /*	allocation of buffer						    */         	 
          	          	  
 Hsig=(double *) malloc(length_sig/2*sizeof(double));
 Gsig=(double *)malloc(length_sig/2*sizeof(double));
 if(Hsig==NULL || Gsig==NULL) return(1);

/*	Wavelet decomposition						*/
 
 while(length_sig>=length_filt)
  {H(sig,Hsig,h,length_sig,length_filt);
   G(sig,Gsig,g,length_sig,length_filt);
   for(i=0;i<length_sig/2;i++)
     {sig[i]=Hsig[i];sig[length_sig/2+i]=Gsig[i];}
   length_sig/=2;
  }       

/*	release memory							*/  

free(Hsig);free(Gsig);

return(0);
} 


/*************************************************************************/
/*	Wavelet reconstruction						 */
/*									 */
/*	spec		spectrum used for reconstruction		 */
/*	length_spec	spectrum length			         	    */
/*	length_filt	filter_length (set by daub_init)		    */
/*	h		smoothing filter (set by daub_init)		    */
/*	g		differencing filter (set by daub_init)		    */
/*									    */
/*	output:		spec overwritten by the signal			    */
/*			reconst: 0 	OK				    */
/*				 1	not enough memory		    */
/****************************************************************************/

int reconst(double *spec,int length_spec,int length_filt,double *h,double *g) 
{int l,length;
 double *Hsig,*Gsig,*buff;
 
 Hsig=(double *)malloc(length_spec/2*sizeof(double));
 Gsig=(double *)malloc(length_spec/2*sizeof(double));
 buff=(double *)malloc(length_spec*sizeof(double));
 
 if(Hsig==NULL || Gsig==NULL || buff==NULL) return(1);
  
 switch(length_filt)
   {case 2 : length=1;
             break;
    case 4 : length=2;
	    break;
    case 6 : length=4;
            break;
    case 8 : length=4;
           break;	    
   }	    
	    	    	     
	     
	     
 while(length<=length_spec/2)
 {
  for(l=0;l<length;l++) {Hsig[l]=spec[l];Gsig[l]=spec[length+l];}
  Had(buff,Hsig,h,length,length_filt);
  Gad(spec,Gsig,g,length,length_filt);
  length*=2;
  for(l=0;l<length;l++) spec[l]+=buff[l];
 }
free(Hsig); free(Gsig); free(buff);
return(0);
}