类EMD的“信号分解方法”及MATLAB实现（第六篇）——LMD

继续完善“类EMD”方法系列，本篇是继EEMD、CEEMD、CEEMDAN、VMD、ICEEMDAN后的第6篇，想要看前几种方法的点击链接可以跳转。

LMD（local mean decomposition，局部均值分解）方法是2005年由Smith等人[1]提出的，本质上是根据信号的包络特征自适应地将一个非线性、非平稳信号按频率递减的顺序逐级分离。LMD的提出也是用来解决EMD分解的端点效应和模态混叠问题，最开始是用来处理脑电数据的。

1. LMD的概念

与EMD衍生的一系列方法不同，经过LMD分解得到的分量被称作“乘积函数(PF)”，即每个PF都是通过包络函数乘以纯调频函数得到的。其中包络函数是PF的瞬时幅值，纯调频函数的频率是PF的瞬时频率。

曲线为待分解曲线

图解一下分解步骤：

（1）找到信号中的全部极值点  。

（2）求出相邻极值点的局部均值点 ，也就是相邻两个极值点的中点。公式为： 

蓝色标记为局部均值点

（3）求出相邻极值点的局部包络  ，也就是相邻两个极值点的幅值绝对值差值的一半。公式为 

绿色标记为局部包络值

（4）将局部均值的折线连线进行平滑处理，得到局部均值函数 ，如下图中的橙色曲线。

橙色曲线为局部均值函数（示意）

（5）将局部包络同样进行平滑处理，得到局部包络估计函数  ，如下图紫色曲线。

局部包络估计函数（示意）

（6）将从原始信号的序列里面分割出来，得到零均值信号。

（7）对  进行解调，通过  除以  可以得到  。反复重复以上（1）~（6）步，直到包络估计函数  ，此时得到的  就是纯调频信号。

（8）将上述迭代获得的全部局部包络估计函数进行乘法运算就可以得到包络信号  。

（9）第一个PF分量就可以写为：

 （10）将原始信号减掉PF1，得到的剩余信号重复（1）~（9）步，直到剩余信号为单调函数为止，此时原始信号就被分解为k个PF和一个单调剩余信号。即：

分解完成！

2. LMD的特点

对于LMD和EMD的区别，这里就直接引用论文了

（1）PF 分量和 IMF 分量的含义不相同。经过 EMD 分解后获得的 IMF 属于调频信号，而利用 LMD 分解后获得的 PF 分量属于调幅调频信号。要想得到 IMF 分量，必须保证原始信号极值点的数量绝对等于过零点的数量，或者两个数值的数值差的结果小于等于 1，因此 IMF 分量并不会显现出不过零点的局部波动；但 PF并不需要满足这个条件。综合上所述可以发现 PF 分量能够更准确的反映原始信号的所有特征信息。
（2）针对局部均值函数 EMD 和 LMD 的求解方法存在明显的差异。EMD 是分别对所有极值点采用三次样条插值获得原始信号的上包络线和下包络线，接下来使用平均值的方法就可以得到局部均值函数，采用这种方法更容易形成过包络或欠包络等缺陷；对于LMD 求解局部均值函数，求取相邻两个极值的平均值，并利用滑动平均算法对其进行平滑处理；LMD 能够避免 EMD 中存在的过包络和欠包络的缺点。因此通过对比可以发现LMD 的分解结果更加准确。
（3）LMD 和 EMD 对瞬时频率的求解思路不同。在 EMD 中，必须求解 Hilbert 才能获得 IMF 的瞬时频率，然后再利用对其瞬时相位的求取倒数，最终获得瞬时频率，当若干个 IMF 中的一个瞬时相位具有突变的时候，求解出来的瞬时频率可能会出现工程实际中难以解释的负值；但对于 LMD 则不会出现频率为负值的情况，因为瞬时频率的值可以直接通过分解后的 PF 分量直接计算得到，这种方法不但简单而且求出的瞬时频率值都属于正值。因此求解瞬时频率的时候 LMD 方法更占优势。
（4）LMD 和 EMD 的整个分解过程计算量有所不同。针对 EMD 的求解过程主要存在两个迭代过程，一个是获取所需要的若干个 IMF 分量，另一个则是将所有的 IMF 从原始信号中分离出来，最终得到一个残余分量；而 LMD 方法就相对于 EMD 复杂一些，其主要包含了三个迭代过程，第一个迭代是利用滑动平均算法对局部均值和局部包络的折线进行平滑处理，得到局部均值函数以及包络估计函数，第二个迭代是通过迭代过程获得一个纯调频函数，第三个迭代就是将所有的 PF 分量全部计算出来。针对计算量来说，LMD的并不占优势，其计算量相比于 EMD 稍大一些。

