Hilbert边际谱程序（matlab）

zhlong

原帖由 zhangnan3509 于 2007-7-11 08:13 发表

                               
登录/注册后可看大图

那太好了！我觉得HHT用在这个数据上不如小波好，想听听你的看法

以前没处理过轴承信号，做的都是转子信号，所以一直没去做这个例子。
版主可以参考一下这篇文献：
《小波变换及经验模式分解方法在电机轴承早期故障诊断中的应用》

zhangnan3509

我只是对于这个数据来说的，看看效果到底怎么样

form

原帖由 zhlong 于 2007-7-11 07:34 发表

                               
登录/注册后可看大图

破凰的这段程序没有问题的，频率序列也可以和边际谱点数对应起来！

两位版主起的都很早啊，呵呵。我想把频率与边际谱对应起来。取了2048点，结果bjp算出来只有400了。模拟信号好像容易对应些

zhlong

边际谱的点数取决于

[E,tt1]=toimage(A,fa,tt,length(tt));

这里面参数的设置。

zhangnan3509

看来还是要研究啊 原以为会告一段落，没想到这是刚刚开始

serendipity_w

我运行破凰的程序怎么提示出错啊？
??? Undefined function or variable 'isscalar'.

Error in ==> C:\MATLAB6p5\work\HHT分析完整版\toimage.m
On line 50  ==>     if isscalar(varargin{1})

chiu

clear;
fs=1000;    %fs为采样频率;
N=1000;     %采样点数
t=1/fs:1/fs:1;
y1=2*sin(60*pi*t);
y2=5*sin(90*pi*t);
y=[y1;y2;zeros(size(y1))]; %IMF集
%%%%%%%%%%%%%求边际谱
[A,fa,tt]=hhspectrum(y);
[E,tt1]=toimage(A,fa,tt);
% E=flipud(E);
for k=1:size(E,1)
    bjp(k)=sum(E(k,:))*1/fs;  
end
%f=(0:N-3)/N*(fs/2);
r=size(E,1);
f=((1:r)/(2*r))*(fs);
plot(f,bjp);
xlabel('频率 / Hz');
ylabel('幅值');
這是小弟對畫圖座標部份做的一些小小更改,不知是否正確,希望與大家討論,並得到指正
toimage的部份frequency resolution仍然維持預設的400(因為小弟分析的數據都是好幾萬點,如果再提高的話,電腦會跑不動),所以frequency resolution會比原來的998稍差一點,由於小弟第一次發文,還
不知道怎麼附圖請見諒
另外 想請問bjp(k)=sum(E(k,:))*1/fs;中乘上1/fs是為甚麼,小弟一時想不到,請大家幫忙了 謝謝

zhlong

isscalar是matlab自带函数，可能你的matlab太老。

%ISSCALAR True if array is a scalar.
%   ISSCALAR(S) returns logical true (1) if S is a 1 x 1 matrix
%   and logical false (0) otherwise.
%
%   See also ISVECTOR, ISNUMERIC, ISLOGICAL, ISCHAR, ISEMPTY.

%   Copyright 1984-2005 The MathWorks, Inc.
%   $Revision: 1.1.6.6 $  $Date: 2005/06/21 19:28:31 $

mofei

我直接使用了破凰的程序：
clear;
fs=1000;    %fs为采样频率;
N=1000;     %采样点数
t=1/fs:1/fs:1;
y1=2*sin(60*pi*t);
y2=5*sin(90*pi*t);
y=[y1;y2;zeros(size(y1))]; %IMF集
%%%%%%%%%%%%%求边际谱
[A,fa,tt]=hhspectrum(y);
[E,tt1]=toimage(A,fa,tt,length(tt));
E=flipud(E);
for k=1:size(E,1)
    bjp(k)=sum(E(k,:))*1/fs;  
end
f=(0:N-3)/N*(fs/2);
plot(f,bjp);
xlabel('频率 / Hz');
ylabel('幅值');
但算出的结果如下图，这是怎么回事呢？

zhangnan3509

你的EMD程序是哪一版的？

mofei

用到的函数都是从法国的Grilling提供的开源程序中下载的。
下载地址：http://perso.ens-lyon.fr/patrick.flandrin/emd.html

zhlong

E=flipud(E);  去掉这一句

zhangnan3509

你现在用的一定不是新的啊，前面的讨论说起过这个问题

mofei

mofei

EMD算法是破凰推荐的网站下载的。难道是matlab的版本的问题吗？