请教HHT运行问题

棒棒糖662

我在运行程序时出现下面的提示信息，频谱图能够显示出来，（附件中是显示出的频谱图）但跟跟实际的好像不相符，刚接触这方面内容麻烦高手帮忙解释一下

??? Error using ==> instfreq
T can not be equal to 1 neither to the last element of X
Error in ==> hhspectrum at 45
  f(i,:)=instfreq(an(i,:)',tt,l)';
Error in ==> Untitled at 17
[A,f,tt]=hhspectrum(imf);
源程序如下：
fs=2000;%采样频率
N=600;%采样点数
t=0:1/fs:(N-1)*1/fs

x1=sin(2*pi*30*t);
x2=0.4*sin(2*pi*65*t+140);
x=x1+x2;
plot(t,x)[local]1[/local]
imf=emd(x);
emd_visu(x,1:length(x),imf,1)
[A,f,tt]=hhspectrum(imf);
[im,tt]=toimage(A,f);
disp_hhs(im);
colormap(flipud(gray))

zhangnan3509

emd_visu(x,1:length(x),imf,1)
这行为什么要加1？去掉可以吗

棒棒糖662

function [f2c,c2f] = emd_visu(x,t,imf,i);

if(nargin==4)
  figure(i)
else
  figure
end

其中i是figure的序列号为可选项，不要也可以

棒棒糖662

棒棒糖662

上图是运行出来的结果，不知道究竟哪里出错了

zhangnan3509

是不是EMD分解图也没有把这两个信号分开？

棒棒糖662

这是EMD分解图，这两个信号倒是能分解出来

zhangnan3509

那就不该是分解的问题了。低频信号在时频图上没看出来

棒棒糖662

??? Error using ==> instfreq
T can not be equal to 1 neither to the last element of X
Error in ==> hhspectrum at 45
  f(i,:)=instfreq(an(i,:)',tt,l)';
Error in ==> Untitled at 17
[A,f,tt]=hhspectrum(imf);
还有上面这个错误提示，我也不知道怎么回事。下面是hhspectrum的内容，instfreq程序是在网上下的，应该没问题
%inputs:
%         - imf : matrix with one IMF per row
%   - t   : time instants
%   - l   : estimation parameter for instfreq
%   - aff : if 1, displays the computation evolution
function [A,f,tt] = hhspectrum(imf,t,l,aff)
if nargin < 2

  t=1:size(imf,2);

end

if nargin < 3

  l=1;

end

if nargin < 4

  aff = 0;

end

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

  A=abs(an(:,l+1:end-l));%%%%%
  [im,tt] = toimage(A,f(i,:));
disp_hhs(im,t,-20);

  if aff

disp(['mode ',int2str(i),' trait?1?7'])
   

  end

end

棒棒糖662

顺便问一下，正常情况下这个频谱图应该是什么样的

zhangnan3509

应该在对应的频率上，是两条线。不一定很直，但是能反应出在这个频率上有信号成分。我做出的结果是这样的

[ 本帖最后由 zhangnan3509 于 2008-3-21 12:50 编辑 ]

棒棒糖662

那低频是35，归一化频率大约是0.0175，会不会是因为与0轴太接近了，所以显示不出来

zhangnan3509

你看看上面的图

棒棒糖662

那到底是我hhspectrum（）这个程序出现错误了还是instfreq有异常呢

zhangnan3509

我按照你的程序作的，除了把plot(t,x)[local]1[/local]变成把plot(t,x)，之后就是在disp_hht中纵坐标乘采样频率。除此之外，没有任何变化。