小波包三层分解C语言源程序

此程序是用于信号处理分析，突出奇异值的前段处理，对信号进行小波包分解，用C语言实现的，这个程序相比MATLAB程序，没有进行边缘的处理，直接根据卷积原理进行编写而成的，如有对信号边缘处理的要求，可以自行改进。程序在DEV C++软件中运行成功，能够进行，程序中有注释部分，欢迎大家改进提高。

1. #include <stdio.h>
   
2. #include <stdlib.h>
   
3. #include <stdio.h>
   
4. #include <stdlib.h>
   
5. #define LENGTH 4096//信号长度
   
6. #define DB_LENGTH 8  //Daubechies小波基紧支集长度
   
7. /******************************************************************
   
8. *　一维卷积函数
   
9. *
   
10. *　说明: 循环卷积,卷积结果的长度与输入信号的长度相同
    
11. *
    
12. *　输入参数: data[],输入信号; core[],卷积核; cov[],卷积结果;
    
13. *            n,输入信号长度; m,卷积核长度.
    
14. *
    
15. ******************************************************************/
    
16. /*void Covlution(double data[], double core[], double cov[], int n, int m)
    
17. {
    
18.     int i = 0;
    
19.     int j = 0;
    
20.     int k = 0;
    
21. 
    
22.     //将cov[]清零
    
23.     for(i = 0; i < n; i++)
    
24.     {
    
25.         cov[i] = 0;
    
26.     }
    
27. 
    
28.     //前m/2+1行
    
29.     i = 0;
    
30.     for(j = 0; j < m/2; j++, i++)
    
31.     {
    
32.         for(k = m/2-j; k < m; k++ )
    
33.         {
    
34.             cov[i] += data[k-(m/2-j)] * core[k];//k针对core[k]
    
35.         }
    
36. 
    
37.         for(k = n-m/2+j; k < n; k++ )
    
38.         {
    
39.             cov[i] += data[k] * core[k-(n-m/2+j)];//k针对data[k]
    
40.         }
    
41.     }
    
42. 
    
43.     //中间的n-m行
    
44.     for( i = m/2; i <= (n-m)+m/2; i++)
    
45.     {
    
46.         for( j = 0; j < m; j++)
    
47.         {
    
48.             cov[i] += data[i-m/2+j] * core[j];
    
49.         }
    
50.     }
    
51. 
    
52.     //最后m/2-1行
    
53.     i = (n - m) + m/2 + 1;
    
54.     for(j = 1; j < m/2; j++, i++)
    
55.     {
    
56.         for(k = 0; k < j; k++)
    
57.         {
    
58.             cov[i] += data[k] * core[m-j-k];//k针对data[k]
    
59.         }
    
60. 
    
61.         for(k = 0; k < m-j; k++)
    
62.         {
    
63.             cov[i] += core[k] * data[n-(m-j)+k];//k针对core[k]
    
64.         }
    
65.     }
    
66. 
    
67. }
    
68. */
    
69. //定义一个线性卷积
    
70. void Covlution(double data[], double core[], double cov[], int n, int m)
    
71. {
    
72.          int i = 0;
    
73.     int j = 0;
    
74.     int t = 0;
    
75. 
    
76.     //将cov[]清零
    
77.     for(j = 0; j < n+m-1; j++)
    
78.     {
    
79.         cov[j] = 0;
    
80.     }
    
81. 
    
82.          for(j=0;j<m+n-1;j++)
    
83.          {
    
84.                  if(j<=m-1)       //前面m行
    
85.                  {
    
86.                         for(i=0,t=j;t>=0;i++,t--)
    
87.                          cov[j]+=data[i]*core[t];
    
88.                  
    
89.                  }
    
90.                  else if(j<=n-1)    //中间n-m行
    
91.                  {
    
92.                          for(i=j-m+1,t=m-1;t>=0;i++,t--)
    
93.                          cov[j]+=data[i]*core[t];
    
94.                  }
    
95.                  else     //后面m行
    
96.                           {
    
97.                                   for(i=j-m+1,t=m-1;i<n;i++,t--)
    
98.                                  cov[j]+=data[i]*core[t];        
    
99.                           }
    
100.          }        
     
101. 
     
102. }
     
