资讯详情

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

继续完善类EMD方法系列，本文是继承EEMD、CEEMD、CEEMDAN、VMD、ICEEMDAN在接下来的第六篇文章中，如果你想看到以前的方法，可以通过点击链接跳转。

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)步骤，直到剩余信号为单调函数，原始信号分解为kPF单调剩余信号。

分解完成！

LMD

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

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

LMD

没有在互联网上找到非常可靠的LMD所以作者决定写一个程序。

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

模态混叠端点效应是个小问题。

PF模态混合在1分量中存在

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

例如，将微量噪声添加到上图中的原始信号中，分解结果就会变成这样：

注：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         dataTep{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.

标签： m11传感器传感器imf18

锐单商城拥有海量元器件数据手册、 IC替代型号，打造 电子元器件IC百科大全！

锐单商城 - 一站式电子元器件采购平台