我是用hspec画谱图的，自己又在后面添加了求边际谱的代码，但感觉有问题
for k=1:size(E)
bjp(k)=sum(E(k,:))*1/fs;%fs为采样频率;
end
figure
plot(bjp);
xlabel('频率 / Hz');
ylabel('幅值');
...
前面这是我做的，你是在hhspectrum后面加的吗？ 我不是这么做的！我用的是nspab.m和mspc.m这两个做的。你也应该有吧
[ 本帖最后由 zhangnan3509 于 2007-4-8 21:50 编辑 ] 原帖由 zhangnan3509 于 2007-4-8 21:45 发表
我不是这么做的！我用的是nspab.m和mspc.m这两个做的。你也应该有吧
你两个函数是自己编的吗？我下载的程序包没有这两个函数啊？
我第二次作的是在hhspectrum后面添加的，感觉按描述的应该是在这个函数基础上来作，那你做的图和我的比较如何？
[ 本帖最后由 zhlong 于 2007-6-1 20:12 编辑 ]
回复 #34 form 的帖子
你自己再用这个试试 我没有用rilling的hhspectrumnspab:
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, = 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-
% = 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)
= 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
=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,:); 请问HHT边际谱和FFT谱有什么区别，差别在哪里？:loveliness: 看了上面的程序，真是有点思路！我记得小波分析功率谱图就是这样画的 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=; %IMF集
%%%%%%%%%%%%%求边际谱
=hhspectrum(y);
=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 编辑 ]
程序结果
回复 #38 破凰 的帖子
强啊！应该是这样的。回复 #38 破凰 的帖子
看了你的程序，我有几点不明白，首先y=; %IMF集这句代表的含义是什么？你好像没有作EMD，哪里会有IMF，还有就是E=flipud(E); ，这句的作用是什么？[ 本帖最后由 zhlong 于 2007-7-27 17:10 编辑 ] 请问toimage函数在哪里找到啊，MATLAB里没有这个函数。。
回复 #41 zhangnan3509 的帖子
一个正弦函数本身就是一个IMF，所以y=就是一个IMF集（当然假定了残余函数为0）。ps：调用了toimage后得到的结果才是真正的Hilbert谱！
flipud是一个使矩阵上下翻转的函数。在Grilling 提供的程序toimage中，频率是从上往下递增，而通常在时频图中频率应是从下往上递增，所以使用flipud将矩阵翻转后，更便于我们阅读时频图。对于边际谱来说，如果不对E翻转，边际谱图中的频率将是从从右往左递增的。
[ 本帖最后由 破凰 于 2007-4-29 16:49 编辑 ] 原帖由 michaeldy 于 2007-4-29 15:13 发表 http://forum.vibunion.com/forum/images/common/back.gif
请问toimage函数在哪里找到啊，MATLAB里没有这个函数。。
这是法国的Grilling提供的开源程序。
下载地址：http://perso.ens-lyon.fr/patrick.flandrin/emd.html
回复 #44 破凰 的帖子
感谢破凰，这条我收下了！这也是目前我在论坛上看见比较完整的边际谱程序。这对于做HHT边际谱又多出了新的办法。非常感谢！这也是此主题帖子一个比较好的探讨结果！[ 本帖最后由 zhangnan3509 于 2007-4-29 17:57 编辑 ]