图像处理之基础---一维小波变换，可多次分解

1、题目：一维小波变换，可多次分解

2、原理：卷积核变为Daubechies正交小波基h[]和g[]的交替形式。增加了多次分解的功能。

3、代码：

[cpp] view plaincopy

1. #include <stdio.h>
2. #include <stdlib.h>
3. #include <math.h>
4. #define LENGTH 4096//信号长度
5. /*****************************************************************
6. *　一维卷积函数
7. *
8. *　说明: 循环卷积,卷积结果的长度与输入信号的长度相同
9. *
10. *　输入参数: data[],输入信号; h[],Daubechies小波基低通滤波器系数;
11. *　           g[],Daubechies小波基高通滤波器系数; cov[],卷积结果;
12. *            n,输入信号长度; m,卷积核长度.
13. *
14. *　李承宇, [email protected]
15. *
16. *  2010-08-22
17. *****************************************************************/
18. void Covlution(double data[], double h[], double g[], double cov[]
19. , int n, int m)
20. {
21. int i = 0;
22. int j = 0;
23. int k = 0;
24. //将cov[]清零
25. for(i = 0; i < n; i++)
26. {
27. cov[i] = 0;
28. }
29. //****************************************************
30. //奇数行用h[]进行卷积
31. //****************************************************
32. //前m/2+1行
33. i = 0;
34. for(j = 0; j < m/2; j+=2, i+=2)
35. {
36. for(k = m/2-j; k < m; k++ )
37. {
38. cov[i] += data[k-(m/2-j)] * h[k];//k针对core[k]
39. }
40. for(k = n-m/2+j; k < n; k++ )
41. {
42. cov[i] += data[k] * h[k-(n-m/2+j)];//k针对data[k]
43. }
44. }
45. //中间的n-m行
46. for( ; i <= (n-m)+m/2; i+=2)
47. {
48. for( j = 0; j < m; j++)
49. {
50. cov[i] += data[i-m/2+j] * h[j];
51. }
52. }
53. //最后m/2-1行
54. //  i = ( (n - m) + m/2 + 1 )/2*2;//**********
55. for(j = 1; j <= m/2; j+=2, i+=2)
56. {
57. for(k = 0; k < j; k++)
58. {
59. cov[i] += data[k] * h[m-j-k];//k针对data[k]
60. }
61. for(k = 0; k < m-j; k++)
62. {
63. cov[i] += h[k] * data[n-(m-j)+k];//k针对core[k]
64. }
65. }
66. //****************************************************
67. //偶数行用g[]进行卷积
68. //****************************************************
69. //前m/2+1行
70. i = 1;
71. for(j = 0; j < m/2; j+=2, i+=2)
72. {
73. for(k = m/2-j; k < m; k++ )
74. {
75. cov[i] += data[k-(m/2-j)] * g[k];//k针对core[k]
76. }
77. for(k = n-m/2+j; k < n; k++ )
78. {
79. cov[i] += data[k] * g[k-(n-m/2+j)];//k针对data[k]
80. }
81. }
82. //中间的n-m行
83. for( ; i <= (n-m)+m/2; i+=2)
84. {
85. for( j = 0; j < m; j++)
86. {
87. cov[i] += data[i-m/2+j] * g[j];
88. }
89. }
90. //最后m/2-1行
91. //  i = ( (n - m) + m/2 + 1 ) ;//*********
92. for(j = 1; j <= m/2; j+=2, i+=2)
93. {
94. for(k = 0; k < j; k++)
95. {
96. cov[i] += data[k] * g[m-j-k];//k针对data[k]
97. }
98. for(k = 0; k < m-j; k++)
99. {
100. cov[i] += g[k] * data[n-(m-j)+k];//k针对core[k]
101. }
102. }
103. }
104. /*****************************************************************
105. *   排序函数
106. *
107. *   将卷积后的结果进行排序，使尺度系数和小波系数分开
108. *****************************************************************/
109. void Sort(double data[], double sort[], int n)
110. {
111. for(int i = 0; i < n; i+=2)
112. {
113. sort[i/2] = data[i];
114. }
115. for(i = 1; i < n; i+=2)
116. {
117. sort[n/2+i/2] = data[i];
118. }
119. }
120. /*****************************************************************
121. *　一维小波变换函数
122. *
123. *　说明: 一维小波变换,可进行多次分解　
124. *
125. *　输入参数: input[],输入信号; output[],小波变换结果，包括尺度系数
126. *　和小波系数两部分; temp[],存放中间结果;h[],Daubechies小波基低通滤
127. *　波器系数;g[],Daubechies小波基高通滤波器系数;n,输入信号长度; m,
128. *　Daubechies小波基紧支集长度; nStep,小波变换分解次数
129. *
130. *　李承宇, [email protected]
131. *
132. *  2010-08-22
133. *****************************************************************/
134. void DWT1D(double input[], double output[], double temp[], double h[],
135. double g[], int n, int m, int nStep)
136. {
137. int i = 0;
138. for(i = 0; i < n; i++)
139. {
140. output[i] = input[i];
141. }
142. for(i = 0; i < nStep; i++)
143. {
144. Covlution(output, h, g, temp, n, m);
145. Sort(temp, output, n);
146. n = n/2;
147. }
148. }
149. void main()
150. {
151. double data[LENGTH];//输入信号
152. double temp[LENGTH];//中间结果
153. double data_output[LENGTH];//一维小波变换后的结果
154. int n = 0;//输入信号长度
155. int m = 6;//Daubechies正交小波基长度
156. int nStep = 6;//分解级数
157. int i = 0;
158. char s[32];//从txt文件中读取一行数据
159. static double h[] = {.332670552950, .806891509311, .459877502118,
160. -.135011020010, -.085441273882, .035226291882};
161. static double g[] = {.035226291882, .085441273882, -.135011020010,
162. -.459877502118, .806891509311, -.332670552950};
163. //读取输入信号
164. FILE *fp;
165. fp=fopen("data.txt","r");
166. if(fp==NULL) //如果读取失败
167. {
168. printf("错误！找不到要读取的文件/"data.txt/"/n");
169. exit(1);//中止程序
170. }
171. while( fgets(s, 32, fp) != NULL )//读取长度n要设置得长一点，要保证读到回车符，这样指针才会定位到下一行？回车符返回的是零值？是，非数字字符经过atoi变换都应该返回零值
172. {
173. //  fscanf(fp,"%d", &data[count]);//一定要有"&"啊！！！最后读了个回车符！适应能力不如atoi啊
174. data[n] = atof(s);
175. n++;
176. }
177. //一维小波变换
178. DWT1D(data, data_output, temp, h, g, n, m, nStep);
179. //一维小波变换后的结果写入txt文件
180. fp=fopen("test.txt","w");
181. //打印一维小波变换后的结果
182. for(i = 0; i < n/pow(2,nStep-1); i++)///pow(2,nStep-1)
183. {
184. printf("%f/n", data_output[i]);
185. fprintf(fp,"%f/n", data_output[i]);
186. }
187. //关闭文件
188. fclose(fp);
189. }

4、测试结果：

输入信号x(i)为：

取f1 = 5, f2 = 10, f0 = 320, n = 512。x(i)如图１所示：

图１　输入信号

各级分解的结果如图2~图7所示，左半部分为尺度系数，右半部分为小波系数：

图2　1级分解结果

图3　2级分解结果

图4　3级分解结果

图5　4级分解结果

图6　5级分解结果

图7　6级分解结果

图8是各级小波系数和第6级尺度系数的完整结果：

图8　第6级尺度系数和各级小波系数的完整结果