一维小波变换，可多次分解

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

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

3、代码：

view plaincopy to clipboardprint?

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. *　李承宇, [e-mail]
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. *　李承宇, [e-mail]
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. } 

#include <stdio.h>#include <stdlib.h>#include <math.h>#define LENGTH 4096//信号长度/******************************************************************　一维卷积函数**　说明: 循环卷积,卷积结果的长度与输入信号的长度相同**　输入参数: data[],输入信号; h[],Daubechies小波基低通滤波器系数;*　 g[],Daubechies小波基高通滤波器系数; cov[],卷积结果;* n,输入信号长度; m,卷积核长度.**　李承宇, [e-mail]** 2010-08-22 *****************************************************************/void Covlution(double data[], double h[], double g[], double cov[] , int n, int m){ int i = 0; int j = 0; int k = 0; //将cov[]清零 for(i = 0; i < n; i++) { cov[i] = 0; } //**************************************************** //奇数行用h[]进行卷积 //**************************************************** //前m/2+1行 i = 0; for(j = 0; j < m/2; j+=2, i+=2) { for(k = m/2-j; k < m; k++ ) { cov[i] += data[k-(m/2-j)] * h[k];//k针对core[k] } for(k = n-m/2+j; k < n; k++ ) { cov[i] += data[k] * h[k-(n-m/2+j)];//k针对data[k] } } //中间的n-m行 for( ; i <= (n-m)+m/2; i+=2) { for( j = 0; j < m; j++) { cov[i] += data[i-m/2+j] * h[j]; } } //最后m/2-1行// i = ( (n - m) + m/2 + 1 )/2*2;//********** for(j = 1; j <= m/2; j+=2, i+=2) { for(k = 0; k < j; k++) { cov[i] += data[k] * h[m-j-k];//k针对data[k] } for(k = 0; k < m-j; k++) { cov[i] += h[k] * data[n-(m-j)+k];//k针对core[k] } } //**************************************************** //偶数行用g[]进行卷积 //**************************************************** //前m/2+1行 i = 1; for(j = 0; j < m/2; j+=2, i+=2) { for(k = m/2-j; k < m; k++ ) { cov[i] += data[k-(m/2-j)] * g[k];//k针对core[k] } for(k = n-m/2+j; k < n; k++ ) { cov[i] += data[k] * g[k-(n-m/2+j)];//k针对data[k] } } //中间的n-m行 for( ; i <= (n-m)+m/2; i+=2) { for( j = 0; j < m; j++) { cov[i] += data[i-m/2+j] * g[j]; } } //最后m/2-1行// i = ( (n - m) + m/2 + 1 ) ;//********* for(j = 1; j <= m/2; j+=2, i+=2) { for(k = 0; k < j; k++) { cov[i] += data[k] * g[m-j-k];//k针对data[k] } for(k = 0; k < m-j; k++) { cov[i] += g[k] * data[n-(m-j)+k];//k针对core[k] } }}/****************************************************************** 排序函数** 将卷积后的结果进行排序，使尺度系数和小波系数分开*****************************************************************/void Sort(double data[], double sort[], int n){ for(int i = 0; i < n; i+=2) { sort[i/2] = data[i]; } for(i = 1; i < n; i+=2) { sort[n/2+i/2] = data[i]; }}/******************************************************************　一维小波变换函数**　说明: 一维小波变换,可进行多次分解　**　输入参数: input[],输入信号; output[],小波变换结果，包括尺度系数*　和小波系数两部分; temp[],存放中间结果;h[],Daubechies小波基低通滤*　波器系数;g[],Daubechies小波基高通滤波器系数;n,输入信号长度; m,*　Daubechies小波基紧支集长度; nStep,小波变换分解次数**　李承宇, [e-mail]** 2010-08-22 *****************************************************************/void DWT1D(double input[], double output[], double temp[], double h[], double g[], int n, int m, int nStep){ int i = 0; for(i = 0; i < n; i++) { output[i] = input[i]; } for(i = 0; i < nStep; i++) { Covlution(output, h, g, temp, n, m); Sort(temp, output, n); n = n/2; }}void main(){ double data[LENGTH];//输入信号 double temp[LENGTH];//中间结果 double data_output[LENGTH];//一维小波变换后的结果 int n = 0;//输入信号长度 int m = 6;//Daubechies正交小波基长度 int nStep = 6;//分解级数 int i = 0; char s[32];//从txt文件中读取一行数据 static double h[] = {.332670552950, .806891509311, .459877502118, -.135011020010, -.085441273882, .035226291882}; static double g[] = {.035226291882, .085441273882, -.135011020010, -.459877502118, .806891509311, -.332670552950}; //读取输入信号 FILE *fp; fp=fopen("data.txt","r"); if(fp==NULL) //如果读取失败 { printf("错误！找不到要读取的文件/"data.txt/"/n"); exit(1);//中止程序 } while( fgets(s, 32, fp) != NULL )//读取长度n要设置得长一点，要保证读到回车符，这样指针才会定位到下一行？回车符返回的是零值？是，非数字字符经过atoi变换都应该返回零值 { // fscanf(fp,"%d", &data[count]);//一定要有"&"啊！！！最后读了个回车符！适应能力不如atoi啊 data[n] = atof(s); n++; } //一维小波变换 DWT1D(data, data_output, temp, h, g, n, m, nStep); //一维小波变换后的结果写入txt文件 fp=fopen("test.txt","w"); //打印一维小波变换后的结果 for(i = 0; i < n/pow(2,nStep-1); i++)///pow(2,nStep-1) { printf("%f/n", data_output[i]); fprintf(fp,"%f/n", data_output[i]); } //关闭文件 fclose(fp);}

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

 

http://blog.csdn.net/lichengyu/article/details/5829785#