搜索
1、Hilbert边际谱我觉得既然已经做出EMD了，也就是得到了IMF。这个时候就是做hilbert幅值谱，然后对它积分就可以了。程序不是很难搞到吧！我是用hspec画谱图的，自己又在后面添加了求边际谱的代码，但感觉有问题fork=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];%信号画出来的图很粗糙，更不用说对实际信号分析了，所以大家看看如何来修正？黄文章中边际谱对实际信号分析是很好的一条曲线我用hhspectrum算了一下谱图，同时求了一下边际谱，边际谱程序基本想法同form。结果也不太好，20HZ处还行，40HZ就有些问题了，见附图你自己再用这个试试我没有用rilling的hhspectrumnspab:functionh1=nspab(data,nyy,minw,maxw,dt)%ThefunctionNSPABgeneratesasmoothedHHTspectrumofdata(n,k)%intime-frequencyspace,where%nspecifiesthelengthoftimeseries,and%kisthenumberofIMFcomponents.%Thefrequency-axisrangeisprefixed.%Negativefrequencysignisreversed.%%MATLABLibraryfunctionHILBERTisusedtocalculatetheHilberttransform.%%Example,[h,xs,w]=nspab(lod78_p',200,0,0.12,1,3224).%%FunctionsCONTOURorIMGcanbeusedtoviewthespectrum,%forexamplecontour(xs,w,h)orimg(xs,w,h).%%Callingsequence-%[h,xs,w]=nspab(data,nyy,minw,maxw,t0,t1)%%Input-%data-2-Dmatrixdata(n,k)ofIMFcomponents%nyy-thefrequencyresolution%minw-theminimumfrequency%maxw-themaximumfrequency%t0-thestarttime%t1-theendtime%Output-%h-2-DmatrixoftheHHTspectrum,where%the1stdimensionspecifiesthenumberoffrequencies,%the2nddimensionspecifiesthenumberoftimevalues%xs-vectorthatspecifiesthetime-axisvalues%w-vectorthatspecifiesthefrequency-axisvalues%Z.Shen(JHU)July2,1995Initial%-----Getdimensions(numberoftimepointsandcomponents)[npt,knb]=size(data);%-----Gettimeinterval%-----ApplyHilbertTransformdata=hilbert(data);a=abs(data);omg=abs(diff(unwrap(angle(data))))/(2*pi*dt);%-----Smoothamplitudeandfrequencyfiltr=fir1(8,.1);fori=1:knba(:,i)=filtfilt(filtr,1,a(:,i));omg(:,i)=filtfilt(filtr,1,omg(:,i));end%-----Limitfrequencyandamplitudefori=1:knbfori1=1:npt-1ifomg(i1,i)=maxw,omg(i1,i)=maxw;a(i1,i)=0;elseifomg(i1,i)=minw,omg(i1,i)=minw;a(i1,i)=0;elseendendendclearfiltrdata%va=var(omg(200:1200))%-----Getlocalfrequencydw=maxw-minw;wmx=maxw;wmn=minw;%-----Constructtheplotingmatrixclearp;h1=zeros(npt-1,nyy+1);p=round(nyy*(omg-wmn)/dw)+1;forj1=1:npt-1fori1=1:knbii1=p(j1,i1);h1(j1,ii1)=h1(j1,ii1)+a(j1,i1);endend%-----Do3-pointto1-pointaveraging[nx,ny]=size(h1);%n1=fix(nx/3);%h=zeros(n1,ny);%fori1=1:n1%h(i1,:)=(h1(3*i1,:)+h1(3*i1-1,:)+h1(3*i1-2,:));%end%clearh1;%-----Do3-pointssmoothinginx-directionfltr=1./3*ones(3,1);forj1=1:nyh1(:,j1)=filtfilt(fltr,1,h1(:,j1));endclearfltr;%-----Definetheresults%w=linspace(wmn,wmx,ny-1)';%xs=linspace(t0,t1,nx)';h1=flipud(rot90(h1));h1=h1(1:ny-1,:);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);fork=1:size(E,1)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后得到的结果才是真正的Hilbert谱！flipud是一个使矩阵上下翻转的函数。在Grilling提供的程序toimage中，频率是从上往下递增，而通常在时频图中频率应是从下往上递增，所以使用flipud将矩阵翻转后，更便于我们阅读时频图。对于边际谱来说，如果不对E翻转，边际谱图中的频率将是从从右往左递增的。2、Hilbert边际谱和FT变换后的幅频谱这得出的Hilbert边际谱和FT变换后的幅频谱为什么会有这么大的区别呢，到底哪个幅值才是真正的实际幅值呢？有参考价值吗？程序如下：loadshujufs=5120;N=4096;a1=a(1:N,1);a2=abs(fft(a1))*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));fork=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边际谱')答：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*pi;b1=linspace(0,0,1000);fori=250:350b1(i)=0.2*sin(t2(i));endfori=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);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但是毕竟是有限次平均逼近存在的残余也很正常我仔细分析了下图中第一项是信号本身第二项是残余噪声（随机噪声一般是高频）第三项是仿真信号中的b1……不过这个分解结果于我来看还是不理想待分解信号中断正弦信号最后没在一个IMF分量中在没有噪音干扰的情况下，采用EMD得到的分解图可以说是正确的。我看到一些文献，总是认为EMD分解就是从高频向低频的分解，实际而言，这种理解是对EMD分解得到单分量IMF的一种误解。如果采用瞬时频率去描述单分量信号，我认为是指信号在某一时刻仅存在一个频率成分。所以对于原信号而言，[0250]、[350750]与[8501000]区间内的信号只含有10Hz的低频正弦信号，所以在EMD分解过程中，这些信号成分应该