声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 13986|回复: 60

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

  [复制链接]
发表于 2012-4-16 17:17 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
在HHT实际振动信号分析中,一旦信号中混有噪声,HHT谱图显示效果非常不好。甚至根本看不出时频特性。
下面给出一个实例,希望大侠们能给出解决方案,谢谢!。(使用的是论坛下载的EMD、hhspectrum、toimage、disp_hhs来生成Hilbert谱)
数据:使用的是http://forum.vibunion.com/thread-105573-1-1.html中的数据(截取前20000个数据)

时间幅值图

时间幅值图

频率幅值图

频率幅值图

IMF1/3

IMF1/3

IMF2/3

IMF2/3

IMF3/3

IMF3/3

Hilbert谱

Hilbert谱

边界谱

边界谱

评分

1

查看全部评分

回复
分享到:

使用道具 举报

发表于 2012-4-16 20:30 | 显示全部楼层
好厉害啊 这些图我都不知道是怎么搞出来的
发表于 2012-4-17 06:04 | 显示全部楼层
用EEMD是否会效果好一些,没试过,有空试试比较一下
 楼主| 发表于 2012-4-17 19:56 | 显示全部楼层
回复 3 # syxqq123 的帖子

这是用EEMD做出来的,Hilbert谱效果还是不怎么好。
EEMD_IMFs1/3:

EEMD_IMFs3_1

EEMD_IMFs3_1

EEMD_IMFs2/3:

EEMD_IMFs3_2

EEMD_IMFs3_2

EEMD_IMFs3/3:

EEMD_IMFs3_3

EEMD_IMFs3_3

EEMD_Hilbert谱:

EEMD_Hilbert谱

EEMD_Hilbert谱

EEMD_边界谱:

EEMD_边界谱

EEMD_边界谱


 楼主| 发表于 2012-4-17 20:18 | 显示全部楼层
回复 3 # syxqq123 的帖子

另外,使用诸福磊等编写《机械故障诊断中的现代信号处理方法》中关于Hilbert谱matlab程序(HHT、InstFreq函数)做出的Hilbert谱图如下图,感觉效果要好一些(相对使用hhspectrum、toimage、disp_hhs函数来生成的Hilbert谱)。
后面是《机械故障诊断中的现代信号处理方法》中关于Hilbert谱matlab程序,现在正在对比这两中方法,看看到底是算法本身的差异,还是作图的差异。望老师也帮忙看看,多谢!
HHT、InstFreq函数做出的HIlbert谱图:
新Hilbert谱.jpg
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');

评分

1

查看全部评分

 楼主| 发表于 2012-4-17 20:27 | 显示全部楼层
回复 3 # syxqq123 的帖子

还有边界谱幅值最大值接近14000,这是怎么回事呢?
发表于 2012-4-18 17:30 | 显示全部楼层
回复 6 # 梦泉 的帖子

不好意思,我也只了解一点,不知道你的问题出在哪里,呵呵
发表于 2012-4-18 17:32 | 显示全部楼层
本帖最后由 syxqq123 于 2012-4-18 17:33 编辑

最近也一直在学习LMD,初步感觉比HHT要好一些,但是程序很难编
到现在还没搞出来
 楼主| 发表于 2012-4-18 22:25 | 显示全部楼层
回复 8 # syxqq123 的帖子

哈,还是谢谢啦!
下面的网址是从论坛发现关于LMD源程序网址,不知你发现没。
http://blog.sina.com.cn/s/blog_574d08530100r1yw.html
发表于 2012-4-19 17:44 | 显示全部楼层
回复 9 # 梦泉 的帖子

感谢,发现了,有一定借鉴意义,
发表于 2012-4-19 21:09 | 显示全部楼层
不知道我的函数问题还是软件工具箱的问题,程序正常,可惜就是结果出不来!??? 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);
这是我的错误提示。。。
发表于 2012-4-23 09:04 | 显示全部楼层
学习时间不长,看不懂这些东西啊
 楼主| 发表于 2012-4-24 14:48 | 显示全部楼层
回复 11 # znas0707 的帖子

把你写的程序贴出来看看
发表于 2012-4-24 19:17 | 显示全部楼层
这个东西我倒是做出来了,可是不理解那个HIlbert谱图所表达的含义啊,请教一下?谢谢!
 楼主| 发表于 2012-4-24 22:14 | 显示全部楼层
回复 14 # 176045173 的帖子

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

评分

1

查看全部评分

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-24 23:41 , Processed in 0.070157 second(s), 23 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表