DB4小波实现重构

贴一个小波变换的C源代码。这个用的是 DB4小波，周期延拓。

#define N0 128
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "string.h"
void db4(double *h,double *g,double *hh,double *gg);
void wd(int N,double *h,double *g,double *c0,double *c,double *d);
void wr(int N,double *h,double *g,double *c, double *d,double *cd);
void main()
{
      double fk[N0],c0[N0],c[N0],d[N0];
      double h[8],g[8],hh[8],gg[8];
      float fk0[N0];
      FILE *fp;
      int i,k,j,n,l,N;
      fp=fopen("EEGdata.dat","rt");
      fscanf(fp,"%d",&N);
      for(k=0;k<N;k++)     fscanf(fp,"%f",&fk0[k]);
      fclose(fp);
      db4(h,g,hh,gg);
      for(k=0;k<N;k++)  {
           c0[k]=fk0[k];
           c[k]=0;
           d[k]=0;
       }
       wd(N,hh,gg,c0,c,d);
       wr(N,hh,gg,c,d,c0);
       for(k=0;k<N;k++)  printf("k=%d  c0=%f  c=%f\n",k,fk0[k],c0[k]);
       return;
}
void wd(int N,double *h,double *g,double *c0,double *c,double *d)
/*   wavelet decomposition  */
{
       int k,n,k2,l;
       double ck,dk;
       for(k=0;k<N;k++) {
            ck=0.0;
            dk=0.0;
            for(l=0;l<8;l++) {
               n=k+l;
              ck+=c0[n%N]*h[l];
              dk+=c0[n%N]*g[l];
             }
             c[k]=ck;
             d[k]=dk;
        }
       for(k=0;k<N/2;k++)  {
           k2=2*k;
           c0[k]=c[k2];
           c0[N/2+k]=d[k2];
       }
       return;
}
void wr(int N,double *h,double *g,double *c,double *d,double *c0)
/*   wavelet reconstruction  */
{
       int k,n,l,k2;
       double ck,cn,dn;
       for(k=0;k<N/2;k++)  {
           k2=2*k;
           c[k2]=c0[k];
           c[k2+1]=0;
           d[k2]=c0[N/2+k];
           d[k2+1]=0;
        }
        for(k=0;k<N;k++)  c0[k]=0.0;
        for(k=0;k<N;k++) {
            ck=0.0;
            for(l=0;l<8;l++) {
                n=k-l;
                cn=c[(N+n)%N];
                dn=d[(N+n)%N];
                ck+=cn*h[l]+dn*g[l];
            }
            c0[k]=ck;
         }
         return;
}
void db4(double *h,double *g,double *hh,double *gg)
/*  Daubechies 4 wavelet  */
{
        int k,isgn;
        h[7]=-0.0105974017850890;
        h[6]= 0.0328830116668852;
        h[5]= 0.0308413818355607;
        h[4]=-0.1870348117190931;
        h[3]=-0.0279837694168599;
        h[2]= 0.6308807679398597;
        h[1]= 0.7148465705529154;
        h[0]= 0.2303778133088964;
       isgn=1;
       for(k=0;k<8;k++) {
           gg[k]=isgn*h[7-k];
           isgn=-isgn;
        }
        for(k=0;k<8;k++) {
            g[k]=gg[7-k];
            hh[k]=h[7-k];
        }
        return;
}
float fun(float x)
{
       float pi=3.1415926;
       float yx=30*exp(-x/40)*sin(2*pi*x/40);
       return(yx);
}