#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);}