CN102697495B - 基于总体平均经验模式分解的二代小波肌电信号消噪方法 - Google Patents

Abstract

本发明涉及一种基于总体平均经验模式分解的二代小波肌电信号消噪方法。本发明采用基于总体平均经验模式分解的二代小波消噪方法，获取人体的肌电信号样本数据，对原始的肌电信号加入白噪声后进行经验模态分解，然后对高频的内蕴模式函数分量进行二代小波分解及阈值处理，小波重构高频分量。最后把经过处理的高频分量与低频分量进行叠加，重构后的信号即为去噪信号。本发明将信号自适应地分解到不同的尺度上，适合非线性、非平稳信号的处理，不仅具有小波分析的全部优点，还有更清晰准确的谱结构，并能改善信号的极值点分布，具有抗混分解能力，既尽可能地保存了有用信号，又有效地消除了噪声，并且能够大幅提高信号信噪比。

Description

基于总体平均经验模式分解的二代小波肌电信号消噪方法

技术领域

本发明属于信号消噪领域，涉及一种基于总体平均经验模式分解的二代小波肌电信号消噪方法。

背景技术

肌电信号（EMG）是人体自主运动过程中神经肌肉活动发出的生物电信号在体表的综合反映，已在临床、运动医学等领域被广泛的研究。同时由于其与肌肉的活动状态和功能状态之间存在着不同程度的关联性，在一定程度上反映神经肌肉的活动，因而成为了人工假肢和功能性神经电刺激的理想控制信号。

但是，由于肌电信号通常是在皮肤表面募集的，是一种非线性非平稳信号,有用能量主要分布在0-500 Hz之间，非常微弱,因而不可避免地会夹杂着生理噪声。同时，在拾取及处理过程中，存在外界或电路内部因素的干扰，使得实际采集的肌电信号是被污染的信号，如果直接用其作为控制源对假肢实施动作控制，就会造成不可挽回的损失。因而探索一种能将夹杂在采集信号中噪声的滤去方法是很有必要的。

目前已有的肌电信号消噪方法主要包括以下几种：

1、传统傅里叶变换滤波

肌电信号表现出明显的非线性、非平稳特性，给信号降噪和特征提取取带来了很大困难。传统傅里叶变换是建立在信号平稳性假设条件下的一种时域和频域的全局性变换，对于分析基于平稳信号或准平稳过程的信号十分有效，但其不能很好地揭示非平稳信号的特征。

2、小波变换

小波变换是近年来出现的新的时频分析方法，在生物医学信号处理中已得到了广泛的应用。小波消噪的方法大概可以分为三大类：第一类方法是基于小波变换模极大值原理的，最初由Mallat提出，即根据信号和噪声在小波变换各个尺度上的不同传播特性，剔除由噪声产生的模极大值点，保留信号所对应的模极大值点，然后利用所余模极大值点重构小波系数，进而恢复信号；第二类方法是对含噪信号作小波变换之后，计算相邻尺度间小波系数的相关性，根据相关性的大小区分小波系数的类型，进行取舍，然后直接重构信号；第三类方法是Donoho提出的阈值方法，该方法认为信号的小波系数是一致分布的，个数较多，但幅值小。基于这一思想，Donoho等人提出了软阈值和硬阈值去噪方法，即在众多小波系数中，把绝对值较小的系数置为零，而让绝对值较大的系数保留或收缩，分别对应于硬阈值和软阈值方法，得到估计小波系数，然后利用估计小波系数直接进行信号重构，即可达到去噪的目的。

小波变换具有良好的时频局部化特性，使其在肌电信号处理中得到广泛应用。但小波分析本质上是一种窗口可调的傅里叶变换，存在小波基选择困难，参数敏感和平稳性假设等问题，且一旦小波基被选定之后，就必须用小波基分析整个信号，缺乏自适应性。

3、经验模态分解（EmpiricalMode Decomposition，EMD）方法则是一种新的具有自适应性的时频分析方法，可根据信号的局部时变特性进行自适应的时频分解，比较适合对非线性、非平稳的肌电信号进行处理，但其重要缺陷就是产生模式混淆现象。

综上所述，已有的肌电信号消噪方法都存在不足之处。

发明内容

本发明的目的就是针对现有技术的不足，提出一种基于总体平均经验模式分解(Ensemble empirical mode decomposition ,EEMD)的二代小波变换肌电信号消噪新方法。总体平均经验模式分解方法本质上是一种借助于噪声的信号处理方法，由于在不连续的信号中，高频成分的极值点间隔分布密集，低频成分的极值点间隔分布稀疏，通过将白噪声加入到被分析的信号中，能够改变信号低频成分的极值点分布特性，使信号在整个频带中极值点间隔分布均匀，也保证了每次对信号进行准确求取上、下包络的局部均值，从而达到在经验模式分解中避免模式混淆的目的。

首先对信号加入白噪声处理后进行经验模态分解（EMD），然后对高频的内蕴模式函数(IMF)分量进行二代小波分解及阈值消噪处理，把经过处理的高频IMF分量与低频IMF分量进行叠加，重构后的信号即为去噪信号。

