边际谱m程序问题

babybear713

用了两种方法求边际谱~为何结果相差这么大呢？毕业在即，请高手指教，不甚感激
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的时频图
[imf,ort,nbits]=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 integer  values-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([1 1000 1 500 0 0.05])
% 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([0 772 0 300]); %%%%%%%%%%%%%%%%%%%%%%
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.接上
%%% 求边际谱
[A,fa,tt]=hhspectrum(imf(1:end-1,:));  % A--瞬时幅值，fa--瞬时频率，tt--
[E,tt1,ff]=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('幅值');

babybear713

另，前一种方法 得到的横坐标不知是否要除以1000才是其真实的频率，若是，那范围只有0~0.5hz，未免也太小了一点吧
SOS:@)

babybear713

function [S,freq] = 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是不是一个随便设定的值，取决于你想得到的精度。

cboboc

原帖由 babybear713 于 2010-4-8 21:49 发表

                               
登录/注册后可看大图

另，前一种方法 得到的横坐标不知是否要除以1000才是其真实的频率，若是，那范围只有0~0.5hz，未免也太小了一点吧
SOS:@)

不是除以1000吧，你的采样频率不是50吗？应该除以50吧

babybear713

这个f 的设定是在画hht频谱图的时候用到的，因为根据采样定理是结构的最大频率应该不超过采样频率的一半，所以自认为它的最大频率为25

xiangyu537

请问各位？ EMD后，如果采样频率是fS ,采样点数N  ，EMD后 只对IMF1做边际谱，所得到的频谱图的频率分辨率是fS /2N吗？  跟EMD的分解层数有关系吗？

chunhuajia

对楼主的问题，我也遇到了，希望高手指点！

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？

babybear713

呼唤版主及各位高人呀

杨德昌

其实您就是高人了！

babybear713

我是诚心求教

summerback

你好，一直困惑为什么
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);