LMD经验模态分解matlab程序

LMD经验模态分解matlab程序——原味的曾经也用滑动平均写过LMD，其实滑动平均的EMD才是原汁原味的居于均值分解。分享给有需要的人，程序写的不好，只是希望提供一种思路。如果谁写了更完美LMD程序，别忘了发我一份，快毕业了，一直没有把LMD写完美，对于我来说始终是个遗憾。来分完美的LMD让我也品尝下，我也无憾了~代码下载地址：此处没有提供测试代码，如需要可以点这里：点我源代码如下：%原始lmd算法，效果很不好，不知道程序哪里写错function[PF,A,SI]=lmd(m)c=m;k=0wucha1=0.001;n_l=nengliang(m);while1k=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);iflength(c_pos)3||n_cn_l/100||length(pos(pf))length(c_pos)||n_pfn_cPF(k+1,:)=c;breakendendfunctionpos=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.*d20&d10)+1;indmax=find(d1.*d20&d10)+1;ifany(d==0)imax=[];imin=[];bad=(d==0);dd=diff([0bad0]);debs=find(dd==1);fins=find(dd==-1);ifdebs(1)==1iflength(debs)1debs=debs(2:end);fins=fins(2:end);elsedebs=[];fins=[];endendiflength(debs)0iffins(end)==miflength(debs)1debs=debs(1:(end-1));fins=fins(1:(end-1));elsedebs=[];fins=[];endendendlc=length(debs);iflc0fork=1:lcifd(debs(k)-1)0ifd(fins(k))0imax=[imaxround((fins(k)+debs(k))/2)];endelseifd(fins(k))0imin=[iminround((fins(k)+debs(k))/2)];endendendendiflength(imax)0indmax=sort([indmaximax]);endiflength(imin)0indmin=sort([indminimin]);endendminmax=cat(2,indmin,indmax);pos=sort(minmax);end%----------zhaochun.mfunction[pf,a,si]=zhaochun(a,h,wucha1)chun_num=0;while1chun_num=chun_num+1t=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));ifkmax13move_k=3;%b3,滑动跨度=3，elseifjiou==0%b=3,当b是偶数时，跨度=b-1move_k=(kmax1-1);elsemove_k=kmax1;%b=3,b是奇数时，跨度=bend[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)breakendendpf=a.*si;endfunction[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才人为无差异，实际操作要考虑采样精度，所以可以设置一个允许范围）whileflag=max_flag||min(abs(diff(x)))=max_error;%这里用到abs，然后再min，因为幅值的差值会出现负值，我们的目的是找差值=0，而不是负数ifflag==0flag=flag+1;%flag标志位，标志平滑的次数，这里这里设置为不超过11次（最大10次）elseflag=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));ifkmax13k=3;%b3,滑动跨度=3，elseifjiou==0%b=3,当b是偶数时，跨度=b-1k=(kmax1-1);elsek=kmax1;%b=3,b是奇数时，跨度=bendendy=zeros(1,l);%初始化序列y,y是中间变量%%%%%%%%%%%%%%%%%%%边界点的处理%%%%%%%%%%%%%%fori=1:(k-1)/2v=0;z=0;forj=1:i*2-1v=v+x(j);endy(i)=v/(i*2-1);forj=l:-1:l+2-i*2%j=l:-1:l+1-(i*2-1)z=z+x(j);endy(l+1-i)=z/(i*2-1);end%%%%%%%%%%%%%%中间点的处理%%%%%%%%%%%%%%%%%%%%%%%%%%fori=1+(k-1)/2:l-(k-1)/2%这个(d+1)/2是跨度的中点单边的边界点数=中点值-1=(d+1)/2-1=(d-1)/2w=x(i);forj=1:(k-1)/2w=w+x(i-j)+x(i+j);endy(i)=w/k;end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%x=y;end%k%调试的时候查看%flag%调试的时候查看endfunctionhpp=bianjie(hh,hp,style)%hh:需要边界处理的序列%hp:需要边界处理序列的极值点的幅度值%style:边界处理类型1：端点全部看作极值点2：左右两端各增加1个幅值为0的极值点%hpp:边界处理后，极值点幅值多处了两个，因为左边右边各人为加上了一个ifstyle==1a=hh(1);b=hh(length(hh));endifstyle==2a=0;b=0;endhpp=link(a,b,hp);endfunctionpos=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.*d20&d10)+1;indmax=find(d1.*d20&d10)+1;ifany(d==0)imax=[];imin=[];bad=(d==0);dd=diff([0bad0]);debs=find(dd==1);fins=find(dd==-1);ifdebs(1)==1iflength(debs)1debs=debs(2:end);fins=fins(2:end);elsedebs=[];fins=[];endendiflength(debs)0iffins(end)==miflength(debs)1debs=debs(1:(end-1));fins=fins(1:(end-1));elsedebs=[];fins=[];endendendlc=length(debs);iflc0fork=1:lcifd(debs(k)-1)0ifd(fins(k))0imax=[imaxround((fins(k)+debs(k))/2)];endelseifd(fins(k))0imin=[iminround((fins(k)+debs(k))/2)];endendendendiflength(imax)0indmax=sort([indmaximax]);endiflength(imin)0indmin=sort([indminimin]);endendminmax=cat(2,indmin,indmax);pos=sort(minmax);endfunction[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;endfunctionc=junzhi(a)l=length(a);al=wke