为了实现以上目的，本发明方法主要包括以下步骤：

步骤一，在含噪原始肌电信号                                               

中多次加入具有均值为零、标准差为常数的白噪声

, 即 

式中，

原始肌电信号，信号长度记为

,

为第

次加入白噪声后的信号，为第

次加入的白噪声。

步骤二，对步骤一所得的信号

分别进行经验模态分解，得到各自的内蕴模式函数（IMF）分量，记为

，与一个余项，记为

。其中

表示第次加入白噪声后分解所得的第

个IMF。

步骤三，利用不相关随机序列的统计均值为0的原理，将步骤二所得的IMF分量进行总体平均运算，得到分解后最终的IMF，即

      

式中，N 为添加白噪声的次数，

为对原始肌电信号进行总体平均经验模态分解后所得的第

个IMF分量。当

越大，对应的白噪声的IMF的和将趋于0。此时总体经验模态分解的结果为

             

       

式中J为分解的最高层数，

为最终的残余分量。

步骤四，对处于高频的前K项IMF分量

进行累加记为信号

,即

=

。对信号进行二代小波分解，得到低频系数

和高频系数

,式中l为二代小波分解的层级，最高层数记为L。

对信号

进行二代小波分解的具体步骤如下：

将信号

分割成相互关联的两部分，即将一个序列分为偶数序列

和奇数序列

，且和

相关性越强，分裂效果越好。

预测：针对原信号数据间的相关性，可用

预测

，预测误差为，也即是高频系数：

，式中

为预测算子，为奇样本的估计值。当信号

的相关性较大时，预测效果将是非常有效的。

更新：通过上述两个步骤产生的偶数序列

在某些整体性质上会与原始数据有所区别，为了对

预测产生更好的数据，需要引入修正算子

进行更新。

步骤五，对步骤四所得到的高频系数

进行阈值消噪处理得到消噪后的高频系数

；

所述的阈值消噪处理采用如下公式：

式中

，其中M为信号

的长度，

，

为实验可调参数，只要选取合适的参数，就可以既保证信号的整体连续性，抑制振荡，又能减少

与

之间存在的偏差，提高信号的精度。

步骤六,将步骤五所得的

与步骤四所得的

进行小波重构得到消噪后的高频IMF分量累加信号

。

步骤七，用步骤六所得的

和后J-K项IMF分量来重构信号，得到消噪后的肌电信号。即

   

式中 

为最终消噪后的肌电信号。 

本发明与已有的诸多肌电信号的消噪算法相比，具有如下特点：

由于肌电信号是非线性、非平稳的随机信号，有用能量主要集中在信号的低频部分，利用本发明作消噪处理，它将信号自适应地分解到不同的尺度上，非常适合非线性、非平稳信号的处理，不仅具有小波分析的全部优点，还有更清晰准确的谱结构，并能改善信号的极值点分布，具有有效的抗混分解能力，既尽可能的保存了有用信号，又有效的消除了噪声，并且能够大幅提高信号的信噪比。

附图说明

图1为本发明的实施流程图；

图2为本发明实施例的原始含噪声的肌电信号图；

图3为本发明实施例经过EEMD分解的结果图；

图4为本发明实施例的处于高频的前3项IMF分量累加信号图；

图5为本发明实施例前3项IMF分量累加信号经二代小波消噪后的信号图；

图6为本发明实施例的消噪结果图。

具体实施方式

下面结合附图对本发明的实施例作详细说明：本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

如图1所示，本实施例包括如下步骤：

步骤一，在含噪原始肌电信号中多次加入具有均值为零、标准差为常数的白噪声

, 即

                          (1)

式中，

是本实施例的含噪声的原始肌电信号，如图2所示，信号长度

=2500,

为第

次加入白噪声后的信号，

为第

次加入的白噪声，添加的白噪声标准差为0.15，添加白噪声的总体次数N为200。

步骤二，对步骤一所得的信号分别进行经验模态分解，得到各自的内蕴模式函数（IMF）分量，记为，与一个余项，记为

。其中

表示第

次加入白噪声后分解所得的第

个IMF。

步骤三，利用不相关随机序列的统计均值为0的原理，将步骤二所得的IMF分量

进行总体平均运算，得到分解后最终的IMF，即

                   (2)

式中，N 为添加白噪声的次数,在本实施例中N=200，

为对原始肌电信号进行总体平均经验模态分解后所得的第

个IMF分量。此时总体经验模态分解的结果为

             

                    （3）

式中J为分解的最高层数,在本实施例中J=8，为最终的残余分量。

本实施例肌电信号进行总体平均经验模式分解的结果如图3所示。

步骤四，对处于高频的前K项IMF分量

进行累加记为信号

,即

= 。对信号

进行二代小波分解，得到低频系数

和高频系数

,式中l为二代小波分解的层级，最高层数记为L。在本实施例中K=3，L=3,信号

如图4所示。

对信号

