[转载]LMD经验模态分解matlab程序——原味的

sovereign

曾经也用滑动平均写过LMD，其实滑动平均的EMD才是原汁原味的居于均值分解。

分享给有需要的人，程序写的不好，只是希望提供一种思路。如果谁写了更完美LMD程序，别忘了发我一份，快毕业了，一直没有把LMD写完美，对于我来说始终是个遗憾。来分完美的LMD让我也品尝下，我也无憾了~

代码下载地址：http://download.csdn.net/source/3102096

此处没有提供测试代码，如需要可以点这里：点我

源代码如下：

%原始lmd算法，效果很不好，不知道程序哪里写错
function[PF,A,SI]=lmd(m)
c=m;
k=0
wucha1=0.001;
n_l=nengliang(m);
while 1
    k=k+1;
    a=1;
    h=c;
    [pf,a,si]=zhaochun(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);
    if length(c_pos)<3 || n_c<n_l/100 || length(pos(pf))<length(c_pos) || n_pf<n_c
        PF(k+1,:)=c;
        break
    end
end

function pos=pos(y)
%功能：找序列极值点位置坐标

%y:输入序列
%pos:极值点的序列位置坐标
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

minmax=cat(2,indmin,indmax);
pos=sort(minmax);
end
%----------zhaochun.m
function [pf,a,si]=zhaochun(a,h,wucha1)
chun_num=0;
while 1
chun_num=chun_num+1
t=1:length(h);
h_pos=position(h);%极值点位置序列
tpoint=t(h_pos);%极值点时间值
hpoint=h(h_pos);%极值点幅度值
hpoint=bianjie(h,hpoint,1);%端点处理后的极值点，多出了2个极值点
[mi,ai]=find_miai(hpoint);%找出极值点之间的均值函数和包络函数
mi_point=junzhi(mi);%所有极值点处的均值序列（幅值）-纵坐标（点序列）
ai_point=junzhi(ai);%所有极值点处的包络序列（幅值）-纵坐标（点序列）
%此处端点处的均值和包络都是 端点和极值点之间的均值和包络值（如果端点视作极值点，这样处理是不合理的，端点都不是极值点，这样处理是正确的）
lmi_point=mi(1);%左端点的均值
rmi_point=mi(length(mi));%右端点的均值
lai_point=ai(1);%左端点的包络
rai_point=ai(length(ai));%右端点的包络
%
mi_point_d=link(lmi_point,rmi_point,mi_point);%连接端点均值及所有极值点处的 均值 （带端点的均值序列）（点序列）
ai_point_d=link(lai_point,rai_point,ai_point);%连接端点包络及所有极值点出的 包络值（带端点的包络序列）（点序列）
%
tpoint_d=link(t(1),t(length(t)),tpoint);
mi_fun=chadian1(tpoint_d,mi,mi_point_d);%包含端点和极值点和普通点的 均值序列-平缓前的均值序列
ai_fun=chadian1(tpoint_d,ai,ai_point_d);%包含端点和极值点和普通点的 包络序列-平滑前的包络序列
%以上是完整的未处理的均值函数mi_fun和包络函数ai_fun

%找出第一平滑的滑动跨度
kmax=max(diff(tpoint_d));%找出时间跨度最大的  相邻几点 间的 距离
kmax1=uint16(kmax/3);
kmax1=double(kmax1);
jiou=uint8(rem(kmax1,2));
if kmax1<3
    move_k=3; % b<3,滑动跨度=3，
elseif jiou==0 % b>=3,当b是偶数时，跨度=b-1
            move_k=(kmax1-1);
   else move_k=kmax1; %b>=3,b是奇数时，跨度=b
end
[mi_move,move_mi_num]=move(mi_fun,move_k);
[ai_move,move_ai_num]=move(ai_fun,move_k);

mi=mi_move;
ai=ai_move;
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;
end
function [x,flag]=move(a,k)
l=length(a);
t=1:l;
% jingdu=t(2)-t(1);
flag=0;
x=a;
max_flag=10;%平滑重复的最大次数设置
max_error=0;%相邻两点间相等允许的最大差异（理论上=0才人为无差异，实际操作要考虑采样精度，所以可以设置一个允许范围）
while flag<=max_flag || min(abs(diff(x)))<=max_error;%这里用到abs，然后再min，因为幅值的差值会出现负值，我们的目的是找差值=0，而不是负数
    if flag==0
        flag=flag+1;%flag标志位，标志平滑的次数，这里这里设置为不超过11次（最大10次）
    else
        flag=flag+1;%flag标志位，标志平滑的次数，这里这里设置为不超过11次（最大10次）
        x_pos=position(x);%极值点位置序列
        tpoint=t(x_pos);%极值点时间值
        tpoint_d=link(t(1),t(l),tpoint);%极值点的时间轴上的位置
        kmax=max(diff(tpoint_d));%找出时间跨度最大的  相邻几点 间的 距离
        kmax1=uint16(kmax/3);
        kmax1=double(kmax1);
        jiou=uint8(rem(kmax1,2));
        if kmax1<3
            k=3; % b<3,滑动跨度=3，
        elseif jiou==0 % b>=3,当b是偶数时，跨度=b-1
                    k=(kmax1-1);
               else k=kmax1; %b>=3,b是奇数时，跨度=b
        end
    end
    y=zeros(1,l);% 初始化序列y,   y是中间变量
