贴一个小波变换算法的C源码,大家指点下!

wujunhoney

1. #define N0 128
   
2. #include "stdio.h"
   
3. #include "stdlib.h"
   
4. #include "math.h"
   
5. #include "string.h"
   
6. void db4(double *h,double *g,double *hh,double *gg);
   
7. void wd(int N,double *h,double *g,double *c0,double *c,double *d);
   
8. void wr(int N,double *h,double *g,double *c, double *d,double *cd);
   
9. void main()
   
10. {
    
11. double fk[N0],c0[N0],c[N0],d[N0];
    
12. double h[8],g[8],hh[8],gg[8];
    
13. float fk0[N0];
    
14. FILE *fp;
    
15. int i,k,j,n,l,N;
    
16. fp=fopen("wdata.dat","rt");
    
17. fscanf(fp,"%d",&N);
    
18. for(k=0;k<N;k++) fscanf(fp,"%f",&fk0[k]);
    
19. fclose(fp);
    
20. db4(h,g,hh,gg);
    
21. for(k=0;k<N;k++) {
    
22. c0[k]=fk0[k];
    
23. c[k]=0;
    
24. d[k]=0;
    
25. }
    
26. wd(N,hh,gg,c0,c,d);
    
27. wr(N,hh,gg,c,d,c0);
    
28. for(k=0;k<N;k++) printf("k=%d c0=%f c=%f\n",k,fk0[k],c0[k]);
    
29. return;
    
30. }
    
31. void wd(int N,double *h,double *g,double *c0,double *c,double *d)
    
32. /* wavelet decomposition */
    
33. {
    
34. int k,n,k2,l;
    
35. double ck,dk;
    
36. for(k=0;k<N;k++) {
    
37. ck=0.0;
    
38. dk=0.0;
    
39. for(l=0;l<8;l++) {
    
40.    n=k+l;
    
41.    ck+=c0[n%N]*h[l];
    
42.    dk+=c0[n%N]*g[l];
    
43. }
    
44. c[k]=ck;
    
45. d[k]=dk;
    
46. }
    
47. for(k=0;k<N/2;k++) {
    
48. k2=2*k;
    
49. c0[k]=c[k2];
    
50. c0[N/2+k]=d[k2];
    
51. }
    
52. return;
    
53. }
    
54. void wr(int N,double *h,double *g,double *c,double *d,double *c0)
    
55. /* wavelet reconstruction */
    
56. {
    
57. int k,n,l,k2;
    
58. double ck,cn,dn;
    
59. for(k=0;k<N/2;k++) {
    
60. k2=2*k;
    
61. c[k2]=c0[k];
    
62. c[k2+1]=0;
    
63. d[k2]=c0[N/2+k];
    
64. d[k2+1]=0;
    
65. }
    
66. for(k=0;k<N;k++) c0[k]=0.0;
    
67. for(k=0;k<N;k++) {
    
68. ck=0.0;
    
69. for(l=0;l<8;l++) {
    
70.    n=k-l;
    
71.    cn=c[(N+n)%N];
    
72.    dn=d[(N+n)%N];
    
73.   ck+=cn*h[l]+dn*g[l];
    
74. }
    
75. c0[k]=ck;
    
76. }
    
77. return;
    
78. }
    
79. void db4(double *h,double *g,double *hh,double *gg)
    
80. /* Daubechies 4 wavelet */
    
81. {
    
82. int k,isgn;
    
83. h[7]=-0.0105974017850890;
    
84. h[6]= 0.0328830116668852;
    
85. h[5]= 0.0308413818355607;
    
86. h[4]=-0.1870348117190931;
    
87. h[3]=-0.0279837694168599;
    
88. h[2]= 0.6308807679398597;
    
89. h[1]= 0.7148465705529154;
    
90. h[0]= 0.2303778133088964;
    
91. isgn=1;
    
92. for(k=0;k<8;k++) {
    
93. gg[k]=isgn*h[7-k];
    
94. isgn=-isgn;
    
95. }
    
96. for(k=0;k<8;k++) {
    
97. g[k]=gg[7-k];
    
98. hh[k]=h[7-k];
    
99. }
    
100. return;
     
101. }
     
102. float fun(float x)
     
103. {
     
104. float pi=3.1415926;
     
105. float yx=30*exp(-x/40)*sin(2*pi*x/40);
     
106. return(yx);
     
107. }

复制代码

//这个用的是 DB4小波，周期延拓，可以实现精确重构的

[ 本帖最后由 风花雪月 于 2008-6-17 15:01 编辑 ]

wr_9864

我回去验证一下
不知道您的计算周期是多少
我编的程序计算周期不满足要求
验证完毕后回来再顶你的贴

[ 本帖最后由 wr_9864 于 2008-3-22 23:44 编辑 ]

cike5

风花雪月

现在很多个人写的程序都是这样，缺乏必要的注释
过一段时间可能自己也忘记了

2365215

没看懂。。。。。。。。。。。。。。。。。。

graduate

一维的，而且没有采用提升算法
不知道我理解的对不？？

navy05

没有看懂呀，如果加上注释就好了