103. 
     
104. 
     
105. /******************************************************************
     
106. *　一维小波变换函数
     
107. *
     
108. *　说明: 一维小波变换,只变换一次
     
109. *
     
110. *　输入参数: input[],输入信号; output[],小波变换结果，包括尺度系数和
     
111. *　小波系数两部分; temp[],存放中间结果;h[],Daubechies小波基低通滤波器系数;
     
112. *　g[],Daubechies小波基高通滤波器系数;n,输入信号长度; m,Daubechies小波基紧支集长度.
     
113. *
     
114. *　李承宇, lichengyu2345@126.com
     
115. *
     
116. *  2010-08-19  
     
117. ******************************************************************/
     
118. void DWT1D_1(double input[], double output0[], double output1[],double temp[], double h[],
     
119.            double g[], int n, int m)
     
120. {
     
121. //  double temp[LENGTH] = {0};//?????????????
     
122.      
     
123.     int i = 0;
     
124. /*
     
125.     //尺度系数和小波系数放在一起
     
126.     Covlution(input, h, temp, n, m);
     
127. 
     
128.     for(i = 0; i < n; i += 2)
     
129.     {
     
130.         output[i] = te mp[i];
     
131.     }
     
132. 
     
133.     Covlution(input, g, temp, n, m);
     
134. 
     
135.     for(i = 1; i < n; i += 2)
     
136.     {
     
137.         output[i] = temp[i];
     
138.     }
     
139. */
     
140.    //尺度系数和小波系数分开
     
141.    Covlution(input, h, temp, n, m);
     
142. 
     
143.     for(i = 0; i < n+m-1; i += 2)
     
144.     {
     
145.         output0[i/2] = temp[i];//尺度系数,进行了2抽值 ，即尺度空间，低频概貌部分
     
146.     }
     
147. 
     
148. Covlution(input, g, temp, n, m);
     
149. 
     
150.     for(i = 1; i < n+m-1; i +=2)
     
151.     {
     
152.     //   output[i] = temp[i];
     
153.            output1[(i-1)/2] = temp[i];//小波系数，已经进行了2抽取，即高频细节部分
     
154.     }
     
155. 
     
156. 
     
157. 
     
158. }
     
159. 
     
160. void DWT1D_2(double input[], double AA2[], double DA2[],double AD2[],double DD2[],double temp0[],double temp1[],double temp[], double h[],
     
161.            double g[], int n, int m)
     
162. {
     
163.                 DWT1D_1(input, temp0, temp1,temp, h, g, n, m);
     
164.                 int a1=(m+n)/2;
     
165.                 DWT1D_1(temp0, AA2, DA2,temp, h, g, a1, m);
     
166.                 int d1=(n+m-4)/2;
     
167.                 DWT1D_1(temp1, AD2, DD2,temp, h, g, d1, m);
     
168.                
     
169. }
     
170. 
     
171. void DWT1D_3(double input[], double AAA3[], double DAA3[],double ADA3[],double DDA3[],double AAD3[], double DAD3[],double ADD3[],double DDD3[],
     
172.                          double temp0[],double temp1[],double temp2[],double temp3[],double temp00[], double temp11[], double temp[], double h[], double g[], int n, int m)
     
173. {
     
174.                 DWT1D_2(input, temp0, temp1,temp2, temp3,temp00, temp11,temp, h, g, n, m);
     
175.                 int aa2=(m+n)/4;
     
176.                 DWT1D_1(temp0, AAA3, DAA3,temp, h, g, aa2, m);
     
177.                 int da2=(n+3*m-8)/4;
     
178.                 DWT1D_1(temp1, ADA3, DDA3,temp, h, g, da2, m);
     
179.                 int ad2=(n+3*m-4)/4;
     
180.                 DWT1D_1(temp2, AAD3, DAD3,temp, h, g, ad2, m);
     
181.                 int dd2=(n+3*m-12)/4;
     
182.                 DWT1D_1(temp3, ADD3, DDD3,temp, h, g, dd2, m);
     
183.                
     
184.                
     
185. }
     
