|
3#
楼主 |
发表于 11.6.2006 22:39:48
|
只看该作者
#define N0 128<br />#include "stdio.h"<br />#include "stdlib.h"<br />#include "math.h"<br />#include "string.h"<br />void db4(double *h,double *g,double *hh,double *gg);<br />void wd(int N,double *h,double *g,double *c0,double *c,double *d);<br />void wr(int N,double *h,double *g,double *c, double *d,double *cd);<br />void main()<br />{<br />double fk[N0],c0[N0],c[N0],d[N0];<br />double h[8],g[8],hh[8],gg[8];<br />float fk0[N0];<br />FILE *fp;<br />int i,k,j,n,l,N;<br />fp=fopen("wdata.dat","rt");<br />fscanf(fp,"%d",&N);<br />for(k=0;kfclose(fp);<br />db4(h,g,hh,gg);<br />for(k=0;kc0[k]=fk0[k];<br />c[k]=0;<br />d[k]=0;<br />}<br />wd(N,hh,gg,c0,c,d);<br />wr(N,hh,gg,c,d,c0);<br />for(k=0;kreturn;<br />}<br />void wd(int N,double *h,double *g,double *c0,double *c,double *d)<br />/* wavelet decomposition */<br />{<br />int k,n,k2,l;<br />double ck,dk;<br />for(k=0;kck=0.0;<br />dk=0.0;<br />for(l=0;l<8;l++) {<br />n=k+l;<br />ck+=c0[n%N]*h[l];<br />dk+=c0[n%N]*g[l];<br />}<br />c[k]=ck;<br />d[k]=dk;<br />}<br />for(k=0;kk2=2*k;<br />c0[k]=c[k2];<br />c0[N/2+k]=d[k2];<br />}<br />return;<br />}<br />void wr(int N,double *h,double *g,double *c,double *d,double *c0)<br />/* wavelet reconstruction */<br />{<br />int k,n,l,k2;<br />double ck,cn,dn;<br />for(k=0;kk2=2*k;<br />c[k2]=c0[k];<br />c[k2+1]=0;<br />d[k2]=c0[N/2+k];<br />d[k2+1]=0;<br />}<br />for(k=0;kfor(k=0;kck=0.0;<br />for(l=0;l<8;l++) {<br />n=k-l;<br />cn=c[(N+n)%N];<br />dn=d[(N+n)%N];<br />ck+=cn*h[l]+dn*g[l];<br />}<br />c0[k]=ck;<br />}<br />return;<br />}<br />void db4(double *h,double *g,double *hh,double *gg)<br />/* Daubechies 4 wavelet */<br />{<br />int k,isgn;<br />h[7]=-0.0105974017850890;<br />h[6]= 0.0328830116668852;<br />h[5]= 0.0308413818355607;<br />h[4]=-0.1870348117190931;<br />h[3]=-0.0279837694168599;<br />h[2]= 0.6308807679398597;<br />h[1]= 0.7148465705529154;<br />h[0]= 0.2303778133088964;<br />isgn=1;<br />for(k=0;k<8;k++) {<br />gg[k]=isgn*h[7-k];<br />isgn=-isgn;<br />}<br />for(k=0;k<8;k++) {<br />g[k]=gg[7-k];<br />hh[k]=h[7-k];<br />}<br />return;<br />}<br />float fun(float x)<br />{<br />float pi=3.1415926;<br />float yx=30*exp(-x/40)*sin(2*pi*x/40);<br />return(yx);<br />}<br /><br />这个用的是 DB4小波,周期延拓,可以实现精确重构的。 |
|