EMD.rar_EMD 程序_decomposition_emd matlab程序_emd程序_经验

function [Imf,nbits]=fun_emd(Sig) if nargin<1 error('At least one input is needed!'); end %Sig=Sig'; c=Sig; SigLen=length(Sig); t=1:SigLen; c=Sig; Imf=[]; nbits=[]; MaxLoop=1000; [IndMax,IndMin]=ExtrPoint(c);%提取信号的极值点 [Index_zer]=ZeroPoint(c);%提取信号的零点 NExtr=length(IndMin)+length(IndMax); NZero=length(Index_zer); Bool=1; NEtd=2; IMFNum=0; while NExtr>2&&Bool>0 h=c; nLoop=1; SD=1; [IndMax,IndMin]=ExtrPoint(c);%提取信号的极值点%%%%% [Index_zer]=ZeroPoint(c);%提取信号的零点%%%%%%%%%%%%%%% NExtr=length(IndMin)+length(IndMax); NZero=length(Index_zer); %Bool=1; %NEtd=2; %IMFNum=0; while (SD>0.3||(abs(NZero-NExtr)>1))&&(NExtr>2)&&nLoop<MaxLoop %镜像拓展序列 LMaxEtd=fliplr(IndMax(1:min(end,NEtd))); LMinEtd=fliplr(IndMin(1:min(end,NEtd))); hlMaxEtd=h(LMaxEtd); hlMinEtd=h(LMinEtd); LMaxEtd=-LMaxEtd+1; LMinEtd=-LMinEtd+1; RMaxEtd=fliplr(IndMax(max(end-NEtd+1,1):end)); RMinEtd=fliplr(IndMin(max(end-NEtd+1,1):end)); hrMaxEtd=h(RMaxEtd); hrMinEtd=h(RMinEtd); RMaxEtd=2*SigLen-RMaxEtd; RMinEtd=2*SigLen-RMinEtd; %计算信号的上包络线 Max_Env=interp1([LMaxEtd,t(IndMax),RMaxEtd],[hlMaxEtd,h(IndMax),hrMaxEtd],t,'spline'); %计算信号的下包络线 Min_Env=interp1([LMinEtd,t(IndMin),RMinEtd],[hlMinEtd,h(IndMin),hrMinEtd],t,'spline'); Mean_Env=(Max_Env+Min_Env)/2;%计算包络线平均值 Preh=h; h=h-Mean_Env; alarm=1.0e-8; SD=sum(((Preh-h).^2)./(Preh.^2+alarm)); [IndMax,IndMin]=ExtrPoint(h); Index_zero=ZeroPoint(h);%%%%%%%%%%%%%%%%%%%%% NExtr=length(IndMin)+length(IndMax); NZero=length(Index_zero);%%%%%%%%%%%%%%%%%%%%%% nLoop=nLoop+1; end %提取到一个本征模分量 IMFNum=IMFNum+1; disp([num2str(IMFNum),'IMFs have been obtained.']); Imf=[Imf;h]; c=c-h; [IndMax,IndMin]=ExtrPoint(c);%提取信号的极值点 Index_zer=ZeroPoint(c);%提取信号的零点 NExtr=length(IndMin)+length(IndMax); NZero=length(Index_zer); nbits=[nbits;nLoop-1]; Bool=1; if (range(c)/range(Sig))<2e-4; Bool=-1; end end Imf=[Imf;c]; [n,m]=size(Imf); xCor=[]; end %提取信号的极大值和极小值 function [Pos_max,Pos_min]=ExtrPoint(Sig) if nargin==0 error('At least one input is require.'); end SigLen=length(Sig); DSig=diff(Sig); DSig1=DSig(1:end-1); DSig2=DSig(2:end); Pos_min=find(DSig1.*DSig2<0&DSig1<0)+1; Pos_max=find(DSig1.*DSig2<0&DSig1>0)+1; Pos_max=sort(Pos_max); Pos_min=sort(Pos_min); end %提取信号的零点 function PZero=ZeroPoint(Sig) if nargin==0 error('At least one inputis requere'); end SigLen=length(Sig); s1=Sig(1:SigLen-1); s2=Sig(2:SigLen); PZero=find(s1.*s2<0); IndZero=find(Sig==0, 1); if ~isempty(IndZero) zero=find(Sig==0); DZero=ifft([0,zero,0]);%%%%%%%%%%%%dfft LZero=find(DZero==1); RZero=find(DZero==-1); IndZero=round((LZero+RZero)/2); PZero=sort([PZero,IndZero]); end end