186. 
     
187. 
     
188. void main()
     
189. {
     
190. 
     
191.     double data[LENGTH];//输入信号
     
192.     double temp0[(LENGTH+DB_LENGTH)/4]; //先定义了一个中间结果
     
193.     double temp1[(LENGTH+DB_LENGTH*3-8)/4];                        
     
194.     double temp2[(LENGTH+DB_LENGTH*3-4)/4];                     
     
195.     double temp3[(LENGTH+DB_LENGTH*3-12)/4];                       
     
196.     double temp00[(LENGTH+DB_LENGTH)/2];                          
     
197.     double temp11[(LENGTH+DB_LENGTH-4)/2];                                                     
     
198.     double temp[LENGTH+DB_LENGTH-1];
     
199.    
     
200.     double AAA3[(LENGTH+DB_LENGTH*5)/8];//一维小波变换后的结果
     
201.     double DAA3[(LENGTH+DB_LENGTH*5-16)/8];
     
202.    
     
203.         double ADA3[(LENGTH+DB_LENGTH*7-8)/8];   
     
204.     double DDA3[(LENGTH+DB_LENGTH*7-24)/8];
     
205.    
     
206.     double AAD3[(LENGTH+DB_LENGTH*7-4)/8];
     
207.         double DAD3[(LENGTH+DB_LENGTH*7-20)/8];
     
208.          
     
209.         double ADD3[(LENGTH+DB_LENGTH*7-12)/8];
     
210.     double DDD3[(LENGTH+DB_LENGTH*7-28)/8];
     
211.    
     
212.           int aaa3=(LENGTH+DB_LENGTH*5)/8;//一维小波变换后的结果数组的长度
     
213.     int daa3=(LENGTH+DB_LENGTH*5-16)/8;
     
214.    
     
215.         int ada3=(LENGTH+DB_LENGTH*7-8)/8;   
     
216.     int dda3=(LENGTH+DB_LENGTH*7-24)/8;
     
217.    
     
218.     int aad3=(LENGTH+DB_LENGTH*7-4)/8;
     
219.         int dad3=(LENGTH+DB_LENGTH*7-20)/8;
     
220.          
     
221.         int add3=(LENGTH+DB_LENGTH*7-12)/8;
     
222.     int ddd3=(LENGTH+DB_LENGTH*7-28)/8;
     
223.    
     
224.     int n = 0;//输入信号长度
     
225.     int m = 8;//Daubechies正交小波基长度
     
226.     int i = 0;
     
227.     char s[32];//从txt文件中读取一行数据
     
228. /*//DB3
     
229.     static double h[] = {.332670552950, .806891509311, .459877502118, -.135011020010,
     
230.                     -.085441273882, .035226291882};
     
231.     static double g[] = {.035226291882, .085441273882, -.135011020010, -.459877502118,
     
232.                     .806891509311, -.332670552950};
     
233. */
     
234. //DB4
     
235. 
     
236. static double h[] = {0.2303778133088964,   0.7148465705529154,         0.6308807679398587,-0.0279837694168599, -0.1870348117190931,     0.0308413818355607, 0.0328830116668852, -0.0105974017850690};//h[],Daubechies小波基低通滤波器系数;
     
237. static double g[] = {-0.0105974017850690, -0.0328830116668852, 0.0308413818355607, 0.1870348117190931,-0.0279837694168599, -0.6308807679398587, 0.7148465705529154, -0.2303778133088964 };//g[],Daubechies小波基高通滤波器系数
     
238.     //读取输入信号
     
239.     FILE *fp;
     
240.         FILE *fp0,*fp1,*fp2,*fp3,*fp4,*fp5,*fp6,*fp7;
     
241.     fp=fopen("heb1.txt","r");
     
242.     if(fp==NULL) //如果读取失败
     
243.     {
     
244.         printf("错误！找不到要读取的文件hea1.txt/n");
     
245.         exit(1);//中止程序
     
246.     }
     
247. 
     
248.     while( fgets(s, 32, fp) != NULL )//读取长度n要设置得长一点，要保证读到回车符，这样指针才会定位到下一行？回车符返回的是零值？是，非数字字符经过atoi变换都应该返回零值
     
249.     {
     
250.     //  fscanf(fp,"%d", &data[count]);//一定要有"&"啊！！！最后读了个回车符！适应能力不如atoi啊
     
251.         data[n] = atof(s);
     
252.         n++;
     
253.     }
     
254. 
     
255.     //一维小波变换
     
256.     //DWT1D(data, data_output, temp, h, g, n, m);
     
257.     DWT1D_3(data, AAA3, DAA3,ADA3,DDA3,AAD3, DAD3,ADD3,DDD3,
     
258.                          temp0, temp1, temp2, temp3, temp00, temp11, temp, h, g, n, m);
     
259. 
     
260.     //一维小波变换后的结果写入txt文件
     
261.     fp0=fopen("data_output_db4_aaa3.txt","w");
     
262.     fp1=fopen("data_output_db4_daa3.txt","w");
     
263.          fp2=fopen("data_output_db4_ada3.txt","w");
     
264.     fp3=fopen("data_output_db4_dda3.txt","w");
     
265.     fp4=fopen("data_output_db4_aad3.txt","w");
     
266.     fp5=fopen("data_output_db4_dad3.txt","w");
     
267.          fp6=fopen("data_output_db4_add3.txt","w");
     
268.     fp7=fopen("data_output_db4_ddd3.txt","w");
     
269.     //打印一维小波变换后的结果
     
270.     for(i = 0; i <aaa3; i++)
     
271.     {        //if(i==0)
     
272.                 printf("%s\n","AAA3");
     
273.         printf("%f\n", AAA3[i]);
     
274.         fprintf(fp0,"%f\n", AAA3[i]);
     
275. //                }
     
276. //        else
     
277. //        {
     
278. //        printf("%f\n", AAA3[i]);
     
279. //        fprintf(fp0,"%f\n", AAA3[i]);
     
280. //        }
     
281.   
     
282.     }
     
283. 
     
284.          for(i = 0; i <daa3; i++)
     
285.     {//if(i==0)
     
286.                 printf("%s\n","DAA3");
     
287.         printf("%f\n", DAA3[i]);
     
288.         fprintf(fp1,"%f\n", DAA3[i]);
     
289. //                }
     
290. //                else
     
291. //                {
     
292. //        printf("%f\n", DAA3[i]);
     
293. //        fprintf(fp1,"%f\n", DAA3[i]);
     
294. //                }
     
295.   
     
296.     }
     
297.    
     
298.         for(i = 0; i <ada3; i++)
     
299.     {//if(i==0)
     
300.                     printf("%s\n","ADA3");
     
301.         printf("%f\n", ADA3[i]);
     
302.         fprintf(fp2,"%f\n", ADA3[i]);
     
303. //            }
     
304. //                else
     
305. //                {
     
306. //        printf("%f\n", ADA3[i]);
     
307. //        fprintf(fp2,"%f\n", ADA3[i]);
     
308. //                }
     
309.   
     
310.     }
     
311.    
     
312.     for(i = 0; i <dda3; i++)
     
313.     {        //if(i==0)
     
314.                 printf("%s\n","DDA3");
     
315.         printf("%f\n", DDA3[i]);
     
316.         fprintf(fp3,"%f\n", DDA3[i]);
     
317. //                }
     
318. //                else
     
319. //                {
     
320. //        printf("%f\n", DDA3[i]);
     
321. //        fprintf(fp3,"%f\n", DDA3[i]);
     
322. //                }
     
323.   
     
324.     }
     
325.    
     
326.     for(i = 0; i <aad3; i++)
     
327.     {        //if(i==0)
     
328.                 printf("%s\n","AAD3");
     
329.         printf("%f\n", AAD3[i]);
     
330.         fprintf(fp4,"%f\n", AAD3[i]);        
     
331. //                }
     
332. //                else
     
333. //                {
     
334. //        printf("%f\n", AAD3[i]);
     
335. //        fprintf(fp4,"%f\n", AAD3[i]);
     
336. //                }
     
337.   
     
338.     }
     
339.         
     
340.         for(i = 0; i <dad3; i++)
     
341.     {        //if(i==0)
     
342.                 printf("%s\n","DAD3");
     
343.         printf("%f\n", DAD3[i]);
     
344.         fprintf(fp5,"%f\n", DAD3[i]);
     
345. //                }
     
346. //                else
     
347. //                {        
     
348. //        printf("%f\n", DAD3[i]);
     
349. //        fprintf(fp5,"%f\n", DAD3[i]);
     
350. //                }
     
351.   
     
352.     }
     
353.    
     
354.     for(i = 0; i <add3; i++)
     
355.     {        //if(i==0)
     
356.             printf("%s\n","ADD3");
     
357.         printf("%f\n", ADD3[i]);
     
358.         fprintf(fp6,"%f\n", ADD3[i]);
     
359. //            }
     
360. //            else
     
361. //            {
     
362. //        printf("%f\n", ADD3[i]);
     
363. //        fprintf(fp6,"%f\n", ADD3[i]);
     
364. //            }
     
365. //               
     
366.   
     
367.     }
     
368.    
     
369.     for(i = 0; i <ddd3; i++)
     
370.     {        //if(i==0)
     
371.                 printf("%s\n","DDD3");
     
372.         printf("%f\n", DDD3[i]);
     
373.         fprintf(fp7,"%f\n", DDD3[i]);
     
374. //                }
     
375. //                else
     
376. //                {
     
377. //        printf("%f\n", DDD3[i]);
     
378. //        fprintf(fp7,"%f\n", DDD3[i]);
     
379. //                }
     
380.   
     
381.     }
     
382.    
     
383.     //关闭文件
     
384.     fclose(fp);
     
385.     fclose(fp0);
     
386.     fclose(fp1);
     
387.     fclose(fp2);
     
388.     fclose(fp3);
     
389.     fclose(fp4);
     
390.     fclose(fp5);
     
391.     fclose(fp6);
     
392.     fclose(fp7);
     
393. }
     
