局部均值分解（matlab实现有注释）

       局部均值分解（local mean decomposition , LMD）方法同经验模态分解方法（EMD）一样，也是一种自适应信号处理方法。LMD通过改变信号分解过程能有效改进EMD方法存在的包络拟合不准确、边界处发散等问题.

      代码如下：

clear all;
clc;
% 这里是仿真信号，可改成自己的信号
fs=2000;
N=2048;
n=0:N-1;
t=n/fs;
x=(1+0.5*t).*sin(2*pi.*20*t)+2*cos(2*pi.*3*t);
%绘制仿真信号和其频谱图
figure(1)
subplot(211)
plot(t,x)
subplot(212)
y2=x;
L=length(y2);
NFFT = 2^nextpow2(L);
Y = fft(y2,NFFT)/L;
f = fs/2*linspace(0,1,NFFT/2);
plot(f,2*abs(Y(1:NFFT/2)))
% 局域均值分析
x=x';
c = x';
N = length(x);

A = ones(1,N);
PF = [];
aii = 2*A;

while(1)

  si = c;
  a = 1;
  
   while(1)
    h = si;
    
      maxVec = [];
      minVec = [];
      
   % 寻找极大值点和极小值点
      for i = 2: N - 1
         if h (i - 1) < h (i) & h (i) > h (i + 1)
            maxVec = [maxVec i]; 		
         end
         if h (i - 1) > h (i) & h (i) < h (i + 1)
            minVec = [minVec i]; 		
         end         
      end
      
   % 检查是否有残余
      if (length (maxVec) + length (minVec)) < 2
         break;
      end
           
  % 原始信号中的两边两个点的判断 
      lenmax=length(maxVec);
      lenmin=length(minVec);
      %先是左边这个点
      if h(1)>0
          if(maxVec(1)<minVec(1))
              yleft_max=h(maxVec(1));
              yleft_min=-h(1);
          else
              yleft_max=h(1);
              yleft_min=h(minVec(1));
          end
      else
          if (maxVec(1)<minVec(1))
              yleft_max=h(maxVec(1));
              yleft_min=h(1);
          else
              yleft_max=-h(1);
              yleft_min=h(minVec(1));
          end
      end
      %然后判断右边这个点
      if h(N)>0
          if(maxVec(lenmax)<minVec(lenmin))
             yright_max=h(N);
             yright_min=h(minVec(lenmin));
          else
              yright_max=h(maxVec(lenmax));
              yright_min=-h(N);
          end
      else
          if(maxVec(lenmax)<minVec(lenmin))
              yright_max=-h(N);
              yright_min=h(minVec(lenmin));
          else
              yright_max=h(maxVec(lenmax));
              yright_min=h(N);
          end
      end
      %使用三次样条插值方法，对极大值向量和极小值向量进行插值
      %spline interpolate
      maxEnv=spline([1 maxVec N],[yleft_max h(maxVec) yright_max],1:N);
      minEnv=spline([1 minVec N],[yleft_min h(minVec) yright_min],1:N);
      
    mm = (maxEnv + minEnv)/2;%得到局部均值函数
    aa = abs(maxEnv - minEnv)/2;%得到包络函数
    
    mmm = mm;
    aaa = aa;

    preh = h;
    h = h-mmm;%从原始信号中分离处局部均值函数
    si = h./aaa;%对分离出的信号进行解调
    a = a.*aaa;    
    
aii = aaa;

    B = length(aii);
    C = ones(1,B);
    bb = norm(aii-C);%返回aii-C的最大奇异值，aii就是那个包络函数
    if(bb < 1000)%如果bb<1000，就得到了纯调频函数
        break;
    end     
    
   end %分解1个Pf分量在这结束
   
  pf = a.*si;%包络函数和纯调频函数相乘，得到PF分量
  
  PF = [PF; pf];
  
  bbb = length (maxVec) + length (minVec);
 % 简单的一个结束分解的条件
      if (length (maxVec) + length (minVec)) < 20
         break;
      end
           
  c = c-pf;

end
m=x'-PF(1,:)-PF(2,:);%如果分解出2个，从原始信号中减去，得到残余分量

figure(2);
subplot(4,1,1),plot(n,x),ylabel('X(t)');
subplot(4,1,2),plot(n,PF(1,:)),ylabel('PF_1(t)');
subplot(4,1,3),plot(n,PF(2,:)),ylabel('PF_2(t)');
subplot(4,1,4),plot(n,m),ylabel('u(t)');

xlabel('Time / S'); 

最后分解出来的图：