镜像延拓程序-EMD端点效应_EMD_端点效应

FS =100; T = 1/FS; L =2000; t = (0:L-1)*T; t1=(0:L-1)*T; f1=0.25; f2=1.05; x1=2*cos(2*pi*f1*t+1/3*pi); [indmin, indmax, indzer] = extr(x1,t); z=x1 [tmin,tmax,zmin,zmax] = boundary_conditions(indmin,indmax,t,x1,z,1) %======================边界延拓镜像法============================================= % 处理边界条件（镜像法） function [tmin,tmax,zmin,zmax] = boundary_conditions(indmin,indmax,t,x,z,nbsym) % 实数情况下，x = z 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); % 如果对称的部分没有足够的极值点 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 % 延拓点上的取值 zlmax = z(lmax); zlmin = z(lmin); zrmax = z(rmax); zrmin = z(rmin); % 完成延拓 tmin = [tlmin t(indmin) trmin]; tmax = [tlmax t(indmax) trmax]; zmin = [zlmin z(indmin) zrmin]; zmax = [zlmax z(indmax) zrmax]; end %=====================极值点和过零点位置提取========================================== function [indmin, indmax, indzer] = extr(x,t) 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) % 考虑信号采样点恰好为0的位置 iz = find( x==0 ); % 信号采样点恰好为0的位置 indz = []; if any(diff(iz)==1) % 出现连0的情况 zer = x == 0; % x=0处为1，其它地方为0 dz = diff([0 zer 0]); % 寻找0与非0的过渡点 debz = find(dz == 1); % 0值起点 finz = find(dz == -1)-1; % 0值终点 indz = round((debz+finz)/2); % 选择中间点作为过零点 else indz = iz; % 若没有连0的情况，该点本身就是过零点 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; % 最大值 % 当连续多个采样值相同时，把最中间的一个值作为极值点，处理方式与连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