说明：研究LMD局域均值分解有3个月左右，能找到的相关文章也基本上看了一遍，觉得是个很好的方法，号称是EMD经验模态分解的改进版。但是网络上一直没有找到该算法的matlab程序，只见文章说的天花乱坠。后来自己写了一个，但是使用有一个问题没有解决，就是分解的时候怎么去除骑行波的问题。先把自己写的程序贡献出来，让大家分析下，一起讨论下，到底LMD的程序怎样写才能如文献中说的那样达到目的。欢迎大家热烈讨论！
程序仍存在不收敛的问题，拿出来分享只是希望高手指点一二，写的不好欢迎拍砖！
文件夹包含，找出纯调频函数，计算瞬时幅值，瞬时频率的函数等
%%%%%%%%%%%主程序%%%%%%%%%%
%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);
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)的极值点个数(目前是用这个作为停止条件的一个参考写入程序)
%3.当前PF分量的能量应该大于残量c的能量(这个有待商榷)
%4.残余能量不能大于信号能量
if
n_pf
if
length(c_pos)<3 ||
n_c
n_pfn_l
PF(k+1,:)=c;
break
end
end
end
%%%%%%%%%%%%%%%%%%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)
chun_num=0;
while 1
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)
%computes envelopes and mean with various interpolations
NBSYM = 2; % 边界延拓点数
DEF_INTERP = 'spline';
if nargin < 2
x = t;
t = 1:length(x);
INTERP = DEF_INTERP;
end
if nargin == 2
if ischar(x)
INTERP = x;
x = t;
t = 1:length(x);
end
end
if ~ischar(INTERP)
error('interp parameter must be ''linear'''',
''cubic'' or ''spline''')
end
if ~any(strcmpi(INTERP,{'linear','cubic','spline'}))
error('interp parameter must be ''linear'''',
''cubic'' or ''spline''')
end
if min([size(x),size(t)]) > 1
error('x and t must be vectors')
end
s = size(x);
if s(1) > 1
x = x';
end
s = size(t);
if s(1) > 1
t = t';
end
if length(t) ~= length(x)
error('x and t must have the same length')
end
lx = length(x);
[indmin,indmax,indzer] = extr(x,t);
%boundary conditions for interpolation
[tmin,tmax,xmin,xmax] =
boundary_conditions(indmin,indmax,t,x,NBSYM);
% definition of envelopes from interpolation
envmax = interp1(tmax,xmax,t,INTERP); envmin = interp1(tmin,xmin,t,INTERP);
if nargout > 2
envmoy =
(envmax + envmin)/2;
end
end
function [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) < indmin(1)%
第一个极值点是极大值
if x(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
if x(1) <
x(indmax(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
%
序列末尾情况与序列开头类似
if indmax(end) < indmin(end)
if x(end) <
x(indmax(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
if x(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);
% 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)));
else
lmin =
fliplr(indmin(1:min(end,nbsym)));
end
if lsym == 1%
这种情况不应该出现，程序直接中止
error('bug')
end
lsym = 1; %
直接关于第一采样点取镜像
tlmin =
2*t(lsym)-t(lmin);
tlmax =
2*t(lsym)-t(lmax);
end %
序列末尾情况与序列开头类似
if trmin(end) < t(lx) | trmax(end)
< t(lx)
if rsym == indmax(end)
rmax =
fliplr(indmax(max(end-nbsym+1,1):end));
else
rmin =
fliplr(indmin(max(end-nbsym+1,1):end));
end
if rsym == 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 = [tlmin t(indmin) trmin];
tmax = [tlmax t(indmax) trmax];
xmin = [xlmin x(indmin) xrmin];
xmax = [xlmax x(indmax) xrmax];
end
%---------------------------------------------------------------------------------------------------
% 极值点和过零点位置提取
function [indmin, indmax, indzer] = extr(x,t);
%extracts the indices corresponding to extrema
if(nargin==1)
t=1:length(x);
end
m = length(x);
if nargout > 2
x1=x(1:m-1);
x2=x(2:m);
indzer = find(x1.*x2<0);
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)-1;
indz = round((debz+finz)/2);
else
indz = iz;
end
indzer = sort([indzer
indz]);
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;
% when two or more consecutive points have the same value we
consider only one extremum in the middle of the 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 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 end
end
%%%%%%%%%%%%%%%%%%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);
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
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局域均值分解算法，目前个人没有时间研究了，希望得到指点。
![matlab各个指令的含义,[MATLAB基础] 求解这段指令的意思，越详细越好，谢谢啦](https://www.labfans.com/bbs/t14414/images/icons/icon1.gif){kind=icon}
