CN105997067A - 基于分数阶傅立叶变换的自适应肌电信号检测处理方法 - Google Patents

Abstract

本发明公开了一种基于分数阶傅立叶变换的自适应肌电信号检测处理方法。该方法首先对信号进行放大和滤去工频干扰等处理，利用智能分段，不断地将信号分段，同时进行自适应分数阶傅立叶变换滤波，在分数阶傅立叶域上将噪声信号滤除，再通过分数阶傅立叶反变换将信号还原至时域，最后求取肌电信号的特征值判断肌电信号所表示的肌肉状态。本发明采用分数阶傅立叶变换滤波算法分析肌电信号，更好地利用了肌电信号的时频信息，更好地分辨出了两种不同状态下的肌电信号的区别。

Description

基于分数阶傅立叶变换的自适应肌电信号检测处理方法

技术领域

本发明属于数字信号处理和医学领域，具体涉及一种基于分数阶傅立叶变换的自适应肌电信号检测处理方法。

背景技术

表面肌电信号是生物神经肌肉系统活动时所产生的生物电信号，它同生物神经肌肉的活动密切相关，其中包含着许多的肌体运动的信息。对表面肌电信号的监测和分析技术已经广泛地应用于临床医学、生物医学、运动医学等多个领域。由于肌电信号是肌肉神经的生理活动产生的，它具有典型非线性、低信噪比、非平稳的特征。

自1807年法国数学家傅立叶率先提出傅里叶分析的数学方法来分析信号频域信息以来，傅立叶变换已经得到了迅速而广泛的推广与应用。但是，由于傅立叶变换是一种全局性的信号的变换，所以难以去表述信号的局部信息，特别是信号的局部时频信息。而在一些非平稳信号的分析与应用中，此种信号的局部时频信息恰恰是其中至关重要、不容忽视的部分。

传统的傅立叶变换是一种线性算子，可看成是将信号从时间轴逆时针旋转π/2到频率轴。而分数阶傅立叶则是可旋转任意角度α的算子，从而可以得到信号在任意角度下的分析，而不是仅仅局限于频率轴上的信号分析方法。在时域上，依据Wigner-Ville分布对于分数阶傅立叶变换的旋转不变性，分数阶傅立叶可以看作是信号在时域平面绕原点旋转任意角度α后所构成的分数阶傅立叶域u上的表示，其与Wigner-Ville变换、短时傅立叶变换之间存在的只是坐标旋转变化关系，并不影响信号的时域分布特性。这个性质就可以在非平稳过程的情况下，利用分数阶傅立叶域来实现时频滤波，从而在减少计算开支的前提下也可以得到更好的滤波效果。

发明内容

本发明的目的在于针对现有技术的不足，提供一种基于分数阶傅立叶变换的自适应肌电信号检测处理方法。

本方法所采用的技术方案，包括以下步骤：

1)利用信号采集装置，对肌电信号进行采集获取，并对采集得到的原始信号进行预处理，得预处理后的肌电噪声混合信号S_SMG+N(t)。

2)将步骤1中采集得到的预处理后的肌电噪声混合信号S_SMG+N(t)，预设定智能分段的初始参数。

3)对信号S_EMG+N(t)进行智能分段。对分段后的小段信号S_i的时频信息F_i进行处理并线性拟合以求得该小段信号最优分数阶傅立叶变换的旋转角度α_iopt和最优分数阶傅立叶变换的阶数p_iopt。

4)对截取后的剩余信号再进行截取且令i＝i+1，重复步骤3，直至得到整段S_EMG+N(t)肌电信号所有信号段的最优阶数p_iopt的集合p_opt(t)。对信号S_EMG+N(t)进行自适应p_opt(t)阶分数阶傅立叶变换和u域窗函数滤波。

5)对步骤4所得去除噪声后的u域上的肌电信号S_EMG(u)进行p_opt(t)阶分数阶傅立叶反变换，得到时域上的去除噪声的肌电信号S_EMG(t)。

6)对步骤5得到的去除噪声后的肌电信号S_EMG(t)提取特征值，判断肌肉所处的状态。

本发明的有益效果在于：选用了分数阶傅立叶变换,相较于传统傅立叶变换，充分利用了输入肌电信号的局部时频信息。并且针对不同的输入信号，利用了拟合算法自适应地调节分数阶傅立叶的u域旋转角度α，从而更好得适应复杂的肌电信号输入。

