做eemd的希尔伯特谱的时候为什么显示的是灰色条块状的呢？

我爱振动

用的是网上的程序，各位有知道什么原因的吗？先谢谢了啊~
程序如下

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);                        

[A,f,t]=hhspectrum(imf);        

[E,t,Cenf]=toimage(A,f);            

disp_hhs(E）
网上的蓝色图形，灰色是我自己程序的图形

hht123

画谱图的程序太简陋了吧，给你一个我的你看看，是从一本书里弄的：机械故障诊断中的现代信号处理方法
程序代码：

1. function [Spectrum,freq]=Book_HHT(Sig)                    %计算HHT的主程序
   
2. if(nargin<1),
   
3.     error('At least one input is needed');
   
4. end
   
5. Imf=emd(Sig);                           %对信号进行经验模式分解
   
6. [nImf,SigLen]=size(Imf);
   
7. t=1:SigLen;
   
8. HHSpectrum=zeros(max(1,nImf-1),SigLen);
   
9. AmpSpectrum=zeros(max(1,nImf-1),SigLen);
   
10. for k=1:nImf-1,
    
11.      Imf0=real(Imf(k,:));
    
12.      [freq,Amp]=InstFreq(Imf0);                      %用hilbert变换估计每个本征模分量的瞬时频率和功率
    
13.      HHSpectrum(k,:)=freq';
    
14.      AmpSpectrum(k,:)=Amp';
    
15. end
    
16. 
    
17. %组合所有本征模分量的瞬时频率和瞬时功率为Hilbert谱
    
18. 
    
19. MaxFreq=max(max(HHSpectrum));
    
20. MaxFreq=ceil(MaxFreq/0.5)*0.5;
    
21. NumFreq=128;
    
22. freq=linspace(0,MaxFreq,NumFreq);
    
23. Spectrum=zeros(NumFreq,SigLen);
    
24. Temp=HHSpectrum;
    
25. Temp=min(round((Temp/MaxFreq)*NumFreq)+1,NumFreq);
    
26. for k=1:nImf-1,
    
27.      for n=1:SigLen,
    
28.          Spectrum(Temp(k,n),n)=Spectrum(Temp(k,n),n)+AmpSpectrum(k,n);
    
29.      end
    
30. end
    
31. %显示信号的Hilbert谱
    
32. %  clf;
    
33. mesh(t,freq,Spectrum);
    
34. colormap jet;
    
35. axis([min(t) max(t) min(freq) max(freq)]);
    
36. shading interp;
    
37. title('Hilbert-Huang spectrum');
    
38. ylabel('Frequency');
    
39.   xlabel('Time');
    
40. %  function [Imf] =Emd (Sig)             %对信号进行经验模式分解
    
41. %  if(nargin<1),
    
42. %       error('At least one input is needed！');
    
43. %  end
    
44. % Sig=Sig';
    
45. % SigLen=length(Sig);
    
46. % t=1:SigLen;
    
47. % c=Sig;
    
48. % Imf = [ ];
    
49. % nLoop=1;
    
50. % MaxLoop=1000;
    
51. % [IndMax,IndMin]=ExtrPoint(c);                          %提取信号的极值点;极大值点和极小值点
    
52. % [Index_zer]=ZeroPoint(c);                             %提取信号的零点
    
53. % NExtr = length(IndMin) + length(IndMax);               %极值点个数
    
54. % NZero=length(Index_zer);                              %零点个数
    
55. % Bool=1;
    
56. % NEtd=2;
    
57. % IMFNum=0;
    
58. % while( NExtr > 2 && Bool > 0),
    
59. %     h=c;
    
60. %     nLoop=1;
    
61. %     SD=1;
    
62. % while((SD>0.3||(abs(NZero-NExtr)>1))&&(NExtr>2)&&nLoop<MaxLoop)    %镜像拓展序列
    
63. %   LMaxEtd=fliplr(IndMax(1:min(end,NEtd)));
    
64. %   LMinEtd=fliplr(IndMin(1:min(end,NEtd)));
    
65. %   hlMaxEtd=h(LMaxEtd);
    
66. %   hlMinEtd=h(LMinEtd);
    
67. %   LMaxEtd=-LMaxEtd+1;
    
68. %   LMinEtd=-LMinEtd+1;
    
69. %   RMaxEtd=fliplr(IndMax(max(end-NEtd+1,1):end));
    
70. %   RMinEtd=fliplr(IndMin(max(end-NEtd+1,1):end));
    
71. %   hrMaxEtd=h(RMaxEtd);
    
72. %   hrMinEtd=h(RMinEtd);
    
73. %   RMaxEtd=2*SigLen-RMaxEtd;
    
74. %   RMinEtd=2*SigLen-RMinEtd;
    
75. % %计算信号的上包络线
    
76. %   Max_Env=interp1([LMaxEtd t(IndMax) RMaxEtd],[hlMaxEtd h(IndMax) hrMaxEtd],t,'spline');
    
77. % %计算信号的下包络线
    
78. %   Min_Env=interp1([LMinEtd t(IndMin) RMinEtd],[hlMinEtd h(IndMin) hrMinEtd],t,'spline');
    
79. %   Mean_Env=(Max_Env+Min_Env)/2;                 %计算上下包络线的平均值
    
80. %   Preh=h;
    
81. %   h=h-Mean_Env;
    
82. %   alarm=1.0e-08;
    
83. %   SD=sum(((Preh-h).^2)./(Preh.^2+alarm));
    
84. %   [IndMax,IndMin]=ExtrPoint(h);                 %提取信号的极值点            
    
85. %   [Index_zer]=ZeroPoint(h);                     %提取信号的零点
    
86. %   NExtr=length(IndMin)+length(IndMax);
    
87. %   NZero=length(Index_zer);
    
88. %   nLoop=nLoop+1;
    
89. % end %提取到一个本征模分量
    
90. % IMFNum=IMFNum+1;
    
91. % disp([num2str(IMFNum),'IMFs have been obtained']);
    
92. % Imf=[Imf;h];
    
93. % c=c-h;
    
94. % [IndMax,IndMin]=ExtrPoint(c);                    %提取信号的极值点      
    
95. % [Index_zer]=ZeroPoint(c);                        %提取信号的零点
    
96. % NExtr=length(IndMin)+length(IndMax);
    
97. % NZero=length(Index_zer);
    
98. %   Bool=1;
    
99. %   if((range(c)/range(Sig))<2e-4),
    
100. %       Bool=-1;
     
101. %   end
     
102. % end
     
103. % Imf=[Imf;c];
     
104. % [n, m]=size(Imf);
     
105. % xCor=[ ];
     
106. % %
     
107. % norm_Sig=(Sig-mean(Sig))/std(Sig);
     
108. % for i=1:n-1,
     
109. %     norm_Imf=(Imf(i,:)-mean(Imf(i,:)))/std(Imf(i,:));
     
110. %     Coef=xcorr(norm_Sig,norm_Imf)/SigLen;
     
111. %     xCor=[xCor max(abs(Coef))];
     
112. % end
     
113. % k=find(xCor>max(xCor)/8);
     
114. % disp([num2str(length(k)),'IMFs have been selected.']);
     
115. % Temp=zeros(length(k)+1,SigLen);
     
116. % Temp1=Sig;
     
117. % for i=1:length(k),
     
118. %     Temp(i,:)=Imf(k(i),:);
     
119. %     Temp1=Temp1-Imf(k(i),:);
     
120. % end
     
121. % Temp(length(k)+1,:)=Temp1; %最后一个为残差分量
     
122. % Imf=Temp;
     
123. 
     
124. function[Pos_max,Pos_min]=ExtrPoint(Sig)                                    %提取信号的极大值和极小值点
     
125. 
     
126. if(nargin==0)
     
127.     error('At least one input is required');
     
128. end
     
129. SigLen=length(Sig);
     
130. DSig=diff(Sig);
     
131. DSig1=DSig(1:end-1);
     
132. DSig2=DSig(2:end);
     
133. Pos_min=find(DSig1.*DSig2<0 & DSig1<0)+1;
     
134. Pos_max=find(DSig1.*DSig2<0 & DSig1>0)+1;
     
135. Pos_max=sort([Pos_max]);
     
136. Pos_min=sort([Pos_min]);
     
137. function [PZero]=ZeroPoint(Sig)                                             %提取信号的过零点   
     
138. if(nargin==0)
     
139.     error('At least one input is required');
     
140. end
     
141. SigLen=length(Sig);
     
142. s1=Sig(1:SigLen-1);
     
143. s2=Sig(2:SigLen);
     
144. PZero=find(s1.*s2<0);
     
145. IndZero=find(Sig==0);
     
146. if(~isempty(IndZero)),
     
147.     zero=find(Sig==0);
     
148.     DZero=diff([0 zero 0]);
     
149.     LZero=find(DZero==1);
     
150.     RZero=find(DZero==-1)-1;
     
151.     IndZero=round((LZero+RZero)/2);
     
152.     PZero=sort([PZero IndZero]);
     
153. end
     
154. function[freq,Amp]=InstFreq(Sig)                                            %计算信号的瞬时频率和功率
     
155.     %  freq:信号的瞬时频率
     
156.     %  Amp :信号的瞬时功率
     
157. if(nargin==0),
     
158.     error('At least one input is required');
     
159. end;
     
160. SigLen=length(Sig);
     
161. x=hilbert(real(Sig));                                                       %计算信号的Hilbert变换
     
162. Amp=abs(x);                                                                 %计算信号的瞬时功率
     
163. freq=(angle(-x(2:SigLen).*conj(x(1:SigLen-1)))+pi)/(2*pi);                  %计算信号的瞬时功率
     
164. freq=[freq(1) freq];
     
165. std_freq=std(freq);
     
166. k=3;
     
167. threshold=k*std_freq+mean(freq);
     
168. warm= freq>threshold;
     
169. freq(warm)=mean(freq);
     
170. threshold=-k*std_freq+mean(freq);
     
171. warm= find(freq<threshold);
     
172. freq(warm)=mean(freq);     

复制代码

endless

求程序，楼主能不能把你的程序贴上来

cwb

hht123 发表于 2014-5-21 18:49
画谱图的程序太简陋了吧，给你一个我的你看看，是从一本书里弄的：机械故障诊断中的现代信号处理方法
程序 ...

怎么这么复杂啊。。。唉，顿时感觉自己学得太少了。。。。

cwb

好像用colormap设置一下就可以了吧

周2周

楼主你这个问题解决了么？

Edinburgh

colormap就行

Edinburgh

你程序中的imf=eemd(data,0.2,100);  对吗？？？？