小波与傅里叶分析基础_浅析傅里叶频谱分析&滤波&小波分析

原始信号由20Hz、50Hz和100Hz三种频率的正弦波组成，如上图所示。

在Figure图上画框截取待分析的数据，如下所示：

所得结果如下所示，

从上图中可以清晰的看出有三种频率，而且从时频图中可以看到三种频率发生在不 同时刻。

点击单边傅里叶频谱，放大频谱如下：

画框可以选定滤波范围，如带通滤波，选择如下框，可得结果：

MATLAB中具体代码如下：

1. function frequency_analysis

3. clc close all

5. Ts = 0.001;

7. Fs = 1/Ts;

9. f1 = 20;

11. f2 = 50;

13. f3 = 100;

15. dt = 0.2;

17. t1 = (0:Ts:dt-Ts) + 0; t2 = (0:Ts:dt-Ts) + dt;

19. t3 = (0:Ts:dt-Ts) + 2*dt;

21. y1 = sin(2*pi*f1*t1);

23. y2 = sin(2*pi*f2*t2);

25. y3 = sin(2*pi*f3*t3);

27. t = [t1 t2 t3];

29. y = [y1 y2 y3];

31. figure

33. plot(t,y)

35. xlim([t(1) t(end)])

37. ylim([min(y) max(y)])

39. xlabel('时间t')

41. ylabel('信号y(t)')

43. title('原始信号')

45. %

47. % 以下为标准化程序

49. % 上面的Ts和Fs要给定，后面会用到

51. set(gcf,'WindowButtonDownFcn',@BtndownFcn);

53. function BtndownFcn(h,evt)

55. temppt = get(gca,'CurrentPoint');

57. startpt.x = temppt(1,1);

59. startpt.y = temppt(1,2);

61. endpt = startpt;

63. height = 0.0001;

65. width = 0.0001;

67. rectangle('Position',[startpt.x, startpt.y, width, height],'Tag','rangerect'); set(h,'WindowButtonMotionFcn',@BtnmoveFcn);

69. set(h,'WindowButtonUpFcn',@BtnupFcn);

71. function BtnmoveFcn(h,evt)

73. temppt = get(gca,'CurrentPoint');

75. endpt.x = temppt(1,1);

77. endpt.y = temppt(1,2);

79. width = abs(endpt.x-startpt.x)+0.00001;

81. height = abs(endpt.y-startpt.y)+0.00001;

83. hrect = findobj('Tag','rangerect');

85. set(hrect,'Position',[min(startpt.x, endpt.x), min(startpt.y, endpt.y), width, height]);

87. end

89. function BtnupFcn(h,evt)

91. set(h,'WindowButtonMotionFcn','');

93. set(h,'WindowButtonUpFcn','');

95. hrect = findobj('Tag','rangerect');

97. delete(hrect);

99. BtnUp_Spectrum_Analysis(startpt.x, endpt.x)

101. end

103. % 频谱分析

105. function BtnUp_Spectrum_Analysis(starttime, endtime)

107. hline = findobj(gca,'type','line');

109. time = get(hline,'xdata');

111. laser = get(hline,'ydata');

113. % 得到截取的分析数据

115. tempx = find(time >= min(starttime, endtime));

117. startindex = tempx(1);

119. tempx = find(time <= max(starttime, endtime));

121. endindex = tempx(end);

123. yt = laser(startindex:endindex);

125. yt = ytmean(yt);

127. t = time(startindex:endindex);

129. % 傅里叶变换

131. [Yf, f] = Spectrum_Calc(yt,Fs);

133. % 小波变换

135. scale = 1:50;

137. cw2 = cwt(yt,scale,'morl');

139. % 作图

141. figure

143. subplot(231)

145. % 截取的频谱分析的数据

147. plot(t, yt)

149. xlim([t(1) t(end)])

151. ylim([min(yt),max(yt)]);

153. title('频谱分析数据')

155. xlabel('时间t')

157. ylabel('截取的数据y(t)')

159. h1 = subplot(234);

161. % 单边傅里叶变换分析频谱

163. plot(f,Yf,'-')

165. title('单边傅里叶频谱')

167. xlabel('频率Hz')

169. ylabel('|Y(f)|')

171. set(h1,'ButtonDownFcn',@BtnDown_Filter_fcn)

173. function BtnDown_Filter_fcn(h,evt)

175. fig_fft = figure;

177. plot(f,Yf,'-')

179. title('单边傅里叶频谱')

181. xlabel('频率Hz')

183. ylabel('|Y(f)|')

185. xlim([0 Fs/2])

187. ylim([0 max(Yf)])

189. % 构造右键菜单

191. filter_flag = 1;

193. hcmenu = uicontextmenu;

195. uimenu(hcmenu, 'Label', '低通滤波', 'Callback', @hcb1);

197. uimenu(hcmenu, 'Label', '高通滤波', 'Callback', @hcb2);

199. uimenu(hcmenu, 'Label', '带通滤波', 'Callback', @hcb3);

201. uimenu(hcmenu, 'Label', '带阻滤波', 'Callback', @hcb4);

203. set(gca,'UIContextMenu',hcmenu);

205. function hcb1(h,evt)

207. filter_flag = 1;

209. end

211. function hcb2(h,evt)

213. filter_flag = 2;

215. end

217. function hcb3(h,evt)

219. filter_flag = 3;

221. end

223. function hcb4(h,evt)

225. filter_flag = 4;

227. end

229. set(fig_fft,'WindowButtonDownFcn',@BtndownFcn);

