调频信号，HHT和fft哪个正确？

zhangnan3509

这到底为什么，我都有点糊涂了！是采样时间的原因吗？

[ 本帖最后由 zhangnan3509 于 2007-5-7 22:07 编辑 ]

rc-hw-0002

原帖由 shenyongjun 于 2007-5-7 17:11 发表
采样时间好像有问题吧？

恩,错了,没改过来.'

吃书的老虎

我的理解好像和你不相同，楼主能发一篇这样的论文名称给我看看吗？

zhangnan3509

《基于HHT的旋转机械故障诊断方法研究》作者胡劲松，杨世锡。该文章的图6效果和楼主的类似，仍然是两端比较高，中间凹下去了，我想这也不能说FFT是错的，毕竟分析的角度不一样，而且我觉得HHT的瞬时频率定义是唯一的，也是独特的，在FFT里面就没有这样的定义，就不能拿来类比说对错。

[ 本帖最后由 zhangnan3509 于 2007-5-8 14:33 编辑 ]

zhangnan3509

楼主，我好好想了一下这个问题，我还是觉得这个无所谓对错，只是算法的问题！FFT可能错了吗？只用这一个例子就推翻了整个傅立叶变换吗？这不现实！是边际谱错了吗？这也不是，边际谱不过是对HHT谱时间的积分，这也不会错！
  我想错也就错在就不能在这个问题上拿来对比！
  楼主应该注意这两幅图的纵坐标不一样，代表的意义就不相同，它们怎么比中间的形状呢？边际谱从统计意义上表示了整组数据每个频率点的积累幅值分布，而傅立叶谱的某点幅值表示在整个信号里有一个含有此频率的三角函数组分，而且幅值越大只是说明在整个数据段上，局部存在的可能性越大。
   再看得到的图形，FFT表示的是整个数据中，能量在一个频率上分布的可能性地描述，而边际谱表示在在每一个频率上幅值的积累，如果想知道具体时间那么就看HHT谱，这个时间-幅值-频率的三维谱。这两个意义并不一样。
  说到瞬时频率，傅立叶变换不强调局部性，而是强调全局性。咱们的HHT才提出一个唯一的瞬时频率的定义。因此拿瞬时频率来衡量傅立叶变换也是不公平的。
   所以我说，这里面没有对错。
   有些文章里面也有FFT谱和边际谱的比较，但是比较的内容和楼主的完全不同。在上面提到的胡劲松教授的文章中，比较的是频率分布情况和在实际故障处理中，那一个更明显的反映故障特征，而不是谱图线的走法。
    我处理以前转子台数据中也不止一次出现过类似问题，但是我没往这个角度去想。
   以上是我个人的一点小看法。希望大家指正

duguzi

HHT边际谱中的频率与傅立叶分析中的频率意义完全不同。在傅立叶表达中，在某一频率W处能量的存在，代表一个正弦或余弦波在整个时间长度上都存在。HHT中，在某一频率W处能量的存在，仅代表在数据的整个时间长度上，很可能有这样一个频率的振动波在局部出现过。

因为HHT和FFT中频率的定义不同，出现上述差异是正常的。

HHT边际谱是HHT时频谱对时间的积分，出现“为什么HHT边际谱调频段(即x2)的最大幅值分布在两端？”这一问题，可以借助正弦信号概率分布密度函数出现类似情况来理解。

如有甚么问题，可以发信给我：xyg_1975@163.com

我很少上网，一般mail都会回复的

lengrong

这是因为没有中心频率的原因。请加上一个中心频率成份再分析一下。

williamzhaos

我现在已经得到了最后分解得到的IMF分量，而且应该还是比较准确的，但是在求时频图和边际谱上却存在着很大的奇异，那位高手能帮忙解决一下，不胜感激！

zhangnan3509

把图贴上来看看，这样更直观一点

williamzhaos

t=[0:1/10000:1];s=cos(10*pi*t)+2*cos(40*pi*t);imf分量图和时频图如下，时频图中的频率对应的还是对的，第一幅对应的是20Hz，第二幅是5Hz。

[ 本帖最后由 williamzhaos 于 2007-11-8 10:04 编辑 ]

zhangnan3509

你的时频图是怎么做的？

williamzhaos

l=length(s);
a=size(c);
b=a(1);
for m=1:b
    H(m,:)=hilbert(c(m,:));
    x(m,:)=real(H(m,:));
    y(m,:)=imag(H(m,:));
    fi(m,:)=atan(y(m,:)./x(m,:));
    g(m,:)=szwf(fi(m,:),t);
    f(m,:)=g(m,:)./(2*pi);
    subplot(b,1,m);
    plot(t,f(m,:));
    xlabel('t(s)');
    ylabel('freq(hz)')
    title('时频图');
end
这个是做完emd分解后，根据IMF分量作出的，就是根据我看得论文中时频图的做法，求出瞬时频率做出的。
流程是先进行hilbert变换（直接调用的），再求其虚部和实部之比的反正切，再求数值微分（这个程序是自己编的，经验证还算可以），最后得出时频图。

[ 本帖最后由 williamzhaos 于 2007-11-8 10:23 编辑 ]

zhangnan3509

t=1:500;
t=t*1/400;
x1=sin(2*pi*20*t);
x2=0.4*sin(2*pi*40*t+140);
z=x1+x2
plot(t,z)
y= linspace(1,2,500);%添加趋势项
z=z+y
imf=emd(z);
emd_visu(z,1:length(z),imf)
[A,f,tt]=hhspectrum(imf);
[im,tt]=toimage(A,f);
disp_hhs(im);
colormap(flipud(gray))

rilling的程序里面有这样函数，你可以改一下前面的要处理的模拟信号，就能得出时频图了

[ 本帖最后由 zhangnan3509 于 2007-11-8 10:35 编辑 ]

yangzj

是不是由于求出的相位没经过处理直接求导
而且求相位应该用atan2

[ 本帖最后由 yangzj 于 2007-11-8 10:40 编辑 ]

williamzhaos

谢谢高手，其实我要处理的东西应该是没有大家说的麻烦，我编这个程序的目的是为了检测实际信号中的谐波信号，其实一直都是数字信号，我给出的初始信号其实只是为了验证程序的正确性，其实不是一个模拟信号。另外，我所编的程序全部都是用的function编的一个一个子程序然后套用的，加上本人又是一个编程的初学者，有些东西还不是太明白，希望高手能根据这个说的再明白些，麻烦了！:handshake