线谱法 时钟分量的提取 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</