matlab程序设计源程序.docx

1、matlab程序设计源程序哦LMD局域均值分解的matlab程序及示例(2011-03-10 15:55:19) 转载标签： matlab原程序lmd局域均值分解局部均值分解实例程序分类： matlab语言学习 说明：研究LMD局域均值分解有3个月左右，能找到的相关文章也基本上看了一遍，觉得是个很好的方法，号称是EMD经验模态分解的改进版。但是网络上一直没有找到该算法的matlab程序，只见文章说的天花乱坠。后来自己写了一个，但是使用有一个问题没有解决，就是分解的时候怎么去除骑行波的问题。先把自己写的程序贡献出来，让大家分析下，一起讨论下，到底LMD的程序怎样写才能如文献中说的那样达到目的。欢

2、迎大家热烈讨论！程序仍存在不收敛的问题，拿出来分享只是希望高手指点一二，写的不好欢迎拍砖！想要M文件的，给出下载地址：文件夹包含，找出纯调频函数，计算瞬时幅值，瞬时频率的函数等%主程序%lmd1原始版%通过emd.m的三次样条+镜像延拓求出上下包络及均值%局域均值函数=（上包络+下包络）/2%局域包络函数=|上包络-下包络|/2%相关文章见一种基于EMD的振动信号时频分析新方法研究-胡劲松，杨世锡functionPF,A,SI=lmd1(m)%最后一个PF是残余分量%A是瞬时赋值%SI是纯调频函数，求它的瞬时频率就是需要的频率c=m;k=0;wucha1=0.001;n_l=nengliang

3、(m);while 1 k=k+1; a=1; h=c; pf,a,si=zhaochun1(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); %停止调节 %1.emd用的是三次样条求包络，要求至少3个极值点，所以这里c的极值点个数也应该至少为3 %2.如果上一个PF的极值点数比下一个PF的极值点数少，说明结果也不正确（这个也可以作为停止条件考虑进去） %上面一句是否可以等价于当前PF的极值点个数一定要大于等于残量（c）的极值点个数

4、（目前是用这个作为停止条件的一个参考写入程序） %3.当前PF分量的能量应该大于残量c的能量(这个有待商榷) %4.残余能量不能大于信号能量 if n_pfn_c if length(c_pos)3 | n_cn_l/10 |k=3%| n_pfn_l PF(k+1,:)=c; break endendend%nengliang函数%function p=nengliang(y)% my=mean(y);% p=(y-my).*(y-my);% p=sum(p);p=sum(abs(y).2);end%找纯函数%function pf,a,si=zhaochun1(a,h,wucha1)chu

5、n_num=0;while 1chun_num=chun_num+1;t=1:length(h);envmin,envmax,envmoy,indmin,indmax,indzer = envelope(t,h,spline);mi=(envmax+envmin)./2;ai=abs(envmax-envmin)./2;a=a.*ai;si=(h-mi)./ai;h=si;ai_funmax=max(ai);ai_funmin=min(ai);if (ai_funmax=1-wucha1) breakendendpf=a.*si;chun_numfunction envmin, envmax,

6、envmoy,indmin,indmax,indzer = envelope(t,x,INTERP)%computes envelopes and mean with various interpolationsNBSYM = 2;% 边界延拓点数DEF_INTERP = spline;if nargin 1error(x and t must be vectors)ends = size(x);if s(1) 1x = x;ends = size(t);if s(1) 1t = t;endif length(t) = length(x)error(x and t must have the

7、same length)endlx = length(x);indmin,indmax,indzer = extr(x,t);%boundary conditions for interpolationtmin,tmax,xmin,xmax = boundary_conditions(indmin,indmax,t,x,NBSYM);% definition of envelopes from interpolationenvmax = interp1(tmax,xmax,t,INTERP);envmin = interp1(tmin,xmin,t,INTERP);if nargout 2 e

8、nvmoy = (envmax + envmin)/2;endendfunction 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) x(

9、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;endelseif x(1) x(indmax(1)% 以第一个极小值为对称中心lmax = fl

10、iplr(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;endend % 序列末尾情况与序列开头类似if indmax(end) indmin(end)if x(end) x(indmin(end)rmax = fli

11、plr(indmax(max(end-nbsym,1):end-1);rmin = fliplr(indmin(max(end-nbsym+1,1):end);rsym = indmax(end);elsermax = fliplr(indmax(max(end-nbsym+1,1):end);rmin = lx,fliplr(indmin(max(end-nbsym+2,1):end);rsym = lx;endend % 将序列根据对称中心，镜像到两边tlmin = 2*t(lsym)-t(lmin);tlmax = 2*t(lsym)-t(lmax);trmin = 2*t(rsym)-

12、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);elselmin = fliplr(indmin(1:min(end,nbsym);endif lsym = 1% 这种情况不应该出现，程序直接中止error(bug)endlsym =

13、1;% 直接关于第一采样点取镜像tlmin = 2*t(lsym)-t(lmin);tlmax = 2*t(lsym)-t(lmax);end % 序列末尾情况与序列开头类似if trmin(end) t(lx) | trmax(end) 2x1=x(1:m-1);x2=x(2:m);indzer = find(x1.*x20);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)

14、-1; indz = round(debz+finz)/2); else indz = iz; end indzer = sort(indzer indz);endendd = diff(x);n = length(d);d1 = d(1:n-1);d2 = d(2:n);indmin = find(d1.*d20 & d10)+1;indmax = find(d1.*d20)+1;% when two or more consecutive points have the same value we consider only one extremum in the middle of th

15、e 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

16、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 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 =

17、sort(indmin imin); endendendend%pos%function poss=pos(y)%一个输出%功能：找序列极值点位置坐标%y:输入序列%pos:极值点的序列位置坐标indmin,indmax=position(y);minmax=cat(2,indmin,indmax);poss=sort(minmax);end%position%返还极大值和极小值位置下标%两个输出function indmin,indmax=position(y)m = length(y);d = diff(y);n = length(d);d1 = d(1:n-1);d2 = d(2:n);

18、indmin = find(d1.*d20 & d10)+1;indmax = find(d1.*d20)+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

19、 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 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 indmi

20、n = sort(indmin imin); endendend%说明每个程序要单独保存成m文件，放在同一个文件夹下调用使用实例：x=(t) (2+cos(90*t).*cos(500*t+1800.*t.*t);fs=5000;t=0:1/fs:0.341;y=x(t);subplot(5,1,1);plot(t,y);xlabel(原始信号);pf,a,si=lmd1(y);subplot(5,1,2);plot(t,pf(1,:);xlabel(PF1);subplot(5,1,3);plot(t,pf(2,:);xlabel(PF2);subplot(5,1,4);plot(t,pf(3,:);xlabel(PF3);subplot(5,1,5);plot(t,pf(4,:);xlabel(残量信号);从上图可以看出，成功将第一个分量和第二个分量分离出来，但是残余分量存在很大的问题，这是因为分解过程中的骑行波没有去处导致的，至于怎样完善LMD局域均值分解算法，目前个人没有时间研究了，希望得到指点。