均值文件转matlab,LMD局域均值分解的matlab程序及示例代码

文件夹包含，找出纯调频函数，计算瞬时幅值，瞬时频率的函数等

%%%%%%%%%%%主程序%%%%%%%%%%

%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_cn_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&