声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

楼主: zhoupingwoo

[HHT] Hilbert边际谱程序(matlab)

  [复制链接]
发表于 2007-4-8 20:17 | 显示全部楼层
原帖由 form 于 2007-3-8 17:08 发表
我是用hspec画谱图的,自己又在后面添加了求边际谱的代码,但感觉有问题
for k=1:size(E)
    bjp(k)=sum(E(k,:))*1/fs;  %fs为采样频率;
end
figure
plot(bjp);
xlabel('频率 / Hz');
ylabel('幅值');
...



前面这是我做的,你是在hhspectrum后面加的吗?
回复 支持 反对
分享到:

使用道具 举报

发表于 2007-4-8 21:45 | 显示全部楼层
我不是这么做的!我用的是nspab.m和mspc.m这两个做的。你也应该有吧

[ 本帖最后由 zhangnan3509 于 2007-4-8 21:50 编辑 ]
回复 支持 0 反对 1

使用道具 举报

发表于 2007-4-10 12:30 | 显示全部楼层
原帖由 zhangnan3509 于 2007-4-8 21:45 发表
我不是这么做的!我用的是nspab.m和mspc.m这两个做的。你也应该有吧

你两个函数是自己编的吗?我下载的程序包没有这两个函数啊?
我第二次作的是在hhspectrum后面添加的,感觉按描述的应该是在这个函数基础上来作,那你做的图和我的比较如何?

[ 本帖最后由 zhlong 于 2007-6-1 20:12 编辑 ]
发表于 2007-4-10 15:30 | 显示全部楼层

回复 #34 form 的帖子

你自己再用这个试试 我没有用rilling的hhspectrum
nspab:
function h1= nspab(data,nyy,minw,maxw,dt)

