HHT在实际振动信号分析中的问题

梦泉

在HHT实际振动信号分析中，一旦信号中混有噪声，HHT谱图显示效果非常不好。甚至根本看不出时频特性。
下面给出一个实例，希望大侠们能给出解决方案，谢谢！。（使用的是论坛下载的ＥＭＤ、hhspectrum、toimage、disp_hhs来生成Ｈｉｌｂｅｒｔ谱）
数据：使用的是http://forum.vibunion.com/thread-105573-1-1.html中的数据（截取前20000个数据）

zhouwenjun

好厉害啊 这些图我都不知道是怎么搞出来的

syxqq123

用EEMD是否会效果好一些，没试过，有空试试比较一下

梦泉

回复 3 # syxqq123 的帖子

这是用ＥＥＭＤ做出来的，Ｈｉｌｂｅｒｔ谱效果还是不怎么好。
EEMD_IMFs１／３：

EEMD_IMFs3_1

EEMD_IMFs２／３：

EEMD_IMFs3_2

EEMD_IMFs３／３：

EEMD_IMFs3_3

ＥＥＭＤ＿Ｈｉｌｂｅｒｔ谱：

EEMD_Hilbert谱

EEMD_边界谱：

EEMD_边界谱

梦泉

回复 3 # syxqq123 的帖子

另外，使用诸福磊等编写《机械故障诊断中的现代信号处理方法》中关于Hilbert谱matlab程序（HHT、InstFreq函数）做出的Hilbert谱图如下图，感觉效果要好一些（相对使用hhspectrum、toimage、disp_hhs函数来生成的Hilbert谱）。
后面是《机械故障诊断中的现代信号处理方法》中关于Hilbert谱matlab程序，现在正在对比这两中方法，看看到底是算法本身的差异，还是作图的差异。望老师也帮忙看看，多谢！
HHT、InstFreq函数做出的HIlbert谱图：

HHT函数：
function [Spectrum,freq]=HHT(Imf)
if (nargin<1)
    error('至少需要一个输入！');
end
[nImf,SigLen]=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,:));
    [freq,Amp]=InstFreq(Imf0);
    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

InstFreq函数：
function [freq,Amp]=InstFreq(Sig);
if (nargin==0)
    error('At least one input is required');
end
SigLen=length(Sig);
x=hilbert(real(Sig));
Amp=abs(x);
freq=(angle(-x(2:SigLen).*conj(x(1:SigLen-1)))+pi)/(2*pi);
freq=[freq(1) freq];
std_freq=std(freq);
k=3;
threshold=k*std_freq+mean(freq);
warm=find(freq>threshold);
freq(warm)=mean(freq);
threshold=-k*std_freq+mean(freq);
warm=find(freq<threshold);
freq(warm)=mean(freq);

Hilbert谱的绘制：
figure('numbertitle','off','name','Hilbert幅值谱');
mesh(t,freq,Spectrum);
colormap jet;
axis([min(t),max(t),min(freq),max(freq)]);
set(gca,'YLim',[min(freq),max(freq)]);
shading interp;
title('Hilbert-Huang spectrum');
ylabel('Frequency/Fs');
xlabel('Time');

梦泉

回复 3 # syxqq123 的帖子

还有边界谱幅值最大值接近14000，这是怎么回事呢？

syxqq123

回复 6 # 梦泉 的帖子

不好意思，我也只了解一点，不知道你的问题出在哪里，呵呵

syxqq123

最近也一直在学习LMD，初步感觉比HHT要好一些，但是程序很难编
到现在还没搞出来

梦泉

回复 8 # syxqq123 的帖子

哈，还是谢谢啦！
下面的网址是从论坛发现关于LMD源程序网址，不知你发现没。
http://blog.sina.com.cn/s/blog_574d08530100r1yw.html

syxqq123

回复 9 # 梦泉 的帖子

感谢，发现了，有一定借鉴意义，

znas0707

不知道我的函数问题还是软件工具箱的问题，程序正常，可惜就是结果出不来！??? Error using ==> chckxy at 51
The data sites should be distinct.
Error in ==> spline at 55
[x,y,sizey,endslopes] = chckxy(x,y);
Error in ==> emd>getspline at 54
s = spline([0 p N+1],[0 x(p) 0],1:N);
Error in ==> emd at 13
      s1 = getspline(x1);
Error in ==> Untitled38 at 37
imf = emd(y);
这是我的错误提示。。。

twb0624

学习时间不长，看不懂这些东西啊

梦泉

回复 11 # znas0707 的帖子

把你写的程序贴出来看看

176045173

这个东西我倒是做出来了，可是不理解那个HIlbert谱图所表达的含义啊，请教一下？谢谢!

梦泉

回复 14 # 176045173 的帖子

copy论坛的一个表达，不知能否回答你提的问题。
做Hilbert谱的过程是这样的，通过Hilber变换将信号转换成分析信号，然后分解信号的幅值和角度，利用求解瞬时频率的.m文件求解角度的导数（亦即瞬时频率）。注意这里求出来的瞬时频率是归一化的，因为调用时频工具箱中求解瞬时频率的函数是时间项变了（有1/采样频率变成了1）。接着就是利用提取出来的频率和幅值做成三维显示（横坐标时间，纵坐标频率，颜色深浅显示幅值大小）这就是所谓的Hilbert谱。