文件夹包含,找出纯调频函数,计算瞬时幅值,瞬时频率的函数等
%%%%%%%%%%%主程序%%%%%%%%%%
%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&