Hilbert边际谱程序（matlab）

form

原帖由 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后面加的吗？

zhangnan3509

我不是这么做的！我用的是nspab.m和mspc.m这两个做的。你也应该有吧

[ 本帖最后由 zhangnan3509 于 2007-4-8 21:50 编辑 ]

form

原帖由 zhangnan3509 于 2007-4-8 21:45 发表
我不是这么做的！我用的是nspab.m和mspc.m这两个做的。你也应该有吧

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

[ 本帖最后由 zhlong 于 2007-6-1 20:12 编辑 ]

zhangnan3509

你自己再用这个试试 我没有用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,:);

bingbing1110

请问HHT边际谱和FFT谱有什么区别，差别在哪里？:loveliness:

wyong0210

看了上面的程序，真是有点思路！我记得小波分析功率谱图就是这样画的

破凰

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 编辑 ]

破凰

wyong0210

强啊！应该是这样的。

zhangnan3509

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

[ 本帖最后由 zhlong 于 2007-7-27 17:10 编辑 ]

michaeldy

请问toimage函数在哪里找到啊，MATLAB里没有这个函数。。

破凰

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

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

破凰

原帖由 michaeldy 于 2007-4-29 15:13 发表

                               
登录/注册后可看大图

请问toimage函数在哪里找到啊，MATLAB里没有这个函数。。

这是法国的Grilling提供的开源程序。
下载地址：http://perso.ens-lyon.fr/patrick.flandrin/emd.html

zhangnan3509

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

[ 本帖最后由 zhangnan3509 于 2007-4-29 17:57 编辑 ]

zhangnan3509

数据我师兄已经全拿走了，我没留备份。因为这些数据对我做的东西没有意义。我得重新做实验，现在还没开始呢:@)