您好,欢迎访问三七文档
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);fori=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);%总的方差fori=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);('时间(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);('时间(t)');ylabel('瞬时频率');title('imf2')%计算HHT的时频谱[A,fa,tt]=hhspectrum(imf);ifsize(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谱fori=1:m-1faa=fa(i,:);[FA,TT1]=meshgrid(faa,tt1);%三维图显示HHT时频图figure(11);surf(FA,TT1,E)title('HHT时频谱三维显示')画边际谱E=flipud(E);=1:size(E,1)bjp(k)=sum(E(k,:))*1/fs;endf=(0:N-3)/N*(fs/2);figure(12)plot(f,bjp);
本文标题:hht变换
链接地址:https://www.777doc.com/doc-2876249 .html