394. 
     
395. 
     
396. /* run this program using the console pauser or add your own getch, system("pause") or input loop */
     
397. 
     
398. /* 尝试不对
     
399.    double A1[(LENGTH+DB_LENGTH)/2];
     
400.         double D1[(LENGTH+DB_LENGTH)/2];
     
401.         int a1=(LENGTH+DB_LENGTH)/2;
     
402.         int d1=(LENGTH+DB_LENGTH)/2;
     
403.         double AA2[(LENGTH+DB_LENGTH)/4];
     
404.         double DA2[(LENGTH+DB_LENGTH)/4];
     
405.         int a2=(LENGTH+DB_LENGTH)/4;
     
406.         int d2=(LENGTH+DB_LENGTH)/4;
     
407.         double AAA3[(LENGTH+DB_LENGTH)/8];
     
408.         double DAA3[(LENGTH+DB_LENGTH)/8];
     
409.         int a3=(LENGTH+DB_LENGTH)/8;
     
410.         int d3=(LENGTH+DB_LENGTH)/8;
     
411.         double AAAA4[(LENGTH+DB_LENGTH)/16];
     
412.         double DAAA4[(LENGTH+DB_LENGTH)/16];
     
413. 
     
414.     Covlution(input, h, temp, n, m);
     
415. 
     
416.     for(i = 0; i < n+m-2; i += 2)
     
417.     {
     
418.         A1[i/2] = temp[i];//一层分解的低频部分,进行了2抽值 ，即尺度空间，低频概貌部分
     
419.     }
     
420.      Covlution(input, g, temp, n, m);
     
421. 
     
422.     for(i =0 ; i < n+m-2; i +=2)
     
423.     {
     
424.            D1[i/2] = temp[i];//一层分解的高频部分，已经进行了2抽取，即高频细节部分
     
425.     }
     
426.    
     
427. 
     
428.       
     
429.     Covlution(A1, h, temp, a1, m);
     
430.     for(i = 0; i < a1+m-2; i += 2)
     
431.     {
     
432.         AA2[i/2] = temp[i];//一层分解的低频部分,进行了2抽值 ，即尺度空间，低频概貌部分
     
433.     }
     
434.    
     
435.     Covlution(A1, g, temp, a1, m);
     
436.     for(i = 0; i < a1+m-2; i += 2)
     
437.     {
     
438.         DA2[i/2] = temp[i];//一层分解的低频部分,进行了2抽值 ，即尺度空间，低频概貌部分
     
439.     }
     
440.    
     
441.    
     
442.    
     
443.      Covlution(AA2, h, temp, a2, m);
     
444.     for(i = 0; i < a2+m-2; i += 2)
     
445.     {
     
446.         AAA3[i/2] = temp[i];//一层分解的低频部分,进行了2抽值 ，即尺度空间，低频概貌部分
     
447.     }
     
448.    
     
449.     Covlution(AA2, g, temp, a2, m);
     
450.     for(i = 0; i < a2+m-2; i += 2)
     
451.     {
     
452.         DAA3[i/2] = temp[i];//一层分解的低频部分,进行了2抽值 ，即尺度空间，低频概貌部分
     
453.     }
     
454.    
     
455.    
     
456.    
     
457.    
     
458.    
     
459.      Covlution(AAA3, h, temp, a3, m);
     
460.     for(i = 0; i < a3+m-2; i += 2)
     
461.     {
     
462.         AAAA4[i/2] = temp[i];//一层分解的低频部分,进行了2抽值 ，即尺度空间，低频概貌部分
     
463.     }
     
464.    
     
465.     Covlution(AAA3, g, temp, a3, m);
     
466.     for(i = 0; i < a3+m-2; i += 2)
     
467.     {
     
468.         DAAA4[i/2] = temp[i];//一层分解的低频部分,进行了2抽值 ，即尺度空间，低频概貌部分
     
469.     }
     
470.    
     
471.      for(i = 0; i <4116; i++)
     
472.     {
     
473.             if(i<=259)
     
474.       output[i] = AAAA4[i];
     
475.       else if(i<=519)
     
476.       output[i] = DAAA4[i-260];
     
477.       else if(i<=1035)
     
478.       output[i] = DAA3[i-520];
     
479.       else if(i<=2064)
     
480.      output[i] = DA2[i-1036];
     
481.      else
     
482.      output[i] = DA2[i-2065];
     
483.    
     
484.     }
     
485. 
     
486. */        
     

复制代码

【必读】版权免责声明

1、本主题所有言论和内容纯属会员个人意见，与本论坛立场无关。2、本站对所发内容真实性、客观性、可用性不做任何保证也不负任何责任，网友之间仅出于学习目的进行交流。3、对提供的数字内容不拥有任何权利，其版权归原著者拥有。请勿将该数字内容进行商业交易、转载等行为，该内容只为学习所提供，使用后发生的一切问题与本站无关。 4、本网站不保证本站提供的下载资源的准确性、安全性和完整性;同时本网站也不承担用户因使用这些下载资源对自己和他人造成任何形式的损失或伤害。 5、本网站所有软件和资料均为网友推荐收集整理而来，仅供学习用途使用，请务必下载后两小时内删除，禁止商用。6、如有侵犯你版权的，请及时联系我们（电子邮箱1370723259@qq.com）指出，本站将立即改正。