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

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级尺度系数和各级小波系数的完整结果