[转帖]小波变换的C源代码

AaronSpark

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小波，周期延拓，可以实现精确重构的。

qiaoyinhu

感谢！

chidan

　　不错，感谢。

zyx168

　　非常感谢!!