EMD.zip_EMD_EMD matlab_EMD 程序_信号分解

%2007新出来的包含复数的emd函数(端点视作极值点) function [imf,ort,nbits] = emd3(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个极值点，如果采用mask信号，不进入主循环 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 % 如果设定了迭代次数，且保留停止条件|#zeros-#extrema|<=1 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 % 当前模式幅度过小，机器精度就可能引起虚假极值点的出现 if (max(abs(m))) < (1e-10)*(max(abs(x))) % IMF的最大值小于信号最大值的1e-10 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')%查找是否存在变量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 % 演示过程 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 % 筛循环 imf(k,:) = m; if display_sifting disp(['mode ',int2str(k),' stored']) end nbits(k) = nbit; % 记录每个IMF的迭代次数 k = k+1; % IMF计数 r = r - m; % 从原信号中减去提取的IMF nbit = 0; % 单轮迭代次数清0 end % 主循环 % 计入残余信号 if any(r) && ~any(mask) imf(k,:) = r; end % 计数正交指数 ort = io(x,imf); % 关闭图形 if display_sifting close end end %--------------------------------------------------------------------------------------------------- % 测试是否存在足够的极值点（3个）进行分解，极值点个数小于3个则返回1，这是整体停止条件 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 %--------------------------------------------------------------------------------------------------- % 计数包络均值和模式幅度估计值，返回包络均值 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_conditions1(indmin,indmax,t,m); 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_conditions1(indmin,indmax,t,m); 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_conditions1(indmin,indmax,t,m); % 边界延拓 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))); % 停止准则（增加了极值点个数大于2） if ~MODE_COMPLEX stop = stop && ~(abs(nzm-nem)>1); % 对于实信号，要求极值点和过零点的个数相差1 end catch stop = 1; envmoy = zeros(1,length(m)); s = NaN; end end %------------------------------------------------------------------------------- % 针对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 %------------------------------------------------------------------------------- % 针对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 % 极值点与过零点个数相差1后，还要达到指定次数才停止 stop_count = stop_count+1; stop = (stop_count == FIXE_H); end catch moyenne = zeros(1,length(m)); stop = 1; end end %------------------------------------------------------------------------------- % 演示分解过程（默认准则） 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),';