附图说明

图1本发明整体实现过程框图；

图2预处理后两种不同状态下的肌电信号图；

图3原始肌电信号时频图；

图4分数阶傅立叶变换后肌电信号时频图；

图5分数阶傅立叶变换滤波算法后两种不同状态肌电信号图；

图6多个数据组滤波前后肌电信号特征值分布图。

具体实施方式

传统的肌电信号模型，可分成两部分，一部分是有用信号另外一部分是噪声信号e，即其数学模型如下：

对于该肌电信号在傅立叶域上的滤波，

由于肌电信号为非平稳过程，其滤波算子f一般不为时不变，也无法表示为时域卷积或传统傅立叶域上的乘法滤波器，利用时变算子f(t)对原始肌电信号进行滤波。为了得到时变算子f(t)，本实施例将经过放大处理后的待滤波信号进行短时拆分，并对其进行p阶分数阶傅立叶变换(FRFT)。分数阶傅立叶变换的定义式可以写成如下形式：

$f_p\left(u\right)=\frac{A_{\mathrm{\α}}}{\sqrt{2\mathrm{\π}}}{e}^{{\mathrm{jC}}_{\mathrm{\α}}{u}^2}{\∫}_{-\mathrm{\∞}}^{+\mathrm{\∞}}\[f\left(t\right)e^{{\mathrm{jC}}_{\mathrm{\α}}{t}^2}\]e^{-{\mathrm{jtuB}}_{\mathrm{\α}}}dt---\left(3\right)$

其中

在短时间内，将待滤波的肌电信号看做为平稳过程，其滤波算子f(t)采用时不变算子f。从而将时间不断递推，最终得到时变的滤波算子f(t)。本实施例是对肌电信号做p阶分数阶傅立叶变换，在分数阶傅立叶域(u域)对信号进行滤波，所以滤波算子为f(u)。

下面，将结合附图1对该处理肌电信号的方法进行解释。

1)利用信号采集装置，对肌电信号进行采集获取，并对采集得到的原始信号进行预处理，得预处理后的肌电噪声混合信号S_SMG+N(t)；

具体为：利用传感器和信号采集电路对肌电信号进行采集并将采集到的数据进行放大、简单滤波和模数转化等预处理，得到预处理后的肌电噪声混合信号。其中简单滤波为1.5倍增益中心频率为50Hz的巴特沃斯带阻滤波器，主要用以滤除50Hz的工频信号干扰。

2)将步骤1中采集得到的预处理后的肌电噪声混合信号S_SMG+N(t)，预设定智能分段的初始参数；

具体为：设智能分段第i段分割时间长度为t_i，先预设t_i初值为总信号长度L的十分之一，即t_i＝L/10。若求得t_i存在小数，则进行四舍五入求整。按时间t_i分割出第i小段信号S_i(t)，简写为S_i。并进行p＝0阶傅立叶变换，得到信号的时频关系F_i。

3)对信号S_EMG+N(t)进行智能分段。对分段后的小段信号S_i的时频信息F_i进行处理并线性拟合以求得该小段信号最优分数阶傅立叶变换的旋转角度α_iopt和最优分数阶傅立叶变换的阶数p_iopt；

具体为：定义瞬时信号绝对值|S_i|大于该小段信号绝对值的最大值|S_i|_max三分之一的肌电信号点称为该段信号的有效点，记为S_ief。

定义有效点数为该小段信号S_i的有效点S_ief个数，记为N_i。

求S_i信号长度L_i与有效信号点数N_i，若N_i大于或等于整段信号长度L_i的2/3，或L_i小于等于100，则将该小段信号S_i的时频关系F_i进行线性拟合，得到拟合斜率k_i，取α_i＝tan^-1(-k_i)，则根据公式3，分数阶傅立叶变换阶数p_i＝2α_i/π。反之，若该小段信号长度t_i为偶数则信号长度再截半，即t_i＝t_i/2，若t_i现为奇数，则信号长度加1再截半,即t_i＝(t_i+1)/2。将截取后的信号记为S_i重新进行本步骤，直至求出S_i的最优分数阶傅立叶变换的旋转角度α_iopt和最优分数阶傅立叶变换的阶数p_iopt。