3. LMD的MATLAB编程实现

在互联网上没有找到十分靠谱的LMD程序，于是笔者决定动手写一个出来。

但是写出来的代码运行却很有问题。

模态混叠端点效应都算小问题了。

PF1分量中存在模态混叠

某些信号分量会出现巨大的局部畸形，且其数量级远原始信号。

比如对上图中的原始信号加入微量噪声后，分解结果就会变成这样：

注：res还没达到单调，如果要单调需要分解到十几个分量，这里强行停止并将最后几个分量重构了

前三个PF算是勉强还能看的话，从PF4分量就开始暴走了。

所以目前的代码状态是不可用。但是处于交流的目的，代码还是贴上来，如果有同学解决了这个问题记得告诉我一声哦。

function pf = kLMD(data,tol,MaxNum)
% LMD分解函数
% 输入：
% data：待分解信号
% tol：迭代停止阈值
% MaxNum：最大PF分量数，设置此数是为了防止分解难以达到停止条件，陷入过多次循环。默认为10
% 输出：
% pf：乘积函数，即分解得到的分量，每一行为一个分量

% 鉴于互联网上没找到特别靠谱的代码，按照方法流程重头写了这个版本。
% 参考论文（即方法提出论文）为：
% Smith, Jonathan S . The local mean decomposition and its application to EEG perception data[J].
% Journal of The Royal Society Interface, 2005, 2(5):443-454.
% 笔者对于该方法的介绍文章（类EMD的“信号分解方法”及MATLAB实现（第六篇）——LMD）为：
% https://zhuanlan.zhihu.com/p/444277130
% 
%  Copyright (c) 2021 Mr.看海 All rights reserved.

if ~exist('MaxNum')  %如果未设置最大pf数，则设置为10
    MaxNum = 10;
end
k = 0;
while 1
    dataTemp{1} = data;

    k = k+1;
    [pf(k,:)] = getpf(dataTemp{k},tol,'off');
    if k == 1
        dataTemp{k+1} = dataTemp{k}(:)-(pf(:));    %将PF从信号中分割出来
    else
        dataTemp{k+1} = dataTemp{k}(:)'-(pf(k,:)); %将PF从信号中分割出来
    end
    
    if sum(diff(dataTemp{k+1})>0)==0 || sum(diff(dataTemp{k+1})<0)==0 %如果残值单调，结束分解
        pf(k+1,:) = dataTemp{k+1};  %残值
        break
    end
    if k>=MaxNum  %达到最大迭代次数，停止分解
        pf(k+1,:) = dataTemp{k+1};  %残值
        break
    end
end


end

function [pf] = getpf(data,tol,figflag)
% 通过迭代过程获得一个纯调频函数
i = 0;
while 1
    i = i + 1;
    if i == 1
        
        [m11,a11{i},h11,s11{i}] = getamhs(data,figflag);     %调用函数，求得m11，a11，s11等
    else 
        [m11,a11{i},h11,s11{i}] = getamhs(s11{i-1},figflag); %调用函数，求得m11，a11，s11等
    end
    
    if sum((a11{i} - 1).^2) < tol %达到阈值，迭代停止
        break
    end
