|
楼主 |
发表于 2013-8-7 11:27
|
显示全部楼层
clc;clear;
pass='D:\yiqulb\';
file='a1-0.5-15.txt';
z=load([pass,file]);
z=z(1:2:6000,:);
data=z';
imf=emd(data);
fs=50;
t=1:size(imf,2);
l=1;
aff = 0;
lt=length(t);
tt=t((l+1):(lt-l));
for i=1:(size(imf,1)-1)
an(i,:)=hilbert(imf(i,:)')';
f(i,:)=instfreq(an(i,:)',tt,l)';%Subscript indices must either be real positive integers or logicals.
A=abs(an(:,l+1:end-l));
if aff
disp(['mode ',int2str(i),' trait?'])
end
end
%对输入信号进行EMD分解
[A,f,tt]=hhspectrum(imf,t,l,aff); %对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); %希尔伯特谱----------------------------------------------------------
%画出边际谱
%N=length(Cenf);%设置频率点数 %完全从理论公式出发。网格化后中心频率很重要,大家从连续数据变为离散的角度去思考,相信应该很容易理解
for k=1:size(E,1)
bjp(k)=sum(E(k,:))*1/fs;
end
H=size(E,1);
f0=(0:H-1)/H*(fs/2);
figure(3);
plot(f0,bjp); % 作边际谱图 进行求取Hilbert谱时频率已经被抽样成具有一定窗长的离散频率,所以此时的频率轴已经是中心频率
xlabel('??/ Hz');
ylabel('??');
[url=]\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\[/url]
以上是我用的程序,采用降低采样频率的方法取了其中3000个数据,针对每个IMF分量画出了其时频图和边际谱,共有10个IMF分量,只取了前七个,得到的图如下:
|
|