线谱法 时钟分量的提取 matlab,LMD局域均值分解的matlab程序及示例

说明：研究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,LMD局域均值分解的matlab程序及示例的全部内容，希望文章能够帮你解决所遇到的问题。

如果觉得生活随笔网站内容还不错，欢迎将生活随笔推荐给好友。