边际谱m程序问题
用了两种方法求边际谱~为何结果相差这么大呢？毕业在即，请高手指教，不甚感激1.
clear all
clc
load('data.txt');
x=data(1:1000,2);
t=1:1000;
fs=50;
figure(10);
plot(1/fs*t,x);
xlabel('t/s');
ylabel('m/s^2');
%对信号作EMD分解，求出各IMF分量图和重构图及各IMF的时频图
=emd(x);
%作出HHT归一化能量谱图（三维图），边际谱图和瞬时能量图
N=1000;
L=size(imf,1);% Number of components in the decomposition
% loop for on each component
S=[];
clear x z m p freq
x=imf'; % now each column is a compenent
z=hilbert(x); % analytic signal
m=abs(z); % module of z
p=angle(z); % phase of z
for i=1:L
freq(:,i)=instfreq(z(:,i)); % instantaneous frequency
ceilfreq(:,i)=ceil(freq(:,i)*N);% to have integervalues-also do a smoothing given the number
for j=1:length(x)-2
S(ceilfreq(j,i),j+1)=m(j+1,i);
end
end
eps=0.00001; % to avoid zero values before the log
E=20*log(S.^2+eps);% Hilbert energ spectrum
% plot S
figure(90);
t=1:length(x);% time sample
f=t/length(x)*25; % frequency %%%%%%%%%
imagesc(t,f,E);
contour(t(1:999),1:50,E);%画出等高线图
surfl(1:999,1:500,S.^2); shading interp; colormap(gray);axis()
% mesh(1:999,1:50,S.^2);
% mesh(1:999,1:50,S); %画出三维图
set(gca,'YDir','normal');
xlabel('时间采样点');
ylabel('Frequency/Hz'); %%%%%%%%%%%%%%%%%%%%
% ylabel('Normalized Frequency');
% axis(); %%%%%%%%%%%%%%%%%%%%%%
S1=sum(S.^2,1);figure(100);plot((1:999)*0.02,S1); %瞬时能量
figure(110);
S2=sum(S.^2,2);plot(1:500,S2,'r') %边际能量谱
figure(120);
S3=sum(S,2);plot(1:500,S3,'g') %边际谱
2.接上
%%% 求边际谱
=hhspectrum(imf(1:end-1,:));% A--瞬时幅值，fa--瞬时频率，tt--
=toimage(A,fa,tt,length(tt));%% ff--中心频率
% disp_hhs(E);colorbar;
% colormap(flipud(gray)) %显示灰度图
for k=1:size(E,1)
bjp(k)=sum(E(k,:))*1/fs;
end
f2=(0:N-3)/N*(fs/2); % fs/(2*N)是时频图上y方向的分辨率
figure(200)
plot(f2,bjp);%%% 作边际谱图
xlabel('频率 / Hz');
ylabel('幅值');
另，前一种方法 得到的横坐标不知是否要除以1000才是其真实的频率，若是，那范围只有0~0.5hz，未免也太小了一点吧
SOS:@)
回复 楼主 babybear713 的帖子
function = HHTspe(imf,N,fs);%%N--------分辨率
%fs-------采样频率
%s--------时间-频率-幅值矩阵
if(nargin<3)
fs=1;
end
L = size(imf,1);
S = [];
clear x z m p freq
x = imf';
z = hilbtmf(x); % 用改进的Hilbert变换，效果好一些。
m = abs(z);
p = angle(z);
for i = 1:L
freq(:,i) = instfreq(z(:,i))*fs; %乘以了采样频率
ceilfreq(:,i) = ceil(freq(:,i)*N);
for j = 1:length(x)-2
S(ceilfreq(j,i),j+1) = m(j+1,i);
end
end
%plot S
figure;
t=(1:length(x))/fs;
f=(1:size(S,1))/N; %
imagesc(t,f,S);
colorbar;
set(gca,'YDir','normal');
xlabel('Time t/s');
ylabel(' Frequency/HZ');
return
再请教，这个分辨率N是如何确定的？？？ （1）第一个图不用除以1000，应该除以25，在您的程序中有一句，f=t/length(x)*25;，我不明白为什么您要乘以25。
（2）这个N的意思，我也一直没有明白，但是您的第一个程序给了我一些启示，我个人理解有两种：（1）N为采样频率 （2）N为采样点的个数也就是length（t），我用试验的方法做了一下，N为采样点数的时候效果好，但是计算时间较长。所以N是不是一个随便设定的值，取决于你想得到的精度。 原帖由 babybear713 于 2010-4-8 21:49 发表 http://www.chinavib.com/forum/images/common/back.gif
另，前一种方法 得到的横坐标不知是否要除以1000才是其真实的频率，若是，那范围只有0~0.5hz，未免也太小了一点吧
SOS:@)
不是除以1000吧，你的采样频率不是50吗？应该除以50吧
回复 地板 杨德昌 的帖子
这个f 的设定是在画hht频谱图的时候用到的，因为根据采样定理是结构的最大频率应该不超过采样频率的一半，所以自认为它的最大频率为25 请问各位？ EMD后，如果采样频率是fS ,采样点数N，EMD后 只对IMF1做边际谱，所得到的频谱图的频率分辨率是fS /2N吗？跟EMD的分解层数有关系吗？ 对楼主的问题，我也遇到了，希望高手指点！回复 楼主 babybear713 的帖子
发现了一点问题，横坐标应该是（1：500）*fs/N,这样好些就差不多了还有两个问题就是
1. 能量谱的话用
for k=1:size(E,1)
bjp(k)=sum(E(k,:).^2)*1/fs;
可以么？又或者
for k=1:size(E,1)
bjp(k)=sum(E(k,:).*E(k,:))*1/fs*1/tspan;
2什么时候后面要加上乘以1/tspan？
回复 9楼 babybear713 的帖子
呼唤版主及各位高人呀回复 10楼 babybear713 的帖子
其实您就是高人了！回复 11楼 杨德昌 的帖子
我是诚心求教 你好，一直困惑为什么f2=(0:N-3)/N*(fs/2); % fs/(2*N)是时频图上y方向的分辨率
中是0：N-3？这代表着什么呢？ f2=(0:N-3)/N*(fs/2); 应该改为(0:399)/400*(fs/2);
页:
[1]