end

a = 1;
for j = 1:i
    a = a.*a11{j};  %包络信号
end
pf = a.*s11{i}; %得到此阶的乘积函数
pf = pf';  %转变为行向量
end

function [mm11,aa11,h11,s11] = getamhs(data,figflag)
data = data(:)';
len = length(data);
%% 1.求极值点
[pksP,locsP] = findpeaks(data);  %找到最大值
[pksN,locsN] = findpeaks(-data); pksN = -pksN; %找到最小值
%% 2.求局部均值点m
ns = [locsP,locsN;pksP,pksN]; %极值点
[B,I] = sort(ns(1,:)); %获取排序索引
nsSort = ns(:,I);  %对极值点m按顺序排列
nsSortT = [[1;data(1)],nsSort,[len;data(len)]]; %将端点作为一个极值点加入
for i = 1:length(nsSortT(1,:))-1
    ms(i) = (nsSortT(2,i) + nsSortT(2,i+1))/2;  
    msLocs(i)  = (nsSortT(1,i) + nsSortT(1,i+1))/2;  
    ai(i) = abs(nsSortT(2,i) - nsSortT(2,i+1))/2;
    msLocs(i)  = (nsSortT(1,i) + nsSortT(1,i+1))/2;  
end
%% 3.求局部均值函数
m11 = movmean([msLocs',ms'],2); %滑动平均
try
    [fm,~] = fit(m11(:,1),m11(:,2),'smoothingspline'); %曲线拟合
catch ME
    m11 = movmean([msLocs',ms'],2); %滑动平均
    [fm,~] = fit(m11(:,1),m11(:,2),'smoothingspline'); %曲线拟合
end
mm11 = fm(1:len);
%% 4.局部包络
a11 = movmean([msLocs',ai'],2); %滑动平均
try
    [fa,~] = fit(a11(:,1),a11(:,2),'smoothingspline'); %曲线拟合
catch ME
    a11 = movmean([msLocs',ai'],2); %滑动平均
    [fa,~] = fit(a11(:,1),a11(:,2),'smoothingspline'); %曲线拟合
end
%[fa,~] = fit(msLocs',ai','smoothingspline'); %曲线拟合
aa11 = fa(1:len);
%% 5.求零均值信号
h11 = data-mm11'; %零均值信号
%% 6.调频信号
s11 = h11(:)./aa11(:);
end

测试用的仿真信号为：

fs = 1e3;
t = 0:1/fs:1-1/fs;
x = 8*t.^2 + cos(4*pi*t+10*pi*t.^2) + ...
    [cos(60*pi*(t(t<=0.5))) cos(200*pi*(t(t>0.5)-10*pi))];
n = rand(1,1000);
sig = x+0.1*n;

其他

EMD、EEMD、CEEMD、CEEMDAN、ICEEMDAN、VMD以及HHT相关的程序笔者封装了画图函数。关于EMD、EEMD、CEEMD、VMD和HHT的相关介绍可以看这里：

Mr.看海：这篇文章能让你明白经验模态分解（EMD）——EMD在MATLAB中的实现方法

Mr.看海：希尔伯特谱、边际谱、包络谱、瞬时频率/幅值/相位——Hilbert分析衍生方法及MATLAB实现

Mr.看海：类EMD的“信号分解方法”及MATLAB实现（第一篇）——EEMD

Mr.看海：类EMD的“信号分解方法”及MATLAB实现（第二篇）——CEEMD

Mr.看海：类EMD的“信号分解方法”及MATLAB实现（第三篇）——CEEMDAN

Mr.看海：类EMD的“信号分解方法”及MATLAB实现（第四篇）——VMD

Mr.看海：类EMD的“信号分解方法”及MATLAB实现（第五篇）——ICEEMDAN

参考

1. ^[1]Smith, Jonathan S . The local mean decomposition and its application to EEG perception data[J]. Journal of The Royal Society Interface, 2005, 2(5):443-454.