经验模态分解+样本熵matlab程序代码

%当X是一个真正的矢量时，计算X的经验模态分解， %每一行包含一个max的经验模态？最后一个是残余向量，默认准则在文献2提出 %本函数用于对原始信号varargin进行分解，imf为分解出的本征模态向量 %函数返回值imf为本征模态向量,ort=一个数，nbits为1*4的向量？ function [imf,ort,nbits] = emd(varargin) [x,t,sd,sd2,tol,MODE_COMPLEX,ndirs,display_sifting,sdt,sd2t,r,imf,k,nbit,NbIt,MAXITERATIONS,FIXE,FIXE_H,MAXMODES,INTERP,mask] = init(varargin{:}); if display_sifting fig_h = figure; end %主循环:至少需要3极值才能继续 while ~stop_EMD(r,MODE_COMPLEX,ndirs) && (k < MAXMODES+1 || MAXMODES == 0) && ~any(mask) % 当前模式 m = r; % 在前一个迭代模式 mp = m; %计算平均值和停止准则 if FIXE [stop_sift,moyenne] = stop_sifting_fixe(t,m,INTERP,MODE_COMPLEX,ndirs); elseif FIXE_H stop_count = 0; [stop_sift,moyenne] = stop_sifting_fixe_h(t,m,INTERP,stop_count,FIXE_H,MODE_COMPLEX,ndirs); else [stop_sift,moyenne] = stop_sifting(m,t,sd,sd2,tol,INTERP,MODE_COMPLEX,ndirs); end % in case the current mode is so small that machine precision can cause %如果当前模式太小,机器精度的原因 % spurious extrema to appear %伪极值出现 if (max(abs(m))) < (1e-16)*(max(abs(x))) if ~stop_sift warning('emd:warning','forced stop of EMD : too small amplitude') else disp('forced stop of EMD : too small amplitude') end break end % 筛选循环 while ~stop_sift && nbit<MAXITERATIONS if(~MODE_COMPLEX && nbit>MAXITERATIONS/5 && mod(nbit,floor(MAXITERATIONS/10))==0 && ~FIXE && nbit > 100) disp(['mode ',int2str(k),', iteration ',int2str(nbit)]) if exist('s','var') disp(['stop parameter mean value : ',num2str(s)]) end [im,iM] = extr(m); disp([int2str(sum(m(im) > 0)),' minima > 0; ',int2str(sum(m(iM) < 0)),' maxima < 0.']) end %筛选 m = m - moyenne; %计算平均值和停止准则 if FIXE [stop_sift,moyenne] = stop_sifting_fixe(t,m,INTERP,MODE_COMPLEX,ndirs); elseif FIXE_H [stop_sift,moyenne,stop_count] = stop_sifting_fixe_h(t,m,INTERP,stop_count,FIXE_H,MODE_COMPLEX,ndirs); else [stop_sift,moyenne,s] = stop_sifting(m,t,sd,sd2,tol,INTERP,MODE_COMPLEX,ndirs); end % display if display_sifting && ~MODE_COMPLEX NBSYM = 2; [indmin,indmax] = extr(mp); [tmin,tmax,mmin,mmax] = boundary_conditions(indmin,indmax,t,mp,mp,NBSYM); envminp = interp1(tmin,mmin,t,INTERP); envmaxp = interp1(tmax,mmax,t,INTERP); envmoyp = (envminp+envmaxp)/2; if FIXE || FIXE_H display_emd_fixe(t,m,mp,r,envminp,envmaxp,envmoyp,nbit,k,display_sifting) else sxp=2*(abs(envmoyp))./(abs(envmaxp-envminp)); sp = mean(sxp); display_emd(t,m,mp,r,envminp,envmaxp,envmoyp,s,sp,sxp,sdt,sd2t,nbit,k,display_sifting,stop_sift) end end mp = m; nbit=nbit+1; NbIt=NbIt+1; if(nbit==(MAXITERATIONS-1) && ~FIXE && nbit > 100) if exist('s','var') warning('emd:warning',['forced stop of sifting : too many iterations... mode ',int2str(k),'. stop parameter mean value : ',num2str(s)]) else warning('emd:warning',['forced stop of sifting : too many iterations... mode ',int2str(k),'.']) end end end % sifting loop imf(k,:) = m; if display_sifting disp(['mode ',int2str(k),' stored']) end nbits(k) = nbit; k = k+1; r = r - m; nbit=0; end %main loop if any(r) && ~any(mask) imf(k,:) = r; end ort = io(x,imf); if display_sifting close end end %--------------------------------------------------------------------------------------------------- % 测试若有足够的极致则继续分解 function stop = stop_EMD(r,MODE_COMPLEX,ndirs) if MODE_COMPLEX for k = 1:ndirs phi = (k-1)*pi/ndirs; [indmin,indmax] = extr(real(exp(i*phi)*r)); ner(k) = length(indmin) + length(indmax); end stop = any(ner < 3); else [indmin,indmax] = extr(r); ner = length(indmin) + length(indmax); stop = ner < 3; end end %--------------------------------------------------------------------------------------------------- % 计算包络线们和振幅估计模式的平均值 %INIERP为插值函数，MODE_COMPLEX复杂模式，ndirs无本金交割利率互换？ function [envmoy,nem,nzm,amp] = mean_and_amplitude(m,t,INTERP,MODE_COMPLEX,ndirs) NBSYM = 2; if MODE_COMPLEX switch MODE_COMPLEX case 1 for k = 1:ndirs phi = (k-1)*pi/ndirs; y = real(exp(-i*phi)*m); [indmin,indmax,indzer] = extr(y); nem(k) = length(indmin)+length(indmax); nzm(k) = length(indzer); [tmin,tmax,zmin,zmax] = boundary_conditions(indmin,indmax,t,y,m,NBSYM); envmin(k,:) = interp1(tmin,zmin,t,INTERP); envmax(k,:) = interp1(tmax,zmax,t,INTERP); end envmoy = mean((envmin+envmax)/2,1); if nargout > 3 amp = mean(abs(envmax-envmin),1)/2; end case 2 for k = 1:ndirs phi = (k-1)*pi/ndirs; y = real(exp(-i*phi)*m); [indmin,indmax,indzer] = extr(y); nem(k) = length(indmin)+length(indmax); nzm(k) = length(indzer); [tmin,tmax,zmin,zmax] = boundary_conditions(indmin,indmax,t,y,y,NBSYM); envmin(k,:) = exp(i*phi)*interp1(tmin,zmin,t,INTERP); envmax(k,:) = exp(i*phi)*interp1(tmax,zmax,t,INTERP); end envmoy = mean((envmin+envmax),1); if nargout > 3 amp = mean(abs(envmax-envmin),1)/2; end end else [indmin,indmax,indzer] = extr(m); nem = length(indmin)+length(indmax); nzm = length(indzer); [tmin,tmax,mmin,mmax] = boundary_conditions(indmin,indmax,t,m,m,NBSYM); envmin = interp1(tmin,mmin,t,INTERP); envmax = interp1(tmax,mmax,t,INTERP); envmoy = (envmin+envmax)/2; if nargout > 3 amp = mean(abs(envmax-envmin),1)/2; end end end %------------------------------------------------------------------------------- % 默认的停止准则 function [stop,envmoy,s] = stop_sifting(m,t,sd,sd2,tol,INTERP,MODE_COMPLEX,ndirs) try [envmoy,nem,nzm,amp] = mean_and_amplitude(m,t,INTERP,MODE_COMPLEX,ndirs); sx = abs(envmoy)./amp; s = mean(sx); stop = ~((mean(sx > sd) > tol | any(sx > sd2)) & (all(nem > 2))); if ~MODE_COMPLEX stop = stop && ~(abs(nzm-nem)>1); end catch stop = 1; envmoy = zeros(1,length(m)); s = NaN; end end %------------------------------------------------------------------------------- % stopping criterion corresponding to option FIX %停止准则对应选项修复 function [stop,moyenne]= stop_sifting_fixe(t,m,INTERP,MODE_COMPLEX,ndirs) try moyenne = mean_and_amplitude(m,t,INTERP,MODE_COMPLEX,ndirs); stop = 0; catch moyenne = zeros(1,length(m)); stop = 1; end end %------------------------------------------------------------------------------- % stopping criterion corresponding to option FIX_H function [stop,moyenne,stop_count]= stop_sifting_fixe_h(t,m,INTERP,stop_count,FIXE_H,MODE_COMPLEX,ndirs) try [moyenne,nem,nzm] = mean_and_amplitude(m,t,INTERP,MODE_COMPLEX,ndirs); if (all(abs(nzm-nem)>1)) stop = 0; stop_count = 0; else stop_count = stop_count+1; stop = (stop_count == FIXE_H); end catch moyenne = zeros(1,length(m)); stop = 1; end end %------------------------------------------------------------------------------- % displays the progression of the decomposition with the default stopping criterion %显示分解使用默认的进程停止准则 function display_emd(t,m,mp,r,envmin,envmax,envmoy,s,sb,sx,sdt,sd2t,nbit,k,display_sifting,stop_sift) subplot(4,1,1) plot(t,mp);hold on; plot(t,envmax,'--k');plot(t,envmin,'--k');plot(t,envmoy,'r'); title(['IMF ',int2str(k),'; iteration ',int2str(nbit),' before sifting']); set(gca,'XTick',[]) hold off subplot(4,1,2) plot(t,sx) hold on plot(t,sdt,'--r') plot(t,sd2t,':k') title('stop parameter') set(gca,'XTick',[]) hold off subplot(4,1,3) plot(t,m) title(['IMF ',int2str(k),'; iteration ',int2str(nbit),' after sifting']); set(gca,'XTick',[]) subplot(4,1,4); plot(t,r-m) title('residue'); disp(['stop parameter mean value : ',num2str(sb),' before sifting and ',num2str(s),' a