emd算法-matlab编写-亲测可用

%ThisisanEMD/EEMDprogram%%functionallmode=eemd(Y,Nstd,NE)%%INPUT:%Y:Inputteddata;%Nstd:ratioofthestandarddeviationoftheaddednoiseandthatofY;%NE:EnsemblenumberfortheEEMD%OUTPUT:%AmatrixofN*(m+1)matrix,whereNisthelengthoftheinput%dataY,andm=fix(log2(N))-1.Column1istheoriginaldata,columns2,3,...%maretheIMFsfromhightolowfrequency,andcomlumn(m+1)isthe%residual(overalltrend).%%NOTE:%ItshouldbenotedthatwhenNstdissettozeroandNEissetto1,the%programdegeneratestoaEMDprogram.%%ReferencescanbefoundintheReferencesection.%%ThecodeispreparedbyZhaohuaWu.Forquestions,pleasereadtheQ&[e-mail]%functionallmode=eemd(Y,Nstd,NE)xsize=length(Y);dd=1:1:xsize;Ystd=std(Y);Y=Y/Ystd;TNM=fix(log2(xsize))-1;TNM2=TNM+2;forkk=1:1:TNM2,forii=1:1:xsize,allmode(ii,kk)=0.0;endendforiii=1:1:NE,fori=1:xsize,temp=randn(1,1)*Nstd;X1(i)=Y(i)+temp;endforjj=1:1:xsize,mode(jj,1)=Y(jj);endxorigin=X1;xend=xorigin;nmode=1;whilenmode=TNM,xstart=xend;iter=1;whileiter=10,[spmax,spmin,flag]=extrema(xstart);upper=spline(spmax(:,1),spmax(:,2),dd);lower=spline(spmin(:,1),spmin(:,2),dd);mean_ul=(upper+lower)/2;xstart=xstart-mean_ul;iter=iter+1;endxend=xend-xstart;nmode=nmode+1;forjj=1:1:xsize,mode(jj,nmode)=xstart(jj);endendforjj=1:1:xsize,mode(jj,nmode+1)=xend(jj);endallmode=allmode+mode;endallmode=allmode/NE;allmode=allmode*Ystd;%Thisisautilityprogramforsignificancetest.%%function[spmax,spmin,flag]=extrema(in_data)%%INPUT:%in_data:Inputteddata,atimeseriestobesifted;%OUTPUT:%spmax:Thelocations(col1)ofthemaximaanditscorresponding%values(col2)%spmin:Thelocations(col1)oftheminimaanditscorresponding%values(col2)%%ReferencescanbefoundintheReferencesection.%%ThecodeispreparedbyZhaohuaWu.Forquestions,pleasereadtheQ&[e-mail]%function[spmax,spmin,flag]=extrema(in_data)flag=1;dsize=length(in_data);spmax(1,1)=1;spmax(1,2)=in_data(1);jj=2;kk=2;whilejjdsize,if(in_data(jj-1)=in_data(jj)&in_data(jj)=in_data(jj+1))spmax(kk,1)=jj;spmax(kk,2)=in_data(jj);kk=kk+1;endjj=jj+1;endspmax(kk,1)=dsize;spmax(kk,2)=in_data(dsize);ifkk=4slope1=(spmax(2,2)-spmax(3,2))/(spmax(2,1)-spmax(3,1));tmp1=slope1*(spmax(1,1)-spmax(2,1))+spmax(2,2);iftmp1spmax(1,2)spmax(1,2)=tmp1;endslope2=(spmax(kk-1,2)-spmax(kk-2,2))/(spmax(kk-1,1)-spmax(kk-2,1));tmp2=slope2*(spmax(kk,1)-spmax(kk-1,1))+spmax(kk-1,2);iftmp2spmax(kk,2)spmax(kk,2)=tmp2;endelseflag=-1;endmsize=size(in_data);dsize=max(msize);xsize=dsize/3;xsize2=2*xsize;spmin(1,1)=1;spmin(1,2)=in_data(1);jj=2;kk=2;whilejjdsize,if(in_data(jj-1)=in_data(jj)&in_data(jj)=in_data(jj+1))spmin(kk,1)=jj;spmin(kk,2)=in_data(jj);kk=kk+1;endjj=jj+1;endspmin(kk,1)=dsize;spmin(kk,2)=in_data(dsize);ifkk=4slope1=(spmin(2,2)-spmin(3,2))/(spmin(2,1)-spmin(3,1));tmp1=slope1*(spmin(1,1)-spmin(2,1))+spmin(2,2);iftmp1spmin(1,2)spmin(1,2)=tmp1;endslope2=(spmin(kk-1,2)-spmin(kk-2,2))/(spmin(kk-1,1)-spmin(kk-2,1));tmp2=slope2*(spmin(kk,1)-spmin(kk-1,1))+spmin(kk-1,2);iftmp2spmin(kk,2)spmin(kk,2)=tmp2;endelseflag=-1;endflag=1;