4)对截取后的剩余信号再进行截取且令i＝i+1，重复步骤3，直至得到整段S_EMG+N(t)肌电信号所有信号段的最优阶数p_iopt的集合p_opt(t)。对信号S_EMG+N(t)进行自适应p_opt(t)阶分数阶傅立叶变换和u域窗函数滤波；

具体为：假设经过步骤3智能分段后将信号S_EMG+N(t)一共分为n段小段信号，分别为S_i，i＝1,2,3,...,n。

定义肌电信号最优阶数集合p_opt(t)为所有p_iopt按时间分布的分段函数，具体数学表示为p_opt(t_i)＝p_iopt，i＝1,2,3,...,n。

对S_EMG+N(t)做p_opt(t)阶傅立叶变换，对每段信号S_i在u域上加窗函数f_i(u)，使每小段有效信号S_i能在u域的窗函数内，在u域上得到经过滤波后的肌电信号S_EMG(u)。

5)对步骤4所得去除噪声后的u域上的肌电信号S_EMG(u)进行p_opt(t)阶分数阶傅立叶反变换，得到时域上的去除噪声的肌电信号S_EMG(t)；

6)对步骤5得到的去除噪声后的肌电信号S_EMG(t)提取特征值，判断肌肉所处的状态。

具体为：求去除噪声后的肌电信号的中值频率f_mid、与信号的标准差σ等特征值，得出不同肌电信号的动作差异，以识别肌肉所处状态。其具体数值公式分别为：

$f_{mid}=\frac{1}{2}{\∫}_0^{+\mathrm{\∞}}PSD\left(f\right)df---\left(4\right)$

$\mathrm{\σ}=\frac{1}{L_S-1}\underset{t=1}{\overset{L_s}{\mathrm{\Σ}}}{S}_{EMG}^2\left(t\right)---\left(5\right)$

其中，PSD(f)为肌电信号的功率谱密度函数，Ls为滤波后肌电信号的长度，S_EMG(t)为滤波后肌电信号的瞬时值，由于一般肌电信号的均值接近0，为了便于计算，这里将肌电信号的均值设为0。

以下通过实例对本发明内容进行进一步解释，用该方法对实际肌电信号进行检测处理。由于实际采集环境的工频噪声干扰等原因，淹没了肌电信号的特征成分，模糊了肌电信号的特征值，如图2所示。对于两种不同状态(正常和病变)的肌肉状态也较难区分(信号1长度为34.154s，信号2长度为38.769s，为了便于观察信号波形，这里都截取两信号第1秒时间内的信号做图)。本方法对原始肌电信号进行采集放大等初步处理后，采用自适应分数阶傅立叶变换滤波算法对初步处理后的肌电信号进行分析处理，这里以信号1为例，该段信号长度L为34154ms，采样频率SP为20kHz。故预设分割时间t₁＝3415ms。采用自适应拟合截取算法得到的分数阶傅立叶域时频特性如图3所示。自适应拟合得到p＝0.42，做p＝0.42的分数阶傅立叶变换后，信号时频特性如图4所示。最后，算法不断截取肌电信号得到p_iopt，最后得到整段的时变的分数阶傅立叶阶数p_opt(t)，滤波后在进行信号特征值分析，可以看出比滤波前有较大的区分度，见图5，这里采取了20个人体肱二头肌肌电信号样本(正常和病变分别为10个，每个人10个段肌电信号数据，每个肌电信号截取10s长度)，共200个样本数据，得出结果如图6，可以看出滤波后的肌电信号对于病变预正常信号有更好的区分度，说明本方法能更好分辨两个不同状态的肌电信号。

Claims (6)

1.基于分数阶傅立叶变换的自适应肌电信号检测处理方法，其特征在于,包括如下步骤：

1)利用信号采集装置，对肌电信号进行采集获取，并对采集得到的原始信号进行预处理，得预处理后的肌电噪声混合信号S_SMG+N(t)；

2)将步骤1中肌电噪声混合信号S_SMG+N(t)，预设定智能分段的初始参数；

3)对信号S_EMG+N(t)进行智能分段；对分段后的小段信号S_i的时频信息F_i进行处理并线性拟合以求得该小段信号最优分数阶傅立叶变换的旋转角度α_iopt和最优分数阶傅立叶变换的阶数p_iopt；

