lmd系列.rar_lmd

function[PF,A,SI]=lmd(m) c=m; k=0 wucha1=0.001; n_l=nengliang(m); while 1 k=k+1; a=1; h=c; [pf,a,si]=zhaochun(a,h,wucha1); c=c-pf; PF(k,:)=pf; A(k,:)=a; SI(k,:)=si; c_pos=pos(c); n_c=nengliang(c); n_pf=nengliang(pf); if length(c_pos)<3 || n_c<n_l/100 || length(pos(pf))<length(c_pos) || n_pf<n_c PF(k+1,:)=c; break end end function pos=pos(y) %功能：找序列极值点位置坐标 %y:输入序列 %pos:极值点的序列位置坐标 m = length(y); d = diff(y); 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; 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 minmax=cat(2,indmin,indmax); pos=sort(minmax); end %----------zhaochun.m function [pf,a,si]=zhaochun(a,h,wucha1) chun_num=0; while 1 chun_num=chun_num+1 t=1:length(h); h_pos=position(h);%极值点位置序列 tpoint=t(h_pos);%极值点时间值 hpoint=h(h_pos);%极值点幅度值 hpoint=bianjie(h,hpoint,1);%端点处理后的极值点，多出了2个极值点 [mi,ai]=find_miai(hpoint);%找出极值点之间的均值函数和包络函数 mi_point=junzhi(mi);%所有极值点处的均值序列（幅值）-纵坐标（点序列） ai_point=junzhi(ai);%所有极值点处的包络序列（幅值）-纵坐标（点序列） %此处端点处的均值和包络都是 端点和极值点之间的均值和包络值（如果端点视作极值点，这样处理是不合理的，端点都不是极值点，这样处理是正确的） lmi_point=mi(1);%左端点的均值 rmi_point=mi(length(mi));%右端点的均值 lai_point=ai(1);%左端点的包络 rai_point=ai(length(ai));%右端点的包络 % mi_point_d=link(lmi_point,rmi_point,mi_point);%连接端点均值及所有极值点处的 均值 （带端点的均值序列）（点序列） ai_point_d=link(lai_point,rai_point,ai_point);%连接端点包络及所有极值点出的 包络值（带端点的包络序列）（点序列） % tpoint_d=link(t(1),t(length(t)),tpoint); mi_fun=chadian1(tpoint_d,mi,mi_point_d);%包含端点和极值点和普通点的 均值序列-平缓前的均值序列 ai_fun=chadian1(tpoint_d,ai,ai_point_d);%包含端点和极值点和普通点的 包络序列-平滑前的包络序列 %以上是完整的未处理的均值函数mi_fun和包络函数ai_fun %找出第一平滑的滑动跨度 kmax=max(diff(tpoint_d));%找出时间跨度最大的 相邻几点 间的 距离 kmax1=uint16(kmax/3); kmax1=double(kmax1); jiou=uint8(rem(kmax1,2)); if kmax1<3 move_k=3; % b<3,滑动跨度=3， elseif jiou==0 % b>=3,当b是偶数时，跨度=b-1 move_k=(kmax1-1); else move_k=kmax1; %b>=3,b是奇数时，跨度=b end [mi_move,move_mi_num]=move(mi_fun,move_k); [ai_move,move_ai_num]=move(ai_fun,move_k); mi=mi_move; ai=ai_move; 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; end function [x,flag]=move(a,k) l=length(a); t=1:l; % jingdu=t(2)-t(1); flag=0; x=a; max_flag=10;%平滑重复的最大次数设置 max_error=0;%相邻两点间相等允许的最大差异（理论上=0才人为无差异，实际操作要考虑采样精度，所以可以设置一个允许范围） while flag<=max_flag || min(abs(diff(x)))<=max_error;%这里用到abs，然后再min，因为幅值的差值会出现负值，我们的目的是找差值=0，而不是负数 if flag==0 flag=flag+1;%flag标志位，标志平滑的次数，这里这里设置为不超过11次（最大10次） else flag=flag+1;%flag标志位，标志平滑的次数，这里这里设置为不超过11次（最大10次） x_pos=position(x);%极值点位置序列 tpoint=t(x_pos);%极值点时间值 tpoint_d=link(t(1),t(l),tpoint);%极值点的时间轴上的位置 kmax=max(diff(tpoint_d));%找出时间跨度最大的 相邻几点 间的 距离 kmax1=uint16(kmax/3); kmax1=double(kmax1); jiou=uint8(rem(kmax1,2)); if kmax1<3 k=3; % b<3,滑动跨度=3， elseif jiou==0 % b>=3,当b是偶数时，跨度=b-1 k=(kmax1-1); else k=kmax1; %b>=3,b是奇数时，跨度=b end end y=zeros(1,l);% 初始化序列y, y是中间变量 %%%%%%%%%%%%%%%%%%%边界点的处理%%%%%%%%%%%%%% for i=1:(k-1)/2 v=0; z=0; for j=1:i*2-1 v=v+x(j); end y(i)=v/(i*2-1); for j=l:-1:l+2-i*2 %j=l:-1:l+1-(i*2-1) z=z+x(j); end y(l+1-i)=z/(i*2-1); end %%%%%%%%%%%%%%中间点的处理 %%%%%%%%%%%%%%%%%%%%%%%%%% for i=1+(k-1)/2:l-(k-1)/2 %这个 (d+1)/2是跨度的中点 单边的边界点数=中点值-1=(d+1)/2-1=(d-1)/2 w=x(i); for j=1:(k-1)/2 w=w+x(i-j)+x(i+j); end y(i)=w/k; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% x=y; end %k % 调试的时候查看 %flag % 调试的时候查看 end function hpp=bianjie(hh,hp,style) %hh:需要边界处理的序列 %hp:需要边界处理序列的极值点的幅度值 %style:边界处理类型1：端点全部看作极值点2：左右两端各增加1个幅值为0的极值点 %hpp:边界处理后，极值点幅值多处了两个，因为左边右边各人为加上了一个 if style==1 a=hh(1); b=hh(length(hh)); end if style==2 a=0; b=0; end hpp=link(a,b,hp); end function pos=position(y) %功能：找序列极值点位置坐标 %y:输入序列 %pos:极值点的序列位置坐标 m = length(y); d = diff(y); 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; 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 minmax=cat(2,indmin,indmax); pos=sort(minmax); end function [mi,ai]=find_miai(ypoint) % Lyp=length(ypoint); al=wkeep(ypoint,Lyp-1,'l'); ar=wkeep(ypoint,Lyp-1,'r'); mi=(al+ar)/2; ai=abs(ar-al)/2; end function c=junzhi(a) l=length(a); al=wkeep(a,l-1,'l'); ar=wkeep(a,l-1,'r'); c=(al+ar)/2; end function d=link(a,b,c) d=[a';c';b']'; %头：尾：中间 end function f_value=chadian1(a,b,c) % chadian1 把端点及极值点处对应的总坐标值插入，原来的均值函数的方波序列中 %输入参数a:点序列（行向量）（包含端点和极值点以 在时间上的位置-横坐标）（n个）（点序列-横坐标） %输入参数b:段序列（行向量）（点序列a 每两点之间的纵坐标的值-纵坐标） （n-1）-（段序列-纵坐标