hlmd.rar_LMD分解_hermite插值_lmd_基于Hermite插值的LMD_局域均值分解

function [PF,sig,rem,as,markup]=hlmd(signal,thrsd) %UNTITLED 此处显示有关此函数的摘要 % 此处显示详细说明 di = signal(2,1)-signal(1,1); as=[]; PF=[]; rem = []; sig = []; x=signal(:,1); y=signal(:,2); %pb=plot(x,y,'EraseMode','background','MarkerSize',5); [signal,exp_wid] = ExpandSingal(signal); TNM=fix(log2(length(signal)))-1; while TNM TNM = TNM - 1; if TNM==0 rem = signal(:,2); end nextsignal = signal; [emax,emin,~,~] = OldExtrema(nextsignal,0.0001); if isempty(emax) || isempty(emin) || length(emax(:,1))<3 || length(emin(:,1))<3 rem = signal(:,2); break; end [evnMax,evnMin]=HermitSmooth(emax,emin,nextsignal(:,1)); env = Envelop(evnMax,evnMin); a=ones(length(signal(:,1)),1); tickNum=0; %tic; while(1) %set(pb(1),'XData',env(:,1),'YData',env(:,2)); %axis([env(1,1),env(end,1),-3,3]); %drawnow; %pause(0.5); tickNum = tickNum +1; mean = Mean(evnMax,evnMin); a = a .* env(:,2); h = nextsignal(:,2)-mean(:,2); s = h ./ env(:,2); if any(s>10) %disp('11'); end nextsignal = [nextsignal(:,1),s]; [emax,emin,flag,markup] = OldExtrema(nextsignal,0.0001); if flag==0 [evnMax,evnMin]=HermitSmooth(emax,emin,nextsignal(:,1)); env = Envelop(evnMax,evnMin); bCompared = abs(env(:,2)-1)<thrsd; end if flag==1 || all(bCompared==1) %|| tickNum>20 %|| any(find(env<0.0001)) %curPoint>prePoint % %得到一个PF分量，重新设置信号 %if ~isempty(emax) && ~isempty(emin) pf = a .* nextsignal(:,2); as(1:length(a),end+1)=a; signal(:,2) = signal(:,2) - pf; PF(1:length(pf),end+1)=pf; sig(1:length(nextsignal(:,2)),end+1)=nextsignal(:,2); %tickNum break; %end %return; end %prePoint = curPoint; end %toc; end as = CutSignal(as,exp_wid); PF = CutSignal(PF,exp_wid); sig = CutSignal(sig,exp_wid); rem = CutSignal(rem,exp_wid); end function [mean] = Mean(extrema_max,extrema_min) mean_y=(extrema_max(:,2)+extrema_min(:,2))/2; mean = [extrema_max(:,1),mean_y]; end function [a] = Envelop(extrema_max,extrema_min) a = abs(extrema_max(:,2)-extrema_min(:,2))/2; a=[extrema_max(:,1),a]; end function [exp_sig,exp_wid]=ExpandSingal(sig) %%将信号以左右端点为轴，镜像延拓1/3段数据长度 [row,col] = size(sig); exp_wid = floor(row/3); start=exp_wid+1; stop=row+exp_wid; exp_sig = zeros(row+exp_wid*2,col); exp_sig(start:stop,:)=sig(:,:); exp_sig(1:start-1,2)=flipud(sig(1:start-1,2)); exp_sig(stop+1:end,2)=flipud(sig(row-exp_wid+1:end,2)); dx=sig(2,1)-sig(1,1); for i=1:exp_wid exp_sig(i,1)=sig(1,1)-(exp_wid+1-i)*dx; exp_sig(row+exp_wid+i,1)=sig(end,1)+i*dx; end end function [cutsig]=CutSignal(sig,wid) len = length(sig); cutsig=sig(wid+1:len-wid,:); end function [rmax,rmin]=SmoothAgain(pos,vmax,vmin,x,dist) xLen = length(x); %构造插值区间和端点 if ~isempty(pos) %可能包络的多个部分会有相交,寻找出相交的位置 gaps=diff(pos); startIndexes = find(gaps>1); if isempty(startIndexes) %表明只有一段相交区域 startIndexes = 1; startPos = pos(startIndexes)-1; stopPos = pos(end)+1; else stopIndexes=[startIndexes;length(pos)]; startIndexes=[1;startIndexes+1]; startPos = pos(startIndexes)-1; stopPos = pos(stopIndexes)+1; end if stopPos(end)>length(vmax) stopPos(end)=length(vmax); end if startPos(1)==0 startPos(1)=1; end len = length(startIndexes); for i=1:len ranges = [startPos(i):stopPos(i)]; if ranges(1)==1 || ranges(end)==xLen vmax(ranges,2)=vmax(ranges,2)+dist; %vmin(ranges,2)=mean_min; continue; end try px=[startPos(i);stopPos(i)]; interp_x=x(ranges); y_max=[vmax(startPos(i),2);vmax(stopPos(i),2)]; y_min=[vmin(startPos(i),2);vmin(stopPos(i),2)]; catch disp('111'); end interp_max = interp1(x(px),y_max,interp_x,'linear'); interp_min = interp1(x(px),y_min,interp_x,'linear'); mean_max = (interp_max+vmax(ranges,2))/2; mean_min = (interp_min+vmin(ranges,2))/2; while min(mean_max-mean_min)<=dist if all(interp_max-mean_max==0) && all(mean_min-interp_min==0) disp(sprintf('max distance: %d',max(mean_min)-min(mean_max))); break; end mean_max = (interp_max+mean_max)/2; mean_min = (interp_min+mean_min)/2; end vmax(ranges,2)=mean_max; vmin(ranges,2)=mean_min; end end rmax=vmax; rmin=vmin; end function [vmax,vmin]=HermitSmooth(origMax,origMin,x) vmax=[x,interp1(origMax(:,1),origMax(:,2),x,'pchip')]; vmin = [x,interp1(origMin(:,1),origMin(:,2),x,'pchip')]; dist = 0.01; pos=find(vmax(:,2)-vmin(:,2)<dist); [vmax,vmin]=SmoothAgain(pos,vmax,vmin,x,dist); end function [vmax,vmin,flag,markup]=OldExtrema(signal,p) %p是精度 flag=0; markup=0; len = length(signal(:,1)); presign=diff(signal(1:len,2)); zeropos=find(abs(presign)<p); presign(zeropos)=0; pres = sign(presign(1,1)); extremapos=[1]; for i=2:len-1 if presign(i)==0 continue; end v=sign(presign(i)); if pres~=0 && v*pres<0 %得到极大值位置 extremapos=[extremapos,i]; pres=v; elseif pres==0 pres=v; end end extremapos=[extremapos,len]; if isempty(extremapos) extremapos=[1,len]; end vmin=[]; vmax=[]; extrema_len = length(extremapos); exist_equal=0; start=0; stop=0; %标记起止点是否有相同值 for i=2:extrema_len if abs(signal(extremapos(i),2)-signal(extremapos(i-1),2))<p exist_equal=1; if i==2 start=1; elseif i==extrema_len stop=1; end elseif signal(extremapos(i),2)>signal(extremapos(i-1),2) if exist_equal==1 exist_equal=0; %vmin=[vmin;signal(extremapos(i-1),:)]; end vmin=[vmin;signal(extremapos(i-1),:)]; if i==extrema_len vmax=[vmax;signal(extremapos(i),:)]; end else if exist_equal==1 exist_equal=0; %vmax=[vmax;signal(extremapos(i-1),:)]; end vmax=[vmax;signal(extremapos(i-1),:)]; if i==extrema_len vmin=[vmin;signal(extremapos(i),:)]; end end end mincol=length(vmin); maxcol=length(vmax); if maxcol<=2 || mincol<=2 flag =1; return; end %处理端点上下包络 %直接设置端点为最近的极值点 if start==1 extrema_max_x=[signal(1,1);vmax(:,1)]; extrema_max_y=[vmax(1,2);vmax(:,2)]; vmax=[extrema_max_x';extrema_max_y']'; extrema_min_x=[signal(1,1);vmin(:,1)]; extrema_min_y=[vmin(1,2);vmin(:,2)]; vmin=[extrema_min_x';extrema_min_y']'; else if vmax(1,1)>vmin(1,1) %左端点是极小值 extrema_max_x=[signal(1,1);vmax(:,1)]; extrema_max_y=[vmax(1,2);vmax(:,2)]; vmax=[extrema_max_x';extrema_max_y']'; vmin(1,2)=vmin(2,2); else extrema_min_x=[signal(1,1);vmin(:,1)]; extrema_min_y=[vmin(1,2);vmin(:,2)]; vmin=[extrema_min_x';extrema_min_y']'; vmax(1,2)=vmax(2,2); end end if stop==1 extrema_min_x=[vmin(:,1);signal(end,1)]; extrema_min_y=[vmin(:,2);vmin(end,2)]; vmin=[extrema_min_x';extrema_min_y'