matlab程序设计源程序.docx
- 文档编号:10327043
- 上传时间:2023-02-10
- 格式:DOCX
- 页数:13
- 大小:47.42KB
matlab程序设计源程序.docx
《matlab程序设计源程序.docx》由会员分享,可在线阅读,更多相关《matlab程序设计源程序.docx(13页珍藏版)》请在冰豆网上搜索。
matlab程序设计源程序
哦
LMD局域均值分解的matlab程序及示例
(2011-03-1015:
55:
19)
转载
标签:
matlab
原程序
lmd
局域均值分解
局部均值分解
实例
程序
分类:
matlab语言学习
说明:
研究LMD局域均值分解有3个月左右,能找到的相关文章也基本上看了一遍,觉得是个很好的方法,号称是EMD经验模态分解的改进版。
但是网络上一直没有找到该算法的matlab程序,只见文章说的天花乱坠。
后来自己写了一个,但是使用有一个问题没有解决,就是分解的时候怎么去除骑行波的问题。
先把自己写的程序贡献出来,让大家分析下,一起讨论下,到底LMD的程序怎样写才能如文献中说的那样达到目的。
欢迎大家热烈讨论!
程序仍存在不收敛的问题,拿出来分享只是希望高手指点一二,写的不好欢迎拍砖!
想要M文件的,给出下载地址:
文件夹包含,找出纯调频函数,计算瞬时幅值,瞬时频率的函数等
%%%%%%%%%%%主程序%%%%%%%%%%
%lmd1原始版
%通过emd.m的三次样条+镜像延拓求出上下包络及均值
%局域均值函数=(上包络+下包络)/2
%局域包络函数=|上包络-下包络|/2
%相关文章见《一种基于EMD的振动信号时频分析新方法研究》-胡劲松,杨世锡
function[PF,A,SI]=lmd1(m)
%最后一个PF是残余分量
%A是瞬时赋值
%SI是纯调频函数,求它的瞬时频率就是需要的频率
c=m;
k=0;
wucha1=0.001;
n_l=nengliang(m);
while1
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)的极值点个数(目前是用这个作为停止条件的一个参考写入程序)
%3.当前PF分量的能量应该大于残量c的能量(这个有待商榷)
%4.残余能量不能大于信号能量
ifn_pf iflength(c_pos)<3||n_c PF(k+1,: )=c; break end end end %%%%%%%%%%%%%%%%%%nengliang函数%%%%%%%%%% functionp=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) chun_num=0; while1 chun_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&&ai_funmin>=1-wucha1) break end end pf=a.*si; chun_num function[envmin,envmax,envmoy,indmin,indmax,indzer]=envelope(t,x,INTERP) %computesenvelopesandmeanwithvariousinterpolations NBSYM=2; %边界延拓点数 DEF_INTERP='spline'; ifnargin<2 x=t; t=1: length(x); INTERP=DEF_INTERP; end ifnargin==2 ifischar(x) INTERP=x; x=t; t=1: length(x); end end if~ischar(INTERP) error('interpparametermustbe''linear'''',''cubic''or''spline''') end if~any(strcmpi(INTERP,{'linear','cubic','spline'})) error('interpparametermustbe''linear'''',''cubic''or''spline''') end ifmin([size(x),size(t)])>1 error('xandtmustbevectors') end s=size(x); ifs (1)>1 x=x'; end s=size(t); ifs (1)>1 t=t'; end iflength(t)~=length(x) error('xandtmusthavethesamelength') end lx=length(x); [indmin,indmax,indzer]=extr(x,t); %boundaryconditionsforinterpolation [tmin,tmax,xmin,xmax]=boundary_conditions(indmin,indmax,t,x,NBSYM); %definitionofenvelopesfrominterpolation envmax=interp1(tmax,xmax,t,INTERP); envmin=interp1(tmin,xmin,t,INTERP); ifnargout>2 envmoy=(envmax+envmin)/2; end end function[tmin,tmax,xmin,xmax]=boundary_conditions(indmin,indmax,t,x,nbsym) %computestheboundaryconditionsforinterpolation(mainlymirrorsymmetry) lx=length(x); %判断极值点个数 if(length(indmin)+length(indmax)<3) error('notenoughextrema') end %插值的边界条件 ifindmax (1) (1)%第一个极值点是极大值 ifx (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 ifx (1) (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 %序列末尾情况与序列开头类似 ifindmax(end) ifx(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 ifx(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); %incasesymmetrizedpartsdonotextendenough%如果对称的部分没有足够的极值点 iftlmin (1)>t (1)|tlmax (1)>t (1)%对折后的序列没有超出原序列的范围 iflsym==indmax (1) lmax=fliplr(indmax(1: min(end,nbsym))); else lmin=fliplr(indmin(1: min(end,nbsym))); end iflsym==1%这种情况不应该出现,程序直接中止 error('bug') end lsym=1; %直接关于第一采样点取镜像 tlmin=2*t(lsym)-t(lmin); tlmax=2*t(lsym)-t(lmax); end %序列末尾情况与序列开头类似 iftrmin(end) ifrsym==indmax(end) rmax=fliplr(indmax(max(end-nbsym+1,1): end)); else rmin=fliplr(indmin(max(end-nbsym+1,1): end)); end ifrsym==lx error('bug') end rsym=lx; trmin=2*t(rsym)-t(rmin); trmax=2*t(rsym)-t(rmax); end %延拓点上的取值 xlmax=x(lmax); xlmin=x(lmin); xrmax=x(rmax); xrmin=x(rmin); %完成延拓 tmin=[tlmint(indmin)trmin]; tmax=[tlmaxt(indmax)trmax]; xmin=[xlminx(indmin)xrmin]; xmax=[xlmaxx(indmax)xrmax]; end %--------------------------------------------------------------------------------------------------- %极值点和过零点位置提取 function[indmin,indmax,indzer]=extr(x,t); %extractstheindicescorrespondingtoextrema if(nargin==1) t=1: length(x); end m=length(x); ifnargout>2 x1=x(1: m-1); x2=x(2: m); indzer=find(x1.*x2<0); ifany(x==0) iz=find(x==0); indz=[]; ifany(diff(iz)==1) zer=x==0; dz=diff([0zer0]); debz=find(dz==1); finz=find(dz==-1)-1; indz=round((debz+finz)/2); else indz=iz; end indzer=sort([indzerindz]); 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; %whentwoormoreconsecutivepointshavethesamevalueweconsideronlyoneextremuminthemiddleoftheconstantarea %当连续多个采样值相同时,把最中间的一个值作为极值点,处理方式与连0类似 ifany(d==0) imax=[]; imin=[]; bad=(d==0); dd=diff([0bad0]); debs=find(dd==1); fins=find(dd==-1); ifdebs (1)==1 iflength(debs)>1 debs=debs(2: end); fins=fins(2: end); else debs=[]; fins=[]; end end iflength(debs)>0 iffins(end)==m iflength(debs)>1 debs=debs(1: (end-1)); fins=fins(1: (end-1)); else debs=[]; fins=[]; end end end lc=length(debs); iflc>0 fork=1: lc ifd(debs(k)-1)>0 ifd(fins(k))<0 imax=[imaxround((fins(k)+debs(k))/2)]; end else ifd(fins(k))>0 imin=[iminround((fins(k)+debs(k))/2)]; end end end end iflength(imax)>0 indmax=sort([indmaximax]); end iflength(imin)>0 indmin=sort([indminimin]); end end end end %%%%%%%%%%%%%%%%%%pos%%%%%%%%%%%% functionposs=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); indmin=find(d1.*d2<0&d1<0)+1; indmax=find(d1.*d2<0&d1>0)+1; ifany(d==0) imax=[]; imin=[]; bad=(d==0); dd=diff([0bad0]); debs=find(dd==1); fins=find(dd==-1); ifdebs (1)==1 iflength(debs)>1 debs=debs(2: end); fins=fins(2: end); else debs=[]; fins=[]; end end iflength(debs)>0 iffins(end)==m iflength(debs)>1 debs=debs(1: (end-1)); fins=fins(1: (end-1)); else debs=[]; fins=[]; end end end lc=length(debs); iflc>0 fork=1: lc ifd(debs(k)-1)>0 ifd(fins(k))<0 imax=[imaxround((fins(k)+debs(k))/2)]; end else ifd(fins(k))>0 imin=[iminround((fins(k)+debs(k))/2)]; end end end end iflength(imax)>0 indmax=sort([indmaximax]); end iflength(imin)>0 indmin=sort([indminimin]); end end end %说明每个程序要单独保存成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局域均值分解算法,目前个人没有时间研究了,希望得到指点。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- matlab 程序设计 源程序
![提示](https://static.bdocx.com/images/bang_tan.gif)