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_cn_l
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局域均值分解算法，目前个人没有时间研究了，希望得到指点。
升级会员 