1、1、Hilbert 边际谱我觉得既然已经做出 EMD 了，也就是得到了 IMF。这个时候就是做 hilbert 幅值谱，然后对它积分就可以了。程序不是很难搞到吧！我是用 hspec 画谱图的，自己又在后面添加了求边际谱的代码for k=1:size(E)bjp(k)=sum(E(k,:)*1/fs; %fs 为采样频率;endfigureplot(bjp);xlabel(频率 / Hz);ylabel(幅值);比如我用两个正弦信号作仿真fs=1000;t=1/fs:1/fs:1;y1=2*sin(40*pi*t);y2=5*sin(80*pi*t);y=y1,y2; % 信号画出来的图很粗糙，
2、更不用说对实际信号分析了，所以大家看看如何来修正？黄文章中边际谱对实际信号分析是很好的一条曲线我用 hhspectrum 算了一下谱图，同时求了一下边际谱，边际谱程序基本想法同 form。结果也不太好，20HZ 处还行， 40HZ 就有些问题了，见附图你自己再用这个试试 我没有用 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 spac
3、e, 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).% F
4、unctions 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
5、 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 freque
6、ncy-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 Transformdata=hilbert(data);a=abs(data);omg=abs(diff(unwrap(angle(data)/(2*pi*dt);%- Smooth amplitude and frequency filtr=fir1(8,.1); f
7、or i=1:knb a(:,i)=filtfilt(filtr,1,a(:,i);omg(:,i)=filtfilt(filtr,1,omg(:,i);end%- Limit frequency and amplitudefor i=1:knbfor i1=1:npt-1if 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;elseendendendclear filtr data%va=var(omg(200:1200)%- Get local frequency
8、dw=maxw - minw;wmx=maxw;wmn=minw;%- Construct the ploting matrixclear p;h1=zeros(npt-1,nyy+1);p=round(nyy*(omg-wmn)/dw)+1;for j1=1:npt-1for i1=1:knbii1=p(j1,i1); h1(j1,ii1)=h1(j1,ii1)+a(j1,i1);endend%- Do 3-point to 1-point averagingnx,ny=size(h1);%n1=fix(nx/3);%h=zeros(n1,ny);%for i1=1:n1%h(i1,:)=(
9、h1(3*i1,:)+h1(3*i1-1,:)+h1(3*i1-2,:);%end%clear h1;%- Do 3-points smoothing in x-directionfltr=1./3*ones(3,1);for j1=1:nyh1(:,j1)=filtfilt(fltr,1,h1(:,j1);endclear fltr;%- Define the results%w=linspace(wmn,wmx,ny-1);%xs=linspace(t0,t1,nx);h1=flipud(rot90(h1);h1=h1(1:ny-1,:);form 求边际谱时所用程序是没有问题的，用的是矩
10、形积分公式。他所得结果不正确的原因是：输入的应是调用了 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
11、)bjp(k)=sum(E(k,:)*1/fs; endf=(0:N-3)/N*(fs/2);plot(f,bjp);xlabel(频率 / Hz);ylabel(幅值);%（完整答案）问：看了你的程序，我有几点不明白，首先 y=y1;y2;zeros(size(y1); %IMF 集这句代表的含义是什么？你好像没有作 EMD，哪里会有 IMF，还有就是 E=flipud(E); ， 这句的作用是什么？答：一个正弦函数本身就是一个 IMF，所以 y=y1;y2;zeros(size(y1)就是一个 IMF 集（当然假定了残余函数为 0） 。ps：调用了 toimage 后得到的结果才是真正的
12、Hilbert 谱！flipud 是一个使矩阵上下翻转的函数。在 Grilling 提供的程序 toimage 中，频率是从上往下递增，而通常在时频图中频率应是从下往上递增，所以使用 flipud 将矩阵翻转后，更便于我们阅读时频图。对于边际谱来说，如果不对 E 翻转，边际谱图中的频率将是从从右往左递增的。2、Hilbert 边际谱和 FT 变换后的幅频谱这得出的 Hilbert 边际谱和 FT 变换后的幅频谱为什么会有这么大的区别呢，到底哪个幅值才是真正的实际幅值呢？有参考价值吗？ 程序如下：load shujufs=5120;N=4096;a1=a(1:N,1);a2=abs(fft(a1
13、)*2/N;f=fs*(0:N/2-1)/N;n=length(f);subplot(211)plot(f,a2(1:n)xlabel(频率 / Hz);ylabel(幅值);title(FT 后的幅频图)imf=emd(a1);A,fa,tt=hhspectrum(imf);E,tt1=toimage(A,fa,tt,length(tt);for k=1:size(E,1)bjp(k)=sum(E(k,:)*1/fs;endf=(0:N-3)/N*(fs/2);subplot(212)plot(f,bjp);xlabel(频率 / Hz);ylabel(幅值);title(Hilbert 边
14、际谱)答：EMD 分解的 IMF 能量和原信号其能量是不相等的。所以边际谱能量不能和 FFT 的能量相比。边际谱能量只能说明某个信号存在，能量相对于其他的大小。边际谱是对 IMF 取包络线。因此它得到的谱能量要大于被取包络的信号能量。FFT 谱和原信号的能量是相等的。所以从能量的大小讲，应该是边际谱能量大于 FFT 的能量。如果是一个谐波取边际谱和 FFT 的话，应该是频率对应的能量边际谱大于 FFT。3、EEMD 的一些问题刚刚接触 EEMD，从台湾中央大学上下载了程序，然后编写程序如下：clear;clct=1:1000;t1=t/100*2*pi;a1=sin(t1);t2=t/10*2
15、*pi;b1=linspace(0,0,1000);for i=250:350b1(i)=0.2*sin(t2(i);endfor i=750:850b1(i)=0.2*sin(t2(i);endx=a1+b1; %x 是原信号subplot(311);plot(t,a1);subplot(312);plot(t,b1);subplot(313);plot(t,x);plot(x)imf=emd(x);emd_visu(x,t,imf)%eemdimf_eemd=eemd(x,0.1,100);figuresubplot(511);plot(imf_eemd(:,1)subplot(512);
16、plot(imf_eemd(:,2)subplot(513);plot(imf_eemd(:,3)subplot(514);plot(imf_eemd(:,4)subplot(515);plot(imf_eemd(:,5)figuresubplot(511);plot(imf_eemd(:,6)subplot(512);plot(imf_eemd(:,7)subplot(513);plot(imf_eemd(:,8)subplot(514);plot(imf_eemd(:,9)subplot(515);plot(imf_eemd(:,10)原信号：EMD:EEMD:感觉 EEMD 的第二项为什么会是频率这么高？和 EMD 差距这么大？是不是我程序有问题？哪位高手请赐教，不胜感激！答：我认为第二项应该是随机噪声造成的虽然集总平均理论上为 0但是毕竟是有限次平均逼近存在的残余也很正常