进行二代小波分解的具体步骤如下：

将信号

分割成相互关联的两部分，即将一个序列分为偶数序列

和奇数序列

，且

和

相关性越强，分裂效果越好。

预测：针对原信号数据间的相关性，可用

预测

，预测误差为

 ，也即是高频系数：

，式中

为预测算子，

为奇样本的估计值。当信号

的相关性较大时，预测效果将是非常有效的。

更新：通过上述两个步骤产生的偶数序列

在某些整体性质上会与原始数据有所区别，为了对

预测产生更好的数据，需要引入修正算子

进行更新。

步骤五，对步骤四所得到的高频系数

进行阈值消噪处理得到消噪后的高频系数

；

所述的阈值消噪处理采用如下公式：

                   （4）

式中

，其中M为信号

的长度,

，

为实验可调参数，只要选取合适的参数，就可以既保证信号的整体连续性，抑制振荡，又能减少

与

之间存在的偏差，提高信号的精度。

在本实施例中M=2500，L=3,a=2。

步骤六,将步骤五所得的

与步骤四所得的进行小波重构得到消噪后的高频IMF分量累加信号

。本实施例前3项IMF分量累加信号经二代小波消噪后的信号

如图5所示。

步骤七，用步骤六所得的

和后J-K项IMF分量来重构信号，得到消噪后的肌电信号。即

   

                               （5）

式中 

为最终消噪后的肌电信号。 本实施例J=8，K=3，消噪后的肌电信号

如图6所示。

Claims (1)

1.基于总体平均经验模式分解的二代小波肌电信号消噪方法，其特征在于该方法包括如下步骤：

步骤(1).在含噪原始肌电信号x(m)中多次加入具有均值为零、标准差为常数的白噪声n_i(m),即：

x_i(m)=x(m)+n_i(m)

式中x(m)原始肌电信号，信号长度记为M,x_i(m)为第i次加入白噪声后的信号，n_i(m)为第i次加入的白噪声；

步骤(2).对步骤(1)所得的信号x_i(m)分别进行经验模态分解，得到各自的内蕴模式函数分量与一个余项，分别记为I_ij(m)和r(m)；其中I_ij(m)表示第i次加入白噪声后分解所得的第j个内蕴模式函数分量；

步骤(3).利用不相关随机序列的统计均值为0的原理，将步骤(2)所得的内蕴模式函数分量I_ij(m)进行总体平均运算，得到分解后最终的内蕴模式函数分量，即

$I_j\left(m\right)=\frac{1}{N}\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{I}_{\mathrm{ij}}\left(m\right)$

式中，N为添加白噪声的次数，I_j(m)为对原始肌电信号加入N次白噪声进行总体平均经验模态分解后所得的第j个内蕴模式函数分量；当N越大，对应的白噪声的内蕴模式函数分量的和将趋于0；此时总体经验模态分解的结果为：

$x\left(m\right)=\underset{j=1}{\overset{J}{\mathrm{\Σ}}}{I}_j\left(m\right)+r\left(m\right)$

式中J为分解的最高层数；

步骤(4).对处于高频的前K项内蕴模式函数分量I_j(m)进行累加记为信号y[m],即

对信号y[m]进行二代小波分解，得到低频系数c^l[m]和高频系数d^l[m],式中l为二代小波分解的层级，最高层数记为L；

步骤(5).对步骤(4)所得到的高频系数d^l[m]进行阈值消噪处理得到消噪后的高频系数d′^l[m]；

所述的阈值消噪处理采用如下公式：

$d^{\&prime;l}\left[m\right]=\left\{\begin{array}{cc}{d}^l\left[m\right]\&times;\frac{{\left|d^l\right[m\left]\right|}^a-{{\mathrm{\&lambda;}}_l}^a}{{\left|d^l\right[m\left]\right|}^a}& \left|d^l\right[m\left]\right|\&GreaterEqual;{\mathrm{\&lambda;}}_l\\ 0& \left|d^l\right[m\left]\right|<{\mathrm{\&lambda;}}_l\end{array}\right.$

式中 ${\mathrm{\&lambda;}}_l=\mathrm{\&sigma;}\sqrt{2{\log }_2\left(M\right)},$ 其中M为信号y[m]的长度， $\mathrm{\&sigma;}=\underset{l=\mathrm{1,2},...,M}{\mathrm{median}}\left(\right|d^L\left[m\right]\left|\right)/0.6745,$ a为实验可调参数；

步骤(6).将步骤(5)所得的d′ⁱ[m]与步骤四所得的cⁱ[m]进行小波重构得到消噪后的内蕴模式函数分量累加信号y′[m]；

步骤(7).用步骤(6)所得的y′[m]和后J-K项内蕴模式函数分量来重构信号，得到消噪后的肌电信号；即

$x^{\&prime;}\left(m\right)=y^{\&prime;}\left[m\right]+\underset{i=K+1}{\overset{J}{\mathrm{\&Sigma;}}}{I}_j\left(m\right)$

式中x'(m)为最终消噪后的肌电信号。

Images