4)对截取后的剩余信号再进行截取且令i＝i+1，重复步骤3，直至得到整段S_EMG+N(t)肌电信号所有信号段的最优阶数p_iopt的集合p_opt(t)；对信号S_EMG+N(t)进行自适应p_opt(t)阶分数阶傅立叶变换和u域窗函数滤波；

5)对步骤4所得去除噪声后的u域上的肌电信号S_EMG(u)进行p_opt(t)阶分数阶傅立叶反变换，得到时域上的去除噪声的肌电信号S_EMG(t)；

6)对步骤5得到的去除噪声后的肌电信号S_EMG(t)提取特征值，判断肌肉所处的状态。

2.根据权利要求1所述的基于分数阶傅立叶变换的自适应肌电信号检测处理方法，其特征是，所述步骤1具体为：利用传感器和信号采集电路对肌电信号进行采集，并将采集到的数据进行放大、简单滤波和模数转化，得到肌电噪声混合信号；其中简单滤波为1.5倍增益中心频率为50Hz的巴特沃斯带阻滤波器，主要用以滤除50Hz的工频信号干扰。

3.根据权利要求1所述的基于分数阶傅立叶变换的自适应肌电信号检测处理方法，其特征是，所述步骤2具体为：设智能分段第i段分割时间长度为t_i，先预设t_i初值为总信号长度L的十分之一，即t_i＝L/10；若求得t_i存在小数，则进行四舍五入求整；按时间t_i分割出第i小段信号S_i(t)，简写为S_i；并进行p＝0阶傅立叶变换，得到信号的时频关系F_i。

4.根据权利要求1所述的基于分数阶傅立叶变换的自适应肌电信号检测处理方法，其特征是，所述步骤3具体为：定义瞬时信号绝对值|S_i|大于该小段信号绝对值的最大值|S_i|_max三分之一的肌电信号点称为该段信号的有效点，记为S_ief；定义有效点数为该小段信号S_i的有效点S_ief个数，记为N_i；

求S_i信号长度L_i与有效信号点数N_i，若N_i大于或等于整段信号长度L_i的2/3，或L_i小于等于100，则将该小段信号S_i的时频关系F_i进行线性拟合，得到拟合斜率k_i，取α_i＝tan^-1(-k_i)，得到分数阶傅立叶变换阶数p_i＝2α_i/π；反之，若该小段信号长度t_i为偶数则信号长度再截半，即t_i＝t_i/2，若t_i现为奇数，则信号长度加1再截半,即t_i＝(t_i+1)/2；将截取后的信号记为S_i重新进行本步骤，直至求出S_i的最优分数阶傅立叶变换的旋转角度α_iopt和最优分数阶傅立叶变换的阶数p_iopt。

5.根据权利要求1所述的基于分数阶傅立叶变换的自适应肌电信号检测处理方法，其特征是，所述步骤4具体为：设经过步骤3智能分段后将信号S_EMG+N(t)一共分为n段小段信号，分别为S_i，i＝1,2,3,...,n；

定义肌电信号最优阶数集合p_opt(t)为所有p_iopt按时间分布的分段函数，具体数学表示为p_opt(t_i)＝p_iopt；

对S_EMG+N(t)做p_opt(t)阶傅立叶变换，对每段信号S_i在u域上加窗函数f_i(u)，使每小段有效信号S_i能在u域的窗函数内，在u域上得到经过滤波后的肌电信号S_EMG(u)。

6.根据权利要求1所述的基于分数阶傅立叶变换的自适应肌电信号检测处理方法，其特征是，所述步骤6具体为：求去除噪声后的肌电信号的中值频率f_mid与信号的标准差σ，得出不同肌电信号的动作差异，以识别肌肉所处状态；其具体数值公式分别为：

$f_{mid}=\frac{1}{2}{\∫}_0^{+\mathrm{\∞}}PSD\left(f\right)df---\left(4\right)$

$\mathrm{\σ}=\frac{1}{L_S-1}\underset{t=1}{\overset{L_s}{\Σ}}{S}_{EMG}^2\left(t\right)---\left(5\right)$

其中，PSD(f)为肌电信号的功率谱密度函数，Ls为滤波后肌电信号的长度，S_EMG(t)为滤波后肌电信号的瞬时值。