做eemd的希尔伯特谱的时候为什么显示的是灰色条块状的呢？
本帖最后由 我爱振动 于 2014-2-26 14:26 编辑用的是网上的程序，各位有知道什么原因的吗？先谢谢了啊~
程序如下
N=1000;n=1:N;fs=1000;t=n/fs;fx=10;fy=50;x=cos(2*pi*fx*t);y=10*cos(2*pi*fy*t);z=x+y;data=z;imf=eemd(data,0.2,100); =hhspectrum(imf); =toimage(A,f); disp_hhs(E）
网上的蓝色图形，灰色是我自己程序的图形
本帖最后由 牛小贱 于 2014-5-22 12:53 编辑
画谱图的程序太简陋了吧，给你一个我的你看看，是从一本书里弄的：机械故障诊断中的现代信号处理方法
程序代码：
function =Book_HHT(Sig) %计算HHT的主程序
if(nargin<1),
error('At least one input is needed');
end
Imf=emd(Sig); %对信号进行经验模式分解
=size(Imf);
t=1:SigLen;
HHSpectrum=zeros(max(1,nImf-1),SigLen);
AmpSpectrum=zeros(max(1,nImf-1),SigLen);
for k=1:nImf-1,
Imf0=real(Imf(k,:));
=InstFreq(Imf0); %用hilbert变换估计每个本征模分量的瞬时频率和功率
HHSpectrum(k,:)=freq';
AmpSpectrum(k,:)=Amp';
end
%组合所有本征模分量的瞬时频率和瞬时功率为Hilbert谱
MaxFreq=max(max(HHSpectrum));
MaxFreq=ceil(MaxFreq/0.5)*0.5;
NumFreq=128;
freq=linspace(0,MaxFreq,NumFreq);
Spectrum=zeros(NumFreq,SigLen);
Temp=HHSpectrum;
Temp=min(round((Temp/MaxFreq)*NumFreq)+1,NumFreq);
for k=1:nImf-1,
for n=1:SigLen,
Spectrum(Temp(k,n),n)=Spectrum(Temp(k,n),n)+AmpSpectrum(k,n);
end
end
%显示信号的Hilbert谱
%clf;
mesh(t,freq,Spectrum);
colormap jet;
axis();
shading interp;
title('Hilbert-Huang spectrum');
ylabel('Frequency');
xlabel('Time');
%function =Emd (Sig) %对信号进行经验模式分解
%if(nargin<1),
% error('At least one input is needed！');
%end
% Sig=Sig';
% SigLen=length(Sig);
% t=1:SigLen;
% c=Sig;
% Imf = [ ];
% nLoop=1;
% MaxLoop=1000;
% =ExtrPoint(c); %提取信号的极值点;极大值点和极小值点
% =ZeroPoint(c); %提取信号的零点
% NExtr = length(IndMin) + length(IndMax); %极值点个数
% NZero=length(Index_zer); %零点个数
% Bool=1;
% NEtd=2;
% IMFNum=0;
% while( NExtr > 2 && Bool > 0),
% h=c;
% nLoop=1;
% SD=1;
% while((SD>0.3||(abs(NZero-NExtr)>1))&&(NExtr>2)&&nLoop<MaxLoop) %镜像拓展序列
% LMaxEtd=fliplr(IndMax(1:min(end,NEtd)));
% LMinEtd=fliplr(IndMin(1:min(end,NEtd)));
% hlMaxEtd=h(LMaxEtd);
% hlMinEtd=h(LMinEtd);
% LMaxEtd=-LMaxEtd+1;
% LMinEtd=-LMinEtd+1;
% RMaxEtd=fliplr(IndMax(max(end-NEtd+1,1):end));
% RMinEtd=fliplr(IndMin(max(end-NEtd+1,1):end));
% hrMaxEtd=h(RMaxEtd);
% hrMinEtd=h(RMinEtd);
% RMaxEtd=2*SigLen-RMaxEtd;
% RMinEtd=2*SigLen-RMinEtd;
% %计算信号的上包络线
% Max_Env=interp1(,,t,'spline');
% %计算信号的下包络线
% Min_Env=interp1(,,t,'spline');
% Mean_Env=(Max_Env+Min_Env)/2; %计算上下包络线的平均值
% Preh=h;
% h=h-Mean_Env;
% alarm=1.0e-08;
% SD=sum(((Preh-h).^2)./(Preh.^2+alarm));
% =ExtrPoint(h); %提取信号的极值点
% =ZeroPoint(h); %提取信号的零点
% NExtr=length(IndMin)+length(IndMax);
% NZero=length(Index_zer);
% nLoop=nLoop+1;
% end %提取到一个本征模分量
% IMFNum=IMFNum+1;
% disp();
% Imf=;
% c=c-h;
% =ExtrPoint(c); %提取信号的极值点
% =ZeroPoint(c); %提取信号的零点
% NExtr=length(IndMin)+length(IndMax);
% NZero=length(Index_zer);
% Bool=1;
% if((range(c)/range(Sig))<2e-4),
% Bool=-1;
% end
% end
% Imf=;
% =size(Imf);
% xCor=[ ];
% %
% norm_Sig=(Sig-mean(Sig))/std(Sig);
% for i=1:n-1,
% norm_Imf=(Imf(i,:)-mean(Imf(i,:)))/std(Imf(i,:));
% Coef=xcorr(norm_Sig,norm_Imf)/SigLen;
% xCor=;
% end
% k=find(xCor>max(xCor)/8);
% disp();
% Temp=zeros(length(k)+1,SigLen);
% Temp1=Sig;
% for i=1:length(k),
% Temp(i,:)=Imf(k(i),:);
% Temp1=Temp1-Imf(k(i),:);
% end
% Temp(length(k)+1,:)=Temp1; %最后一个为残差分量
% Imf=Temp;
function=ExtrPoint(Sig) %提取信号的极大值和极小值点
if(nargin==0)
error('At least one input is required');
end
SigLen=length(Sig);
DSig=diff(Sig);
DSig1=DSig(1:end-1);
DSig2=DSig(2:end);
Pos_min=find(DSig1.*DSig2<0 & DSig1<0)+1;
Pos_max=find(DSig1.*DSig2<0 & DSig1>0)+1;
Pos_max=sort();
Pos_min=sort();
function =ZeroPoint(Sig) %提取信号的过零点
if(nargin==0)
error('At least one input is required');
end
SigLen=length(Sig);
s1=Sig(1:SigLen-1);
s2=Sig(2:SigLen);
PZero=find(s1.*s2<0);
IndZero=find(Sig==0);
if(~isempty(IndZero)),
zero=find(Sig==0);
DZero=diff();
LZero=find(DZero==1);
RZero=find(DZero==-1)-1;
IndZero=round((LZero+RZero)/2);
PZero=sort();
end
function=InstFreq(Sig) %计算信号的瞬时频率和功率
%freq:信号的瞬时频率
%Amp :信号的瞬时功率
if(nargin==0),
error('At least one input is required');
end;
SigLen=length(Sig);
x=hilbert(real(Sig)); %计算信号的Hilbert变换
Amp=abs(x); %计算信号的瞬时功率
freq=(angle(-x(2:SigLen).*conj(x(1:SigLen-1)))+pi)/(2*pi); %计算信号的瞬时功率
freq=;
std_freq=std(freq);
k=3;
threshold=k*std_freq+mean(freq);
warm= freq>threshold;
freq(warm)=mean(freq);
threshold=-k*std_freq+mean(freq);
warm= find(freq<threshold);
freq(warm)=mean(freq);
求程序，楼主能不能把你的程序贴上来 hht123 发表于 2014-5-21 18:49
画谱图的程序太简陋了吧，给你一个我的你看看，是从一本书里弄的：机械故障诊断中的现代信号处理方法
程序 ...
怎么这么复杂啊。。。唉，顿时感觉自己学得太少了。。。。 好像用colormap设置一下就可以了吧
楼主你这个问题解决了么？ colormap就行 你程序中的imf=eemd(data,0.2,100);对吗？？？？
页:
[1]