[分享]个人收集的一些关于小波分析的matlab程序

simon21

1. clear;clc;
   
2. f=50;
   
3. T=0.001;
   
4. n=1:50;
   
5. y=sin(2*pi*f*n*T);
   
6. noise=[zeros(1,17),0.2*randn(1,15),zeros(1,18)];
   
7. y_noise=y+noise;
   
8. figure(1);
   
9. subplot(2,1,1);
   
10. plot(y);
    
11. title('signal');
    
12. subplot(2,1,2);
    
13. plot(y_noise);
    
14. title('noise & signal');
    
15. [yl2,yh2]=dwt(y,'db2');
    
16. [yl10,yh10]=dwt(y,'db10');
    
17. figure(2);
    
18. subplot(2,1,1);
    
19. plot(yl2);
    
20. title('db2 low frequency signal');
    
21. subplot(2,1,2);
    
22. plot(yh2);
    
23. title('db2 high frequency signal');
    
24. figure(3);
    
25. subplot(2,1,1);
    
26. plot(yl10);
    
27. title('db10 low frequency signal');
    
28. subplot(2,1,2);
    
29. plot(yh10);
    
30. title('db10 high frequency signal');
    
31. [yl2,yh2]=dwt(y_noise,'db2');
    
32. [yl10,yh10]=dwt(y_noise,'db10');
    
33. figure(4);
    
34. subplot(2,1,1);
    
35. plot(yl2);
    
36. title('db2 low frequency noise & signal');
    
37. subplot(2,1,2);
    
38. plot(yh2);
    
39. title('db2 high frequency noise & signal');
    
40. figure(5);
    
41. subplot(2,1,1);
    
42. plot(yl10);
    
43. title('db10 low frequency noise & signal');
    
44. subplot(2,1,2);
    
45. plot(yh10);
    
46. title('db10 high frequency noise & signal');

复制代码

[ 本帖最后由 simon21 于 2007-4-4 07:08 编辑 ]

simon21

1. clear;
   
2. clc; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
3. %% 1.原始信号 f=50; % 信号频率
   
4. fs=800; % 采样频率
   
5. N=128; % 采样点 % 信号赋值
   
6. n=1:N;
   
7. y=sin(2*pi*f*n/fs); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
8. %% 2.噪声
   
9. noise=0.4*rand(1,128); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   
10. %% 3.染噪信号
    
11. y_noise=y+noise; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
12. %% 4.硬消噪采用cycle_spinning技术 % 累加量
    
13. z5=zeros(1,N); % 平移变换频移法
    
14. for i=1:N;
    