231. function BtndownFcn(h,evt)

233. if

235. strcmp(get(h,'SelectionType'),'normal')

237. temppt = get(gca,'CurrentPoint');

239. startpt.x = temppt(1,1);

241. startpt.y = temppt(1,2);

243. endpt = startpt;

245. height = 0.0001;

247. width = 0.0001;

249. rectangle('Position',

251. [startpt.x, startpt.y, width, height],'Tag','rangerect_fft'); set(h,'WindowButtonMotionFcn',@BtnmoveFcn);

253. set(h,'WindowButtonUpFcn',@BtnupFcn);

255. end

257. function BtnmoveFcn(h,evt)

259. temppt = get(gca,'CurrentPoint');

261. endpt.x = temppt(1,1);

263. endpt.y = temppt(1,2);

265. width = abs(endpt.xstartpt.x)+0.00001;

267. height = abs(endpt.ystartpt.y)+0.00001;

269. hrect = findobj(h,'Tag','rangerect_fft');

271. set(hrect,'Position', [min(startpt.x, endpt.x), min(startpt.y, endpt.y), width, height]);

273. end

275. function BtnupFcn(h,evt)

277. set(h,'WindowButtonMotionFcn','');

279. set(h,'WindowButtonUpFcn','');

281. hrect = findobj(h,'Tag','rangerect_fft');

283. delete(hrect);

285. Filter_Analysis(startpt.x, endpt.x);

287. function Filter_Analysis(x1,x2)

289. lowfre = min(x1,x2);

291. highfre = max(x1,x2);

293. if

295. lowfre <= 0 || lowfre >= Fs/2

297. lowfre = 0.01;

299. end

301. if

303. highfre <= 0 || highfre >= Fs/2

305. highfre= Fs/2-0.01;

307. end

309. W1 = lowfre/(Fs/2);

311. W2 = highfre/(Fs/2);

313. filter_str = '(低通)';

315. switch filter_flag

317. case 1

319. [b,a] = butter(5,min(W1,W2)); % 低通，画的框的左边为截止频率

321. filter_str = '(低通)';

323. case 2

325. [b,a] = butter(5,max(W1,W2),'high'); % 高通，画的框的右边为截止频率

327. filter_str = '(高通)';

329. case 3

331. [b,a] = butter(5,[W1 W2]); % 带通

333. filter_str = '(带通)';

335. case 4

337. [b,a] = butter(5,[W1 W2],'stop'); % 带阻

339. filter_str = '(带阻)';

341. end

343. y = filter(b,a,yt);

345. figure

347. subplot(2,1,1)

349. plot(t,yt)

351. xlabel('时间t')

353. ylabel('信号y')

355. xlim([t(1) t(end)])

357. ylim([min(yt),max(yt)]);

359. title('原始信 号')

361. subplot(2,1,2)

363. plot(t,y)

365. xlabel('时间t')

367. ylabel('信号y')

369. xlim([t(1) t(end)])

371. ylim([min(y),max(y)]);

373. title(sprintf('滤波信号%s%.2fHz~%.2fHz',filter_str,lowfre,highfre))

375. end

377. end

379. end

381. end

383. subplot(1,3,[2,3]) % 频率轴化为频率

385. [X,Y] = meshgrid(t,5/(2*pi)./scale*Fs);

387. mesh(X,Y,abs(cw2))

389. view(0,90)

391. title('时频图')

393. xlabel('时间')

395. ylabel('频率')

397. xlim([t(1) t(end)])

399. set(gca,'ylim',[0,max(max(Y))])

401. set(gca,'YScale','log')

403. set(gca,'YTick', [1:9,10:10:90,100:100:900,1000,2000])

405. function [Yf, f] = Spectrum_Calc(yt,Fs)

407. L = length(yt);

409. NFFT = 2^nextpow2(L);

411. Yf = fft(yt,NFFT)/L;

413. Yf = 2*abs(Yf(1:NFFT/2+1));

415. f = Fs/2*linspace(0,1,NFFT/2+1);

417. end

419. end

421. end

423. end

425. 其中选框的代码可以标准化，可以用来在用户交互中用户选择范围，代码整理如下

427. set(gcf,'WindowButtonDownFcn',@BtndownFcn);

429. function BtndownFcn(h,evt)

431. temppt = get(gca,'CurrentPoint');

433. startpt.x = temppt(1,1);

435. startpt.y = temppt(1,2);

437. endpt = startpt;

439. height = 0.0001;

441. width = 0.0001;

443. rectangle('Position',[startpt.x, startpt.y, width, height],'Tag','rangerect'); set(h,'WindowButtonMotionFcn',@BtnmoveFcn);

445. set(h,'WindowButtonUpFcn',@BtnupFcn);

447. function BtnmoveFcn(h,evt)

449. temppt = get(gca,'CurrentPoint');

451. endpt.x = temppt(1,1);

453. endpt.y = temppt(1,2);

455. width = abs(endpt.x-startpt.x)+0.00001;

457. height = abs(endpt.y-startpt.y)+0.00001;

459. hrect = findobj('Tag','rangerect');

461. set(hrect,'Position',[min(startpt.x, endpt.x), min(startpt.y, endpt.y), width, height]);

463. end

465. function BtnupFcn(h,evt)

467. set(h,'WindowButtonMotionFcn','');

469. set(h,'WindowButtonUpFcn','');

471. hrect = findobj('Tag','rangerect');

473. delete(hrect);

475. % ProsessFcn(startpt,endpt);

477. % 选完框后对选择的 范围处理

479. end

481. end