第1页 / 共11页
第2页 / 共11页
第3页 / 共11页
第4页 / 共11页
第5页 / 共11页
第6页 / 共11页
第7页 / 共11页
第8页 / 共11页
LMD的方法改进.docx
说明：研究 LMD 局域均值分解有 3 个月左右，能找到的相关文章也基本上看了一遍，觉得是
个很好的方法，号称是 EMD 经验模态分解的改进版。但是网络上一直没有找到该算法的
matlab 程序，只见文章说的天花乱坠。后来自己写了一个，但是使用有一个问题没有解决，
就是分解的时候怎么去除骑行波的问题。先把自己写的程序贡献出来，让大家分析下，一起
讨论下，到底 LMD 的程序怎样写才能如文献中说的那样达到目的。欢迎大家热烈讨论！
程序仍存在不收敛的问题，拿出来分享只是希望高手指点一二，写的不好欢迎拍砖！
想要 M 文件的，给出下载地址：http://download.csdn.net/source/3090503
文件夹包含，找出纯调频函数，计算瞬时幅值，瞬时频率的函数等
%%%%%%%%%%%主程序%%%%%%%%%%
%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_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)
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
