|
回复 12 # qqzhouqianziyi 的帖子
这样吧,我把以前做的两个信号叠加然后EMD分解后求HT谱和边际谱的程序附上吧。你可以参照。
clear;
fs=1000;
tspan=2;
t=1/fs:1/fs:tspan;
N=length(t);
x=sin(2*pi*20*t);
y=sin(2*pi*60*t+140);
z=x+y;
plot(t,z);
imf=emd(z);
cemd_visu(z,1:length(z),imf);
%%%imf1的Hilbert变换
xn1=hilbert(imf(1,:));
xr1=real(xn1);
xi1=imag(xn1);
A1=sqrt(xr1.^2+xi1.^2);
figure,subplot(2,1,1);plot(t,A1);
xlabel('时间(t)');ylabel('瞬时振幅');title('imf1')
%%%imf2的Hilbert变换
xn2=hilbert(imf(2,:));
xr2=real(xn2);
xi2=imag(xn2);
A2=sqrt(xr2.^2+xi2.^2);
subplot(2,1,2);plot(t,A2);
xlabel('时间(t)');ylabel('瞬时振幅');title('imf2')
%%%imf1的瞬时相位
P1=atan2(xi1,xr1);
figure,subplot(2,1,1);plot(t,P1);
xlabel('时间(t)'); ylabel('瞬时相位');title('imf1')
%%%imf2的瞬时相位
P2=atan2(xi2,xr2);
subplot(2,1,2);plot(t,P2);
xlabel('时间(t)'); ylabel('瞬时相位');title('imf2')
%%%%imf1瞬时频率
xh1=unwrap(P1);
fs=1000;
xhd1=fs*diff(xh1)/(2*pi);
figure,subplot(2,1,1);plot(t(1:1999),xhd1);
xlabel('时间(t)');ylabel('瞬时频率');title('imf1')
%%%%imf2瞬时频率
xh2=unwrap(P2);
fs=1000;
xhd2=fs*diff(xh2)/(2*pi);
subplot(2,1,2);plot(t(1:1999),xhd2);
xlabel('时间(t)');ylabel('瞬时频率');title('imf2')
>> [A,f,t]=hhspectrum(imf);
>> [E,t,cenf]=toimage(A,f);
>> disp_hhs(E,[],1000);
for k=1:size(E,1)
bjp(k)=sum(E(k,:))*1/tspan*1/fs;
end
H=size(E,1);
f=(0:H-1)/H*(fs/2);
figure ();
plot(f,bjp);
xlabel('频率 / Hz');
ylabel('幅值');title('边际谱')
出来的HT谱和边际谱图如下:
|
-
-
|