from http://blog.sina.com.cn/s/blog_6163bdeb0102e21a.html
最近看了看小波,小波理論半斤八兩,小波應用半斤八兩,先知道用吧,以后用的時(shí)候慢慢體會(huì )它的精髓吧??春芏嗟胤接眯〔?,基本就是簡(jiǎn)單的分解,然后重構各近似、細節信號,號稱(chēng)多分辨率分析,故障診斷中還常常配合Hilbert譜分析具體實(shí)現故障的提取?,F介紹一個(gè)簡(jiǎn)單例子吧!
3個(gè)正弦波(10、50、100Hz)混合,再摻和上噪聲,結果如下
小波分解,再重構各細節、近似信號如下
最后簡(jiǎn)單對一個(gè)細節信號d3做Hilbert譜分析
會(huì )發(fā)現,信號的倍頻通過(guò)小波分析會(huì )很明顯。
感覺(jué)小波的問(wèn)題是,不知道選哪個(gè)細節或近似,貌似就是靠目測。好吧,以后多用用,積累經(jīng)驗吧!
代碼如下
clc
clear all
close all
% 當前延拓模式是補零
oldmode = dwtmode('zpd');
ts = 0.001;
fs = 1/ts;
t=0:ts:1;
N = length(t);
x = sin(2*pi*10*t) + sin(2*pi*50*t) + sin(2*pi*100*t) +
0.1*randn(1, length(t));
figure
plot(t,x);
xlabel('時(shí)間 t/s');
ylabel('幅值 A');
% 一維小波分解
[c,l] = wavedec(x,5,'db3');
% 重構第1~5層逼近系數
a5 = wrcoef('a',c,l,'db3',5);
a4 = wrcoef('a',c,l,'db3',4);
a3 = wrcoef('a',c,l,'db3',3);
a2 = wrcoef('a',c,l,'db3',2);
a1 = wrcoef('a',c,l,'db3',1);
% 顯示逼近系數
figure
subplot(5,1,1);plot(a5);ylabel('a5');
subplot(5,1,2);plot(a4);ylabel('a4');
subplot(5,1,3);plot(a3);ylabel('a3');
subplot(5,1,4);plot(a2);ylabel('a2');
subplot(5,1,5);plot(a1);ylabel('a1');
xlabel('時(shí)間 t/s');
% 重構第1~5層細節系數
d5 = wrcoef('d',c,l,'db3',5);
d4 = wrcoef('d',c,l,'db3',4);
d3 = wrcoef('d',c,l,'db3',3);
d2 = wrcoef('d',c,l,'db3',2);
d1 = wrcoef('d',c,l,'db3',1);
% 顯示細節系數
figure
subplot(5,1,1);plot(d5);ylabel('d5');
subplot(5,1,2);plot(d4);ylabel('d4');
subplot(5,1,3);plot(d3);ylabel('d3');
subplot(5,1,4);plot(d2);ylabel('d2');
subplot(5,1,5);plot(d1);ylabel('d1');
xlabel('時(shí)間 t/s');
% 第1層細節信號的包絡(luò )譜
yh = hilbert(d3);
aabs =
abs(yh);
aabs = aabs - mean(aabs);
aangle =
unwrap(angle(yh));
af =
diff(aangle)/2/pi;
% NFFT = 2^nextpow2(N);
NFFT =
2^nextpow2(1024*4);
f = fs*linspace(0,1,NFFT);
YH = fft(yh,
NFFT)/N;
A = fft(aabs,
NFFT)/N;
figure
plot(f,abs(YH))
title('原始復信號的Hilbert譜')
xlabel('頻率f (Hz)')
ylabel('|YH(f)|')
figure
subplot(221)
plot(t, aabs')
title('包絡(luò )的絕對值')
legend('包絡(luò )分析結果', '真實(shí)包絡(luò )')
subplot(222)
plot(t, aangle)
title('調制信號的相位')
subplot(223)
plot(t(1:end-1), af*fs)
title('調制信號的瞬時(shí)頻率')
subplot(224)
plot(f,abs(A))
title('包絡(luò )的頻譜')
xlabel('頻率f (Hz)')
ylabel('|A(f)|')
% 恢復延拓模式
dwtmode(oldmode);
聯(lián)系客服