clc;clear;

load('IR007_0.mat'); fs=12000; N=1024;

t=1/12000:1/12000:N/12000; x=X105_DE_time(1:1024,1)'; %EMD分解，展示各imf分量

imf=emd(x); x1=imf(1,:); x2=imf(2,:); x3=imf(3,:); x4=imf(4,:); x5=imf(5,:); x6=imf(6,:); x7=imf(7,:); x8=imf(8,:); x9=imf(9,:); figure(1);

subplot(4,1,1);plot(t,x1),title('imf1'); subplot(4,1,2);plot(t,x2),title('imf2'); subplot(4,1,3);plot(t,x3),title('imf3'); subplot(4,1,4);plot(t,x4),title('imf4'); xlabel('时间(time)t/s'),ylabel('幅值/mm'); figure(2);

subplot(5,1,1);plot(t,x5),title('imf5'); subplot(5,1,2);plot(t,x6),title('imf6'); subplot(5,1,3);plot(t,x7),title('imf7'); subplot(5,1,4);plot(t,x8),title('imf8'); subplot(5,1,5);plot(t,x9),title('res');

xlabel('时间(time)t/s'),ylabel('幅值/mm');

%计算imf的行列维数

[m,n]=size(imf);

for i=1:m-1

%-------------------------------------------------------------------- %计算各IMF的方差贡献率

%定义：方差为平方的均值减去均值的平方 %均值的平方

%imfp2=mean(c(i,:),2).^2 %平方的均值

%imf2p=mean(c(i,:).^2,2)

%各个IMF的方差

mse(i)=mean(imf(i,:).^2,2)-mean(imf(i,:),2).^2; end;

mmse=sum(mse);%总的方差 for i=1:m-1

%方差百分比，也就是方差贡献率

mseb(i)=mse(i)/mmse*100; end;

figure(3);

bar([1:m-1],mseb);%可以改进，最好能做出柱状图 title('各imf分量的方差贡献率'); xlabel('方差'),ylabel('百分比');

%原始信号的Hilbert变换 hx=hilbert(x);

xr=real(hx);xi=imag(hx); %计算瞬时振幅

A=sqrt(xr.^2+xi.^2);

figure(4);plot(t,A);title('瞬时振幅') %计算瞬时相位

P=angle(hx);

figure(5);plot(t,P);title('瞬时相位') %计算瞬时频率

dt=diff(t); dx=diff(P); sp=dx./dt;

figure(6);plot(t(1:N-1),sp);title('瞬时频率')

%imf1的Hilbert变换 xn1=hilbert(imf(1,:)); xr1=real(xn1); xi1=imag(xn1); %imf1的瞬时振幅

A1=sqrt(xr1.^2+xi1.^2); figure(7);

subplot(2,1,1);plot(t,A1);

xlabel('时间（t）');ylabel('瞬时振幅');title('imf1') %imf2的Hilbert变换 xn2=hilbert(imf(2,:)); xr2=real(xn2); xi2=imag(xn2); %imf2的瞬时振幅 A2=sqrt(xr2.^2+xi2.^2);

subplot(2,1,2);plot(t,A2);

xlabel('时间（t）');ylabel('瞬时振幅');title('imf2') %imf1的瞬时相位

P1=atan2(xi1,xr1); figure(8);

subplot(2,1,1);plot(t,P1);

xlabel('时间（t）'); ylabel('瞬时相位');title('imf1') %imf2的瞬时相位

P2=atan2(xi2,xr2);

subplot(2,1,2);plot(t,P2);www.hunanwang.net xlabel('时间（t）'); ylabel('瞬时相位');title('imf2') %imf1瞬时频率

xh1=unwrap(P1); fs=1000;

xhd1=fs*diff(xh1)/(2*pi); figure(9);

subplot(2,1,1);plot(t(1:1023),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:1023),xhd2); www.visa158.com xlabel('时间（t）');ylabel('瞬时频率');title('imf2')

%计算HHT的时频谱

[A, fa, tt] = hhspectrum(imf); if size(imf,1) > 1

[A,fa,tt] = hhspectrum(imf(1:end-1, :)); else

[A,fa,tt] = hhspectrum(imf); end

[E,tt1]=toimage(A,fa,tt,length(tt));

disp_hhs(E,[],fs) %二维图形显示HHT视频谱，E是求得的HHT谱

for i=1:m-1 faa=fa(i,:);

[FA,TT1]=meshgrid(faa,tt1);%三维图显示HHT时频图 figure(11); surf(FA,TT1,E)

title('HHT时频谱三维显示')www.sm136.com end

%画边际谱