15. z=circshift(y_noise.',i-1).'; % 源信号右平移
    
16. [z1,z2]=lwt(z,'db3'); % 小波正变换
    
17. z2=zeros(1,N/2); % 高频分量全部为零（主要噪声，硬消噪）
    
18. z3=ilwt(z1,z2,'db3'); % 小波反变换
    
19. z4=circshift(z3.',-(i-1)).'; % 变换后信号左平移
    
20. z5=z5+z4/N; % 平均
    
21. end;
    
22. 
    
23. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
24. %% 5.显示 error=norm(y-z5)/norm(y); % 相对误差 figure(1);
    
25. subplot(2,1,1)
    
26. plot(y,'r');
    
27. legend('源信号');
    
28. subplot(2,1,2);
    
29. plot(y_noise);
    
30. legend('染噪信号');
    
31. figure(2);
    
32. subplot(2,1,1)
    
33. plot(y,'r');
    
34. legend('源信号');
    
35. title(error);
    
36. subplot(2,1,2);
    
37. plot(z5);
    
38. legend('消噪后信号');

复制代码

[ 本帖最后由 simon21 于 2007-4-4 07:08 编辑 ]

cinda

　　2代小波示意程序 (三楼)中

　　……

　　f1(:,T+1)=f1(:,1); % 补列

　　f2(T+1,:)=f2(1,:); % 补行

　　……

　　补行是把分裂的偶数的第一行补到偶数的最后，但是补列是怎么进行的?比如图像是256×256，分裂后奇数为128×256，偶数也为128×256，进行f1(:,T+1)=f1(:,1); 补列之后，相当于把原来的第129列替换成第1列这样补行和补列后，奇数和偶数的分裂结果为128×256和129×256还请大侠不吝赐教，小生在此谢过了

　　[em24]

　　

[此贴子已经被作者于2006-5-24 16:09:37编辑过]

crosstree

多谢楼主，都是好东西~

cike5

其实我一直在努力,但是水平有限,还是先向各位学习,谢谢大家!

plutocn

很不错。。。是需要有人来总结和整理。

峰之巅

虽然不搞小波，但是为楼主的这种无私所敬佩！

zswseu

楼主真是太好了！！狂赞！！
博学多才！

wudan_1982

楼主真是大好人！能加楼主的QQ么？我的QQ573618825！谢谢楼主！:handshake

wudan_1982

我也是研究小波的，由于我身边没有人和我研究的一样，有好多问题想咨询楼主，在这里先谢谢楼主了。:@D

zhaorongman

matlab内容太多了,我是想把从实际中测得的信号,用matlab的分析功能进行分析,各位高手请指点,多谢谢!

tesionli

多谢楼主的无私奉献，在此感谢了/

ijkl56

顶一下

simon21

使用小波包变换分析信号的MATLAB程序

1. %t=0.001:0.001:1;
   
2. t=1:1000;
   
3. s1=sin(2*pi*50*t*0.001)+sin(2*pi*120*t*0.001)+rand(1,length(t));
   
4. 
   
5. for t=1:500;
   
6. s2(t)=sin(2*pi*50*t*0.001)+sin(2*pi*120*t*0.001)+rand(1,length(t));
   
7. end
   
8. 
   
9. for t=501:1000;
   
10. s2(t)=sin(2*pi*200*t*0.001)+sin(2*pi*120*t*0.001)+rand(1,length(t));
    
11. end
    
12. 
    
13. 
    
14. subplot(9,2,1)
    
15. plot(s1)
    
16. title('原始信号')
    
17. ylabel('S1')
    
18. 
    
19. subplot(9,2,2)
    
20. plot(s2)
    
21. title('故障信号')
    
22. ylabel('S2')
    
23. 
    
24. wpt=wpdec(s1,3,'db1','shannon');
    
25. 
    
26. %plot(wpt);
    
27. s130=wprcoef(wpt,[3,0]);
    
28. s131=wprcoef(wpt,[3,1]);
    
29. s132=wprcoef(wpt,[3,2]);
    
30. s133=wprcoef(wpt,[3,3]);
    
31. s134=wprcoef(wpt,[3,4]);
    
32. s135=wprcoef(wpt,[3,5]);
    
33. s136=wprcoef(wpt,[3,6]);
    
34. s137=wprcoef(wpt,[3,7]);
    
35. s10=norm(s130);
    
36. s11=norm(s131);
    
37. s12=norm(s132);
    
38. s13=norm(s133);
    
39. s14=norm(s134);
    
40. s15=norm(s135);
    
41. s16=norm(s136);
    
42. s17=norm(s137);
    
43. 
    
44. st10=std(s130);
    
45. st11=std(s131);
    
46. st12=std(s132);
    
47. st13=std(s133);
    
48. st14=std(s134);
    
49. st15=std(s135);
    
50. st16=std(s136);
    
51. st17=std(s137);
    
52. 
    
53. disp('正常信号的特征向量');
    
54. snorm1=[s10,s11,s12,s13,s14,s15,s16,s17]
    
55. std1=[st10,st11,st12,st13,st14,st15,st16,st17]
    
56. 
    
57. subplot(9,2,3);plot(s130);
    
58. ylabel('S130');
    
59. subplot(9,2,5);plot(s131);
    
60. ylabel('S131');
    
61. subplot(9,2,7);plot(s132);
    
62. ylabel('S132');
    
63. subplot(9,2,9);plot(s133);
    
64. ylabel('S133');
    
65. subplot(9,2,11);plot(s134);
    
66. ylabel('S134');
    
67. subplot(9,2,13);plot(s135);
    
68. ylabel('S135');
    
69. subplot(9,2,15);plot(s136);
    
70. ylabel('S136');
    
71. subplot(9,2,17);plot(s137);
    
72. ylabel('S137');
    
73. 
    
74. 
    
75. wpt=wpdec(s2,3,'db1','shannon');
    
76. 
    
77. %plot(wpt);
    
78. s230=wprcoef(wpt,[3,0]);
    
79. s231=wprcoef(wpt,[3,1]);
    
80. s232=wprcoef(wpt,[3,2]);
    
81. s233=wprcoef(wpt,[3,3]);
    
82. s234=wprcoef(wpt,[3,4]);
    
83. s235=wprcoef(wpt,[3,5]);
    
84. s236=wprcoef(wpt,[3,6]);
    
85. s237=wprcoef(wpt,[3,7]);
    
86. 
    
87. s20=norm(s230);
    
88. s21=norm(s231);
    
89. s22=norm(s232);
    
90. s23=norm(s233);
    
91. s24=norm(s234);
    
92. s25=norm(s235);
    
93. s26=norm(s236);
    
94. s27=norm(s237);
    
95. 
    
96. st20=std(s230);
    
97. st21=std(s231);
    
98. st22=std(s232);
    
99. st23=std(s233);
    
100. st24=std(s234);
     
101. st25=std(s235);
     
102. st26=std(s236);
     
103. st27=std(s237);
     
104. 
     
105. disp('故障信号的特征向量');
     
106. snorm2=[s20,s21,s22,s23,s24,s25,s26,s27]
     
107. std2=[st20,st21,st22,st23,st24,st25,st26,st27]
     
108. 
     
109. 
     
110. subplot(9,2,4);plot(s230);
     
111. ylabel('S230');
     
112. subplot(9,2,6);plot(s231);
     
113. ylabel('S231');
     
114. subplot(9,2,8);plot(s232);
     
115. ylabel('S232');
     
116. subplot(9,2,10);plot(s233);
     
117. ylabel('S233');
     
118. subplot(9,2,12);plot(s234);
     
119. ylabel('S234');
     
120. subplot(9,2,14);plot(s235);
     
121. ylabel('S235');
     
122. subplot(9,2,16);plot(s236);
     
123. ylabel('S236');
     
124. subplot(9,2,18);plot(s237);
     
125. ylabel('S237');
     
126. 
     
127. %fft
     
128. 
     
129. figure
     
130. y1=fft(s1,1024);
     
131. py1=y1.*conj(y1)/1024;
     
132. y2=fft(s2,1024);
     
133. py2=y2.*conj(y2)/1024;
     
134. 
     
135. y130=fft(s130,1024);
     
136. py130=y130.*conj(y130)/1024;
     
137. y131=fft(s131,1024);
     
138. py131=y131.*conj(y131)/1024;
     
139. y132=fft(s132,1024);
     
140. py132=y132.*conj(y132)/1024;
     
141. y133=fft(s133,1024);
     
142. py133=y133.*conj(y133)/1024;
     
143. y134=fft(s134,1024);
     
144. py134=y134.*conj(y134)/1024;
     
145. y135=fft(s135,1024);
     
146. py135=y135.*conj(y135)/1024;
     
147. y136=fft(s136,1024);
     
148. py136=y136.*conj(y136)/1024;
     
149. y137=fft(s137,1024);
     
150. py137=y137.*conj(y137)/1024;
     
151. 
     
152. y230=fft(s230,1024);
     
153. py230=y230.*conj(y230)/1024;
     
154. y231=fft(s231,1024);
     
155. py231=y231.*conj(y231)/1024;
     
156. y232=fft(s232,1024);
     
157. py232=y232.*conj(y232)/1024;
     
158. y233=fft(s233,1024);
     
159. py233=y233.*conj(y233)/1024;
     
160. y234=fft(s234,1024);
     
161. py234=y234.*conj(y234)/1024;
     
162. y235=fft(s235,1024);
     
163. py235=y235.*conj(y235)/1024;
     
164. y236=fft(s236,1024);
     
165. py236=y236.*conj(y236)/1024;
     
166. y237=fft(s237,1024);
     
167. py237=y237.*conj(y237)/1024;
     
168. 
     
169. f=1000*(0:511)/1024;
     
170. subplot(1,2,1);
     
171. plot(f,py1(1:512));
     
172. ylabel('P1');
     
173. title('原始信号的功率谱')
     
174. subplot(1,2,2);
     
175. plot(f,py2(1:512));
     
176. ylabel('P2');
     
177. title('故障信号的功率谱')
     
178. 
     
179. figure
     
180. 
     
181. subplot(4,2,1);
     
182. plot(f,py130(1:512));
     
183. ylabel('P130');
     
184. title('S130的功率谱')
     
185. subplot(4,2,2);
     
186. plot(f,py131(1:512));
     
187. ylabel('P131');
     
188. title('S131的功率谱')
     
189. subplot(4,2,3);
     
190. plot(f,py132(1:512));
     
191. ylabel('P132');
     
192. subplot(4,2,4);
     
193. plot(f,py133(1:512));
     
194. ylabel('P133');
     
195. subplot(4,2,5);
     
196. plot(f,py134(1:512));
     
197. ylabel('P134');
     
198. subplot(4,2,6);
     
199. plot(f,py135(1:512));
     
200. ylabel('P135');
     
201. subplot(4,2,7);
     
202. plot(f,py136(1:512));
     
203. ylabel('P136');
     
204. subplot(4,2,8);
     
205. plot(f,py137(1:512));
     
206. ylabel('P137');
     
207. figure
     
208. 
     
209. subplot(4,2,1);
     
210. plot(f,py230(1:512));
     
211. ylabel('P230');
     
212. title('S230的功率谱')
     
213. subplot(4,2,2);
     
214. plot(f,py231(1:512));
     
215. ylabel('P231');
     
216. title('S231的功率谱')
     
217. subplot(4,2,3);
     
218. plot(f,py232(1:512));
     
219. ylabel('P232');
     
220. subplot(4,2,4);
     
221. plot(f,py233(1:512));
     
222. ylabel('P233');
     
223. subplot(4,2,5);
     
224. plot(f,py234(1:512));
     
225. ylabel('P234');
     
226. subplot(4,2,6);
     
227. plot(f,py235(1:512));
     
228. ylabel('P235');
     
229. subplot(4,2,7);
     
230. plot(f,py236(1:512));
     
231. ylabel('P236');
     
232. subplot(4,2,8);
     
233. plot(f,py237(1:512));
     
234. ylabel('P237');
     
235. figure
     
236. %plottree(wpt)

复制代码

simon21

基于小波消噪的雷达回波检测，可以检测雷达回波的有无及其准确的位置

1. close all;clear;clc;
   
2. signalw=1200;%信号宽度
   
3. m=120000; %数据采样点数
   
4. pfa=0.001;%虚警概率
   
5. dectectw=signalw/4;%滑窗宽度
   
6. a=wgn(1,m,13,50,'dbm','real');%噪声
   
7. b=gate(a,m,dectectw,pfa)%门限
   
8. signal=real(fmlin(signalw,0,0.25));%信号
   
9. 
   
10. %加高斯白噪声信号
    
11. mix=zeros(1,m);
    
12. signalpos=40000;
    
13. for l=signalpos:signalpos+signalw-1
    
14.     mix(1,l)=signal(l-(signalpos-1),1);
    
15. end
    
16. mix=mix+a;%加高斯白噪声信号
    
17. 
    
18. %检测信号有无及位置
    
19. g=0;
    
20. for l=1:dectectw
    
21.     g=g+mix(l)*(mix(l))';
    
22.    
    
23. end
    
24.    g0=g;
    
25.    for k=(dectectw+1):m
    
26.     g=g+mix(k)*(mix(k))'-mix(k-dectectw)*(mix(k-dectectw))';
    
27.     if (g>b)
    
28.         over=1
    
29.         index=k
    
30.         break
    
31.     end
    
32. end%检测信号有无及位置
    
33. 
    
34. pnoise=mean(abs(a).^2);
    
35. psignal=mean(abs(signal).^2);
    
36. snr=10*log10(psignal/pnoise) %信噪比
    
37. 
    
38. [mix1,c,l]=wden(mix,'rigrsure','s','sln',3,'db2');
    
39. [a1,c,l]=wden(a,'rigrsure','s','sln',3,'db2');
    
40. % figure(2)
    
41. % subplot(2,1,1)
    
42. % plot(a)
    
43. % subplot(2,1,2)
    
44. % plot(a1)
    
45. [signal1,c,l]=wden(signal,'rigrsure','s','sln',3,'db2');
    
46. % figure(3)
    
47. % subplot(2,1,1)
    
48. % plot(signal)
    
49. % subplot(2,1,2)
    
50. % plot(signal1)
    
51. 
    
52. %检测信号有无及位置
    
53. dectectw1=signalw*1/4;
    
54. b1=gate(a1,m,dectectw1,pfa)
    
55. g=0;
    
56. for l=1:dectectw1
    
57.     g=g+mix1(l)*(mix1(l))';
    
58.    
    
59. end
    
60.    g0=g;
    
61.    for k=(dectectw1+1):m
    
62.     g=g+mix1(k)*(mix1(k))'-mix1(k-dectectw1)*(mix1(k-dectectw1))';
    
63.     if (g>b1)
    
64.         over=1
    
65.         index=k
    
66.         break
    
67.     end
    
68. end%检测信号有无及位置
    
69. 
    
70. pnoise1=(sum(abs(mix1).^2)-sum(abs(signal1).^2))/m;
    
71. psignal1=mean(abs(signal1).^2);
    
72. snr1=10*log10(psignal1/pnoise1) %信噪比
    
73. 
    
74. subplot(2,1,1)
    
75. plot(mix)
    
76. axis([0 60000 -4 4])
    
77. TITLE('消噪前信号')
    
78. str={['SNR=',num2str(snr),'dB']};
    
79. text(45000,3,str)
    
80. subplot(2,1,2)
    
81. plot(mix1)
    
82. axis([0 60000 -4 4])
    
83. TITLE('消噪后信号')
    
84. str={['SNR=',num2str(snr1),'dB']};
    
85. text(45000,2.5,str)

复制代码