% The function NSPAB generates a smoothed HHT spectrum of data(n,k)  
% in time-frequency space, where  
% n specifies the length of time series, and
% k is the number of IMF components.
% The frequency-axis range is prefixed.
% Negative frequency sign is reversed.
%
% MATLAB Library function HILBERT is used to calculate the Hilbert transform.
%
% Example, [h,xs,w] = nspab(lod78_p',200,0,0.12,1,3224).
%
% Functions CONTOUR or IMG can be used to view the spectrum,
%    for example contour(xs,w,h) or img(xs,w,h).
%
% Calling sequence-
% [h,xs,w] = nspab(data,nyy,minw,maxw,t0,t1)
%
% Input-
%        data        - 2-D matrix data(n,k) of IMF components
%        nyy        - the frequency resolution
%        minw        - the minimum frequency
%        maxw        - the maximum frequency
%        t0        - the start time
%        t1        - the end time
% Output-
%        h        - 2-D matrix of the HHT spectrum, where
%                  the 1st dimension specifies the number of frequencies,
%                  the 2nd dimension specifies the number of time values
%        xs        - vector that specifies the time-axis values
%        w        - vector that specifies the frequency-axis values

% Z. Shen (JHU)                July 2, 1995 Initial

%----- Get dimensions (number of time points and components)
[npt,knb] = size(data);

%----- Get time interval


%----- Apply Hilbert Transform
data=hilbert(data);
a=abs(data);
omg=abs(diff(unwrap(angle(data))))/(2*pi*dt);

%----- Smooth amplitude and frequency   
filtr=fir1(8,.1);
for i=1:knb  
    a(:,i)=filtfilt(filtr,1,a(:,i));
   omg(:,i)=filtfilt(filtr,1,omg(:,i));
end


%----- Limit frequency and amplitude
for i=1:knb
   for i1=1:npt-1
      if omg(i1,i) >=maxw,
         omg(i1,i)=maxw;
         a(i1,i)=0;
      elseif omg(i1,i)<=minw,
         omg(i1,i)=minw;
         a(i1,i)=0;
      else
      end
   end
end
clear filtr data

%va=var(omg(200:1200))

%----- Get local frequency
dw=maxw - minw;
wmx=maxw;
wmn=minw;

%----- Construct the ploting matrix
clear p;
h1=zeros(npt-1,nyy+1);
p=round(nyy*(omg-wmn)/dw)+1;
for j1=1:npt-1
   for i1=1:knb
      ii1=p(j1,i1);       
      h1(j1,ii1)=h1(j1,ii1)+a(j1,i1);
   end
end

%----- Do 3-point to 1-point averaging
[nx,ny]=size(h1);
%n1=fix(nx/3);
%h=zeros(n1,ny);
%for i1=1:n1
   %h(i1,:)=(h1(3*i1,:)+h1(3*i1-1,:)+h1(3*i1-2,:));
%end
%clear h1;

%----- Do 3-points smoothing in x-direction
fltr=1./3*ones(3,1);
for j1=1:ny
  h1(:,j1)=filtfilt(fltr,1,h1(:,j1));
end
clear fltr;


%----- Define the results
%w=linspace(wmn,wmx,ny-1)';
%xs=linspace(t0,t1,nx)';
h1=flipud(rot90(h1));
h1=h1(1:ny-1,:);

评分

1

查看全部评分

发表于 2007-4-11 22:38 | 显示全部楼层
请问HHT边际谱和FFT谱有什么区别,差别在哪里?:loveliness:
发表于 2007-4-22 21:53 | 显示全部楼层
看了上面的程序,真是有点思路!我记得小波分析功率谱图就是这样画的
发表于 2007-4-27 21:51 | 显示全部楼层
form求边际谱时所用程序是没有问题的,用的是矩形积分公式。
他所得结果不正确的原因是:输入的应是调用了toimage后的结果,而不是调用了hhspectrum后的结果。
下面给一段程序,大家可以去试下。边际谱的分析结果是完全正确的。
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('幅值');

%(完整答案)

[ 本帖最后由 zhlong 于 2007-7-11 07:59 编辑 ]
回复 支持 2 反对 0

使用道具 举报

发表于 2007-4-27 22:07 | 显示全部楼层

程序结果

ms.jpg
发表于 2007-4-28 13:42 | 显示全部楼层

回复 #38 破凰 的帖子

强啊!应该是这样的。
发表于 2007-4-29 12:54 | 显示全部楼层

回复 #38 破凰 的帖子

看了你的程序,我有几点不明白,首先y=[y1;y2;zeros(size(y1))]; %IMF集这句代表的含义是什么?你好像没有作EMD,哪里会有IMF,还有就是E=flipud(E); ,  这句的作用是什么?

[ 本帖最后由 zhlong 于 2007-7-27 17:10 编辑 ]
发表于 2007-4-29 15:13 | 显示全部楼层
请问toimage函数在哪里找到啊,MATLAB里没有这个函数。。
发表于 2007-4-29 16:38 | 显示全部楼层

回复 #41 zhangnan3509 的帖子

一个正弦函数本身就是一个IMF,所以y=[y1;y2;zeros(size(y1))]就是一个IMF集(当然假定了残余函数为0)。
ps:调用了toimage后得到的结果才是真正的Hilbert谱!
flipud是一个使矩阵上下翻转的函数。在Grilling 提供的程序toimage中,频率是从上往下递增,而通常在时频图中频率应是从下往上递增,所以使用flipud将矩阵翻转后,更便于我们阅读时频图。对于边际谱来说,如果不对E翻转,边际谱图中的频率将是从从右往左递增的。

[ 本帖最后由 破凰 于 2007-4-29 16:49 编辑 ]

评分

1

查看全部评分

发表于 2007-4-29 16:41 | 显示全部楼层
回复 支持 1 反对 0

使用道具 举报

发表于 2007-4-29 17:55 | 显示全部楼层

回复 #44 破凰 的帖子

感谢破凰,这条我收下了!这也是目前我在论坛上看见比较完整的边际谱程序。这对于做HHT边际谱又多出了新的办法。非常感谢!这也是此主题帖子一个比较好的探讨结果!

[ 本帖最后由 zhangnan3509 于 2007-4-29 17:57 编辑 ]
发表于 2007-4-29 22:29 | 显示全部楼层

回复 #47 michaeldy 的帖子

数据我师兄已经全拿走了,我没留备份。因为这些数据对我做的东西没有意义。我得重新做实验,现在还没开始呢:@)
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

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

GMT+8, 2024-3-29 01:50 , Processed in 0.106279 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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