lmd.rar_decomposition_lmd_mode decomposition_should5nl_slipped6e

function [pf,a,si]=zhaochun1(a,h,wucha1) chun_num=0; while 1 chun_num=chun_num+1; t=1:length(h); [envmin,envmax,envmoy,indmin,indmax,indzer] = envelope(t,h,'spline'); mi=(envmax+envmin)./2; ai=abs(envmax-envmin)./2; a=a.*ai; si=(h-mi)./ai; h=si; ai_funmax=max(ai); ai_funmin=min(ai); if (ai_funmax<=1+wucha1&&ai_funmin>=1-wucha1) break end end pf=a.*si; chun_num function [envmin, envmax,envmoy,indmin,indmax,indzer] = envelope(t,x,INTERP) %computes envelopes and mean with various interpolations NBSYM = 2; % 边界延拓点数 DEF_INTERP = 'spline'; if nargin < 2 x = t; t = 1:length(x); INTERP = DEF_INTERP; end if nargin == 2 if ischar(x) INTERP = x; x = t; t = 1:length(x); end end if ~ischar(INTERP) error('interp parameter must be ''linear'''', ''cubic'' or ''spline''') end if ~any(strcmpi(INTERP,{'linear','cubic','spline'})) error('interp parameter must be ''linear'''', ''cubic'' or ''spline''') end if min([size(x),size(t)]) > 1 error('x and t must be vectors') end s = size(x); if s(1) > 1 x = x'; end s = size(t); if s(1) > 1 t = t'; end if length(t) ~= length(x) error('x and t must have the same length') end lx = length(x); [indmin,indmax,indzer] = extr(x,t); %boundary conditions for interpolation [tmin,tmax,xmin,xmax] = boundary_conditions(indmin,indmax,t,x,NBSYM); % definition of envelopes from interpolation envmax = interp1(tmax,xmax,t,INTERP); envmin = interp1(tmin,xmin,t,INTERP); if nargout > 2 envmoy = (envmax + envmin)/2; end end function [tmin,tmax,xmin,xmax] = boundary_conditions(indmin,indmax,t,x,nbsym) % computes the boundary conditions for interpolation (mainly mirror symmetry) lx = length(x); % 判断极值点个数 if (length(indmin) + length(indmax) < 3) error('not enough extrema') end % 插值的边界条件 if indmax(1) < indmin(1)% 第一个极值点是极大值 if x(1) > x(indmin(1))% 以第一个极大值为对称中心 lmax = fliplr(indmax(2:min(end,nbsym+1))); lmin = fliplr(indmin(1:min(end,nbsym))); lsym = indmax(1); else% 如果第一个采样值小于第一个极小值，则将认为该值是一个极小值，以该点为对称中心 lmax = fliplr(indmax(1:min(end,nbsym))); lmin = [fliplr(indmin(1:min(end,nbsym-1))),1]; lsym = 1; end else if x(1) < x(indmax(1))% 以第一个极小值为对称中心 lmax = fliplr(indmax(1:min(end,nbsym))); lmin = fliplr(indmin(2:min(end,nbsym+1))); lsym = indmin(1); else% 如果第一个采样值大于第一个极大值，则将认为该值是一个极大值，以该点为对称中心 lmax = [fliplr(indmax(1:min(end,nbsym-1))),1]; lmin = fliplr(indmin(1:min(end,nbsym))); lsym = 1; end end % 序列末尾情况与序列开头类似 if indmax(end) < indmin(end) if x(end) < x(indmax(end)) rmax = fliplr(indmax(max(end-nbsym+1,1):end)); rmin = fliplr(indmin(max(end-nbsym,1):end-1)); rsym = indmin(end); else rmax = [lx,fliplr(indmax(max(end-nbsym+2,1):end))]; rmin = fliplr(indmin(max(end-nbsym+1,1):end)); rsym = lx; end else if x(end) > x(indmin(end)) rmax = fliplr(indmax(max(end-nbsym,1):end-1)); rmin = fliplr(indmin(max(end-nbsym+1,1):end)); rsym = indmax(end); else rmax = fliplr(indmax(max(end-nbsym+1,1):end)); rmin = [lx,fliplr(indmin(max(end-nbsym+2,1):end))]; rsym = lx; end end % 将序列根据对称中心，镜像到两边 tlmin = 2*t(lsym)-t(lmin); tlmax = 2*t(lsym)-t(lmax); trmin = 2*t(rsym)-t(rmin); trmax = 2*t(rsym)-t(rmax); % in case symmetrized parts do not extend enough% 如果对称的部分没有足够的极值点 if tlmin(1) > t(1) | tlmax(1) > t(1)% 对折后的序列没有超出原序列的范围 if lsym == indmax(1) lmax = fliplr(indmax(1:min(end,nbsym))); else lmin = fliplr(indmin(1:min(end,nbsym))); end if lsym == 1% 这种情况不应该出现，程序直接中止 error('bug') end lsym = 1; % 直接关于第一采样点取镜像 tlmin = 2*t(lsym)-t(lmin); tlmax = 2*t(lsym)-t(lmax); end % 序列末尾情况与序列开头类似 if trmin(end) < t(lx) | trmax(end) < t(lx) if rsym == indmax(end) rmax = fliplr(indmax(max(end-nbsym+1,1):end)); else rmin = fliplr(indmin(max(end-nbsym+1,1):end)); end if rsym == lx error('bug') end rsym = lx; trmin = 2*t(rsym)-t(rmin); trmax = 2*t(rsym)-t(rmax); end % 延拓点上的取值 xlmax =x(lmax); xlmin =x(lmin); xrmax =x(rmax); xrmin =x(rmin); % 完成延拓 tmin = [tlmin t(indmin) trmin]; tmax = [tlmax t(indmax) trmax]; xmin = [xlmin x(indmin) xrmin]; xmax = [xlmax x(indmax) xrmax]; end %--------------------------------------------------------------------------------------------------- % 极值点和过零点位置提取 function [indmin, indmax, indzer] = extr(x,t); %extracts the indices corresponding to extrema if(nargin==1) t=1:length(x); end m = length(x); if nargout > 2 x1=x(1:m-1); x2=x(2:m); indzer = find(x1.*x2<0); if any(x == 0) iz = find( x==0 ); indz = []; if any(diff(iz)==1) zer = x == 0; dz = diff([0 zer 0]); debz = find(dz == 1); finz = find(dz == -1)-1; indz = round((debz+finz)/2); else indz = iz; end indzer = sort([indzer indz]); end end d = diff(x); n = length(d); d1 = d(1:n-1); d2 = d(2:n); indmin = find(d1.*d2<0 & d1<0)+1; indmax = find(d1.*d2<0 & d1>0)+1; % when two or more consecutive points have the same value we consider only one extremum in the middle of the constant area % 当连续多个采样值相同时，把最中间的一个值作为极值点，处理方式与连0类似 if any(d==0) imax = []; imin = []; bad = (d==0); dd = diff([0 bad 0]); debs = find(dd == 1); fins = find(dd == -1); if debs(1) == 1 if length(debs) > 1 debs = debs(2:end); fins = fins(2:end); else debs = []; fins = []; end end if length(debs) > 0 if fins(end) == m if length(debs) > 1 debs = debs(1:(end-1)); fins = fins(1:(end-1)); else debs = []; fins = []; end end end lc = length(debs); if lc > 0 for k = 1:lc if d(debs(k)-1) > 0 if d(fins(k)) < 0 imax = [imax round((fins(k)+debs(k))/2)]; end else if d(fins(k)) > 0 imin = [imin round((fins(k)+debs(k))/2)]; end end end end if length(imax) > 0 indmax = sort([indmax imax]); end if length(imin) > 0 indmin = sort([indmin imin]); end end end end