下面将程序列出来,希望有大神能指导下。- clear all;
- clc;
- %示例程序
- N=1024;
- n=0:N-1;
- fs=12800;
- t=n/fs;
- x=(sin(2*pi*50*t)+sin(2*pi*150*t)).*(t>=0&t<0.5)+(sin(2*pi*50*t)+sin(2*pi*250*t)).*(t>=0.5&t<=1);
- data=x;
- imf=emd(data); %对输入信号进行EMD分解
- [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(3);
- plot(Cenf(1,:)*fs,bjp); % 作边际谱图 进行求取Hilbert谱时频率已经被抽样成具有一定窗长的离散频率,所以此时的频率轴已经是中心频率
- xlabel('频率 / Hz');
- ylabel('幅值');
- %绘制瞬时包络图和瞬时频率图
- 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分量瞬时频率');
复制代码 下面是仿真的结果,还有我理想的结果(黑白的)。
求大神指教!
|