LCD.zip_LCD 分解_LCD 局部特征_LCD数据分解_局部分解_局部分解尺度

function [ISC,chun_num]=zhaochun3(a,h,wucha1) dg=0.01;; chun_num=0; while 1 chun_num=chun_num+1; t=1:length(h); [H] = envelope(t,h,'spline');%求出H，tx,xt; %H=MAPH(H,11); ISC=h-H; [yes]=panduan(ISC); P=[]; if yes==1 break else if chun_num>3 break else h=ISC; end end end function [H] = 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 %if (length(indmin) + length(indmax) < 3) %error('not enough extrema') %break %else [tmin,tmax,xmin,xmax] = boundary_conditions(indmin,indmax,t,x,NBSYM); % definition of envelopes from interpolation tx=sort([tmin,tmax]); xt=[]; if tmin(1)<tmax(1) if length(tmin)==length(tmax) k=1; for i=1:length(tmin) xt(k)=xmin(i); xt(k+1)=xmax(i); k=k+2; end else k=1; for i=1:length(tmax) xt(k)=xmin(i); xt(k+1)=xmax(i); k=k+2; end xt(k)=xmin(i+1); end else if length(tmin)==length(tmax) k=1; for i=1:length(tmax) xt(k)=xmax(i); xt(k+1)=xmin(i); k=k+2; end else k=1; for i=1:length(tmin) xt(k)=xmax(i); xt(k+1)=xmin(i); k=k+2; end xt(k)=xmax(i+1); end end H=[]; L=[]; lk=1; A1=xt(3)+((tx(4)-tx(3))/(tx(5)-tx(3)))*(xt(5)-xt(3)); A2=xt(4)+((tx(5)-tx(4))/(tx(6)-tx(4)))*(xt(6)-xt(4)); L1=0.5*A1+0.5*xt(4); L2=0.5*A2+0.5*xt(5); L(lk)=L1; lk=lk+1; xishu=(L2-L1)/(xt(5)-xt(4)); for i=1:tx(5) H(i)=L1+xishu*(x(i)-xt(4)); end for k=1:length(tx)-9 A1=xt(3+k)+((tx(4+k)-tx(3+k))/(tx(5+k)-tx(3+k)))*(xt(5+k)-xt(3+k)); A2=xt(4+k)+((tx(5+k)-tx(4+k))/(tx(6+k)-tx(4+k)))*(xt(6+k)-xt(4+k)); L1=0.5*A1+0.5*xt(4+k); L2=0.5*A2+0.5*xt(5+k); L(lk)=L1; lk=lk+1; xishu=(L2-L1)/(xt(5+k)-xt(4+k)); for i=tx(k+4):tx(k+5) H(i)=L1+xishu*(x(i)-xt(4+k)); end end k=length(tx)-8; A1=xt(3+k)+((tx(4+k)-tx(3+k))/(tx(5+k)-tx(3+k)))*(xt(5+k)-xt(3+k)); A2=xt(4+k)+((tx(5+k)-tx(4+k))/(tx(6+k)-tx(4+k)))*(xt(6+k)-xt(4+k)); L1=0.5*A1+0.5*xt(4+k); L2=0.5*A2+0.5*xt(5+k); L(lk)=L1; lk=lk+1; L(lk)=L2; xishu=(L2-L1)/(xt(5+k)-xt(4+k)); for i=tx(k+4):length(x) H(i)=L1+xishu*(x(i)-xt(4+k)); end H=0; env = interp1(tx(4:length(tx)-3),L,t,INTERP); H=env;%用三次样条代替线性变换 %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);