#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);
}
- 1
- 2
前往页