基为db4的小波变换（一维数组）

#define N0 128
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
void db4(double *h,double *g);
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);
float fun(float x);
void main()
{
      double fk[N0],c0[N0],c[N0],d[N0];
      double h[8],g[8];
      float fk0[N0];
      int k,N;

      N=N0;
      for(k=0;k<N;k++)    fk0[k]=fun(k);

      db4(h,g);
      for(k=0;k<N;k++)  {
          c0[k]=fk0[k];
          c[k]=0;
          d[k]=0;
      }
      wd(N,h,g,c0,c,d);
      wr(N,h,g,c,d,c0);
      for(k=0;k<N;k++)  printf("c=%f  c0=%f err=%f\n",c0[k],fk0[k],(c0[k]-fk0[k])*1.0e11);】

     system("pause");
      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,c0n;
      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++)  {
          c0[k]=c[2*k];
          c0[N/2+k]=d[2*k];
      }
      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++)  {
          c[2*k]=c0[k];
          c[2*k+1]=0;
          d[2*k]=c0[N/2+k];
          d[2*k+1]=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)
/*  Daubechies 4 wavelet  */
{
      int k,isgn;
/*   GIVED BY DAUBECHIES */
/*   g=(-1)**n*h(1-n)    */

      h[0]= 0.2303778133088964;
      h[1]= 0.7148465705529154;
      h[2]= 0.6308807679398597;
      h[3]=-0.0279837694168599;
      h[4]=-0.1870348117190931;
      h[5]= 0.0308413818355607;
      h[6]= 0.0328830116668852;
      h[7]=-0.0105974017850890;

      isgn=-1;
      for(k=0;k<8;k++) {
          g[k]=isgn*h[7-k];
          isgn=-isgn;
      }
      return;
}

float fun(float x)
{
      float pi=3.1415926;
      float yx=30*exp(-x/40)*sin(2*pi*x/40);
      return(yx);
}