HHT_hht_HHT边际谱_HHT谱_边际谱_希尔伯特变换

clc; clear all; close all; N= 1024; fs=1000000; t=1/fs:1/fs:N/fs; s2=xlsread('SY-3_3_129858.csv' ,'al3:a1036' ); % s2=xlsread('SY-3_3_129880.csv' ,'al3:a1036' ); s2=wden(s2,'heursure' ,'s' ,'one' ,3,'db8' ); data=s2'; imf=emd(data); % 对输入信号进行EMD分解 % X1=imf(1,:); % X2=imf(2,:); % X3=imf(3,:); % X4=imf(4,:); % X=X1+X2+X3+X4; % [A,f,t]=hhspectrum(X); % [E,t,Cenf]=toimage(A,f); [A,f,t]=hhspectrum(imf); % 对IMF分量求取瞬时频率与振幅：A：是每个IMF的振幅向量,f:每个IMF对应的瞬时频率，t：时间序列号 [E,t,Cenf]=toimage(A,f); % 将每个IMF信号合成求取Hilbert谱，E：对应的振幅值，Cenf：每个网格对应的中心频率。这里横轴为时间，纵轴为频率。即时频图(用颜色表示第三维值的大小)和三维图（三维坐标系：时间，中心频率，振幅） % cemd_visu(data,1:length(data),imf); % 显示每个IMF分量及残余信号 disp_hhs(E); % 希尔伯特谱 %% 画出边际谱 % colormap(flipud(gray)); % 黑白显示 %N=length(Cenf); % 设置频率点数。完全从理论公式出发。网格化后中心频率很重要，大家从连续数据变为离散的角度去思考，相信应该很容易理解 for k=1:size(E,1) bjp(k)=sum(E(k,:))*1/fs; end figure; plot(Cenf(1,:)*fs/1000,bjp); % 作边际谱图 进行求取Hilbert谱时频率已经被抽样成具有一定窗长的离散频率，所以此时的频率轴已经是中心频率 grid on; xlabel('频率 / kHz'); ylabel('幅值'); axis([0 150 0 8*10^(-6)]) % 绘制瞬时包络图和瞬时频率图 % figure; % subplot(221),plot(t/N,A(1,:));xlabel('时间 t/s');ylabel('幅值');title('imf1分量瞬时包络'); % subplot(222),plot(t/N,f(1,:)*fs);xlabel('时间 t/s');ylabel('频率');title('imf1分量瞬时频率'); % subplot(223),plot(t/N,A(2,:));xlabel('时间 t/s');ylabel('幅值');title('imf2分量瞬时包络'); % subplot(224),plot(t/N,f(2,:)*fs);xlabel('时间 t/s');ylabel('频率');title('imf2分量瞬时频率');