%%%%%%%%%%%%%%%%%%%边界点的处理%%%%%%%%%%%%%%
    for i=1:(k-1)/2
        v=0;
        z=0;
        for j=1:i*2-1
            v=v+x(j);
        end
        y(i)=v/(i*2-1);
        for j=l:-1:l+2-i*2           %j=l:-1:l+1-(i*2-1)
            z=z+x(j);
        end
        y(l+1-i)=z/(i*2-1);
    end
%%%%%%%%%%%%%%中间点的处理 %%%%%%%%%%%%%%%%%%%%%%%%%%
    for i=1+(k-1)/2:l-(k-1)/2 %这个 (d+1)/2是跨度的中点 单边的边界点数=中点值-1=(d+1)/2-1=(d-1)/2
        w=x(i);
        for j=1:(k-1)/2
            w=w+x(i-j)+x(i+j);
        end
        y(i)=w/k;
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                 
    x=y;
end
%k   %  调试的时候查看
%flag %  调试的时候查看
end

function hpp=bianjie(hh,hp,style)
%hh:需要边界处理的序列
%hp:需要边界处理序列的极值点的幅度值
%style:边界处理类型1：端点全部看作极值点2：左右两端各增加1个幅值为0的极值点
%hpp:边界处理后，极值点幅值多处了两个，因为左边右边各人为加上了一个
if style==1
    a=hh(1);
    b=hh(length(hh));
end
if style==2
    a=0;
    b=0;
end
hpp=link(a,b,hp);
end

function pos=position(y)
%功能：找序列极值点位置坐标

%y:输入序列
%pos:极值点的序列位置坐标
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

minmax=cat(2,indmin,indmax);
pos=sort(minmax);
end

function [mi,ai]=find_miai(ypoint)
%
Lyp=length(ypoint);
al=wkeep(ypoint,Lyp-1,'l');
ar=wkeep(ypoint,Lyp-1,'r');
mi=(al+ar)/2;
ai=abs(ar-al)/2;
end

function c=junzhi(a)
l=length(a);
al=wkeep(a,l-1,'l');
ar=wkeep(a,l-1,'r');
c=(al+ar)/2;
end

function d=link(a,b,c)
d=[a';c';b']';
%头：尾：中间
end

function f_value=chadian1(a,b,c)
% chadian1 把端点及极值点处对应的总坐标值插入，原来的均值函数的方波序列中
%输入参数a:点序列（行向量）（包含端点和极值点以 在时间上的位置-横坐标）（n个）（点序列-横坐标）
%输入参数b:段序列（行向量）（点序列a 每两点之间的纵坐标的值-纵坐标） （n-1）-（段序列-纵坐标）
%输入参数c:点序列（行向量） （包含端点和极值点 在对应时间上的幅值-纵坐标）（n个）-（点序列-纵坐标）
%输入参数d:原始数据的采样精度
%输出参数f_value（行向量）: 点序列a 插入 段序列的值 之后，以c的精度的 值（对应于 横坐标，纵坐标的值）
%精度是0.001
l=length(b);
al=wkeep(a,l,'l');
ar=wkeep(a,l,'r');
d={[]};%d={0}这样是为了初始化 元胞数组
for i=1:l                                 %采样精度0.001
    d{i}=ones(1,(ar(i)-al(i)-1))*b(i);%ones函数参数要为整数，unint16就是数据强制类型转换，
end                                                           %这里没有使用到单独为uint16((ar(i)-al(i))*1000)-1)这个自变量赋值，所以只是个中间变量，对数据不会产生污染
y=c(1);
for i=1:l
    y=link(y,c(i+1),d{i});
end
f_value=y;
end
end
%------
function p=nengliang(y)
% my=mean(y);
% p=(y-my).*(y-my);
% p=sum(p);
p=sum(abs(y).^2);
end
%--------

转自：http://blog.sina.com.cn/s/blog_3f87c5820100vout.html