基于分数阶傅立叶变换的自适应肌电信号检测处理方法
【专利摘要】本发明公开了一种基于分数阶傅立叶变换的自适应肌电信号检测处理方法。该方法首先对信号进行放大和滤去工频干扰等处理，利用智能分段，不断地将信号分段，同时进行自适应分数阶傅立叶变换滤波，在分数阶傅立叶域上将噪声信号滤除，再通过分数阶傅立叶反变换将信号还原至时域，最后求取肌电信号的特征值判断肌电信号所表示的肌肉状态。本发明采用分数阶傅立叶变换滤波算法分析肌电信号，更好地利用了肌电信号的时频信息，更好地分辨出了两种不同状态下的肌电信号的区别。
【专利说明】
基于分数阶傅立叶变换的自适应肌电信号检测处理方法
技术领域
[0001] 本发明属于数字信号处理和医学领域，具体涉及一种基于分数阶傅立叶变换的自 适应肌电信号检测处理方法。
【背景技术】
[0002] 表面肌电信号是生物神经肌肉系统活动时所产生的生物电信号，它同生物神经肌 肉的活动密切相关，其中包含着许多的肌体运动的信息。对表面肌电信号的监测和分析技 术已经广泛地应用于临床医学、生物医学、运动医学等多个领域。由于肌电信号是肌肉神经 的生理活动产生的，它具有典型非线性、低信噪比、非平稳的特征。
[0003]自1807年法国数学家傅立叶率先提出傅里叶分析的数学方法来分析信号频域信 息以来，傅立叶变换已经得到了迅速而广泛的推广与应用。但是，由于傅立叶变换是一种全 局性的信号的变换，所以难以去表述信号的局部信息，特别是信号的局部时频信息。而在一 些非平稳信号的分析与应用中，此种信号的局部时频信息恰恰是其中至关重要、不容忽视 的部分。
[0004] 传统的傅立叶变换是一种线性算子，可看成是将信号从时间轴逆时针旋转V2到 频率轴。而分数阶傅立叶则是可旋转任意角度α的算子，从而可以得到信号在任意角度下的 分析，而不是仅仅局限于频率轴上的信号分析方法。在时域上，依据Wigner-Ville分布对于 分数阶傅立叶变换的旋转不变性，分数阶傅立叶可以看作是信号在时域平面绕原点旋转任 意角度α后所构成的分数阶傅立叶域u上的表示，其与Wigner-Vi lie变换、短时傅立叶变换 之间存在的只是坐标旋转变化关系，并不影响信号的时域分布特性。这个性质就可以在非 平稳过程的情况下，利用分数阶傅立叶域来实现时频滤波，从而在减少计算开支的前提下 也可以得到更好的滤波效果。
【发明内容】
[0005] 本发明的目的在于针对现有技术的不足，提供一种基于分数阶傅立叶变换的自适 应肌电信号检测处理方法。
[0006] 本方法所采用的技术方案，包括以下步骤：
[0007] 1)利用信号采集装置，对肌电信号进行采集获取，并对采集得到的原始信号进行 预处理，得预处理后的肌电噪声混合信号s SMG+N(t)。
[0008] 2)将步骤1中采集得到的预处理后的肌电噪声混合信号SSMG+N(t)，预设定智能分段 的初始参数。
[0009] 3)对信号SEMG+N(t)进行智能分段。对分段后的小段信号31的时频信息？,进行处理 并线性拟合以求得该小段信号最优分数阶傅立叶变换的旋转角度<^_和最优分数阶傅立 叶变换的阶数Pi 〇Pt。
[0010] 4)对截取后的剩余信号再进行截取且令i = i + l，重复步骤3,直至得到整段SEMG+N ⑴肌电信号所有信号段的最优阶数p_t的集合Pcipt(t)。对信号SEMC+N(t)进行自适应 Pcipt(t) 阶分数阶傅立叶变换和U域窗函数滤波。
[0011] 5)对步骤4所得去除噪声后的u域上的肌电信号SEMG(u)进行Pcipt(t)阶分数阶傅立 叶反变换，得到时域上的去除噪声的肌电信号S EMC(t)。
[0012] 6)对步骤5得到的去除噪声后的肌电信号SEMC(t)提取特征值，判断肌肉所处的状 ??τ 〇
[0013] 本发明的有益效果在于:选用了分数阶傅立叶变换，相较于传统傅立叶变换，充分 利用了输入肌电信号的局部时频信息。并且针对不同的输入信号，利用了拟合算法自适应 地调节分数阶傅立叶的u域旋转角度α，从而更好得适应复杂的肌电信号输入。
【附图说明】
[0014] 图1本发明整体实现过程框图；
[0015] 图2预处理后两种不同状态下的肌电信号图；
[0016] 图3原始肌电信号时频图；
[0017] 图4分数阶傅立叶变换后肌电信号时频图；
[0018] 图5分数阶傅立叶变换滤波算法后两种不同状态肌电信号图；
[0019] 图6多个数据组滤波前后肌电信号特征值分布图。
【具体实施方式】
[0020] 传统的肌电信号模型，可分成两部分，一部分是有用信号ρ(Λ)另外一部分是噪声 信号e，即其数学模型如下：
[0021 ] EMG -(p(j(；) + e (1)
[0022] 对于该肌电信号在傅立叶域上的滤波，
[0023] x = f(EMG) (2)
[0024] 由于肌电信号为非平稳过程，其滤波算子f 一般不为时不变，也无法表示为时域卷 积或传统傅立叶域上的乘法滤波器，利用时变算子f(t)对原始肌电信号进行滤波。为了得 到时变算子f(t)，本实施例将经过放大处理后的待滤波信号进行短时拆分，并对其进行p阶 分数阶傅立叶变换(FRFT)。分数阶傅立叶变换的定义式可以写成如下形式：
[0025]
[0026]
[0027] 在短时间内，将待滤波的肌电信号看做为平稳过程，其滤波算子f(t)采用时不变 算子f。从而将时间不断递推，最终得到时变的滤波算子f(t)。本实施例是对肌电信号做p阶 分数阶傅立叶变换，在分数阶傅立叶域(u域)对信号进行滤波，所以滤波算子为f(u)。
[0028] 下面，将结合附图1对该处理肌电信号的方法进行解释。
[0029] 1)利用信号采集装置，对肌电信号进行采集获取，并对采集得到的原始信号进行 预处理，得预处理后的肌电噪声混合信号SSMG+N(t);
[0030]具体为:利用传感器和信号采集电路对肌电信号进行采集并将采集到的数据进行 放大、简单滤波和模数转化等预处理，得到预处理后的肌电噪声混合信号。其中简单滤波为 1.5倍增益中心频率为50Hz的巴特沃斯带阻滤波器，主要用以滤除50Hz的工频信号干扰。
[0031] 2)将步骤1中采集得到的预处理后的肌电噪声混合信号SSMC+N(t)，预设定智能分段 的初始参数；
[0032] 具体为:设智能分段第i段分割时间长度为tl，先预设^初值为总信号长度L的十分 之一，gPtl = L/10。若求得。存在小数，则进行四舍五入求整。按时间^分割出第i小段信号 SKt)，简写为&。并进行p = 0阶傅立叶变换，得到信号的时频关系K。
[0033] 3)对信号SEMG+N(t)进行智能分段。对分段后的小段信号51的时频信息？1进行处理 并线性拟合以求得该小段信号最优分数阶傅立叶变换的旋转角度<^_和最优分数阶傅立 叶变换的阶数Pi〇Pt;
[0034] 具体为:定义瞬时信号绝对值| Si |大于该小段信号绝对值的最大值| Si | max三分之 一的肌电信号点称为该段信号的有效点，记为Sirf。
[0035] 定义有效点数为该小段信号Si的有效点Sief个数，记为Ni。
[0036] 求51信号长度U与有效信号点数化，若K大于或等于整段信号长度U的2/3,或LVJ、 于等于100,则将该小段信号&的时频关系Fi进行线性拟合，得到拟合斜率匕，取 k〇,则根据公式3,分数阶傅立叶变换阶数?1 = 2〇1/31。反之，若该小段信号长度^为偶数则 信号长度再截半，即^ = "/2,若"现为奇数，则信号长度加1再截半，8卩^ = (^+1)/2。将截 取后的信号记为31重新进行本步骤，直至求出最优分数阶傅立叶变换的旋转角度a lcipt 和最优分数阶傅立叶变换的阶数Piopt。
[0037] 4)对截取后的剩余信号再进行截取且令i = i + l，重复步骤3,直至得到整段SEMG+N ⑴肌电信号所有信号段的最优阶数p_t的集合Pcipt(t)。对信号SEMC+N(t)进行自适应 Pcipt(t) 阶分数阶傅立叶变换和u域窗函数滤波；
[0038] 具体为:假设经过步骤3智能分段后将信号SEMG+N(t) -共分为η段小段信号，分别为 Si，i -1，2,3，···，n〇
[0039] 定义肌电信号最优阶数集合p〇pt(t)为所有?1_按时间分布的分段函数，具体数学 ^小为P〇pt(ti) - pi〇pt，i -1，2,3，· · ·，n〇
[0040] 对SEMG+N(t)做Pcipt(t)阶傅立叶变换，对每段信号Si在U域上加窗函数fi(U)，使每小 段有效信号Si能在u域的窗函数内，在u域上得到经过滤波后的肌电信号S EMC(u)。
[0041] 5)对步骤4所得去除噪声后的u域上的肌电信号SEMG(u)进行Pcipt(t)阶分数阶傅立 叶反变换，得到时域上的去除噪声的肌电信号SEMC(t);
[0042] 6)对步骤5得到的去除噪声后的肌电信号SEMG(t)提取特征值，判断肌肉所处的状 ??τ 〇
[0043] 具体为:求去除噪声后的肌电信号的中值频率fmid、与信号的标准差σ等特征值，得 出不同肌电信号的动作差异，以识别肌肉所处状态。其具体数值公式分别为：
[0046] 兵甲，PSLK f) 73肌电1目亏的功準堵想、皮凼数，Ls73滤波后肌电信号的长度，Semg (t) 为滤波后肌电信号的瞬时值，由于一般肌电信号的均值接近0,为了便于计算，这里将肌电
[0044]
[0045] 信号的均值设为0。
[0047]以下通过实例对本
【发明内容】
进行进一步解释，用该方法对实际肌电信号进行检测 处理。由于实际采集环境的工频噪声干扰等原因，淹没了肌电信号的特征成分，模糊了肌电 信号的特征值，如图2所示。对于两种不同状态(正常和病变)的肌肉状态也较难区分(信号1 长度为34.154s，信号2长度为38.769s，为了便于观察信号波形，这里都截取两信号第1秒时 间内的信号做图）。本方法对原始肌电信号进行采集放大等初步处理后，采用自适应分数阶 傅立叶变换滤波算法对初步处理后的肌电信号进行分析处理，这里以信号1为例，该段信号 长度L为34154ms，采样频率SP为20kHz。故预设分割时采用自适应拟合截取 算法得到的分数阶傅立叶域时频特性如图3所示。自适应拟合得到p = 0.42，做p = 0.42的分 数阶傅立叶变换后，信号时频特性如图4所示。最后，算法不断截取肌电信号得到?1^，最后 得到整段的时变的分数阶傅立叶阶数P〇 Pt(t)，滤波后在进行信号特征值分析，可以看出比 滤波前有较大的区分度，见图5,这里采取了 20个人体肱二头肌肌电信号样本(正常和病变 分别为10个，每个人10个段肌电信号数据，每个肌电信号截取l〇s长度），共200个样本数据， 得出结果如图6,可以看出滤波后的肌电信号对于病变预正常信号有更好的区分度，说明本 方法能更好分辨两个不同状态的肌电信号。
【主权项】
1. 基于分数阶傅立叶变换的自适应肌电信号检测处理方法，其特征在于，包括如下步 骤： 1) 利用信号采集装置，对肌电信号进行采集获取，并对采集得到的原始信号进行预处 理，得预处理后的肌电噪声混合信号SSMG+N(t); 2) 将步骤1中肌电噪声混合信号SSMC+N(t)，预设定智能分段的初始参数； 3) 对信号SEMC+N(t)进行智能分段;对分段后的小段信号S1的时频信息巧进行处理并线性 拟合以求得该小段信号最优分数阶傅立叶变换的旋转角度Ci lcipt和最优分数阶傅立叶变换 的阶数Piopt; 4) 对截取后的剩余信号再进行截取且令i = i + l，重复步骤3,直至得到整段SEMG+N(t)肌 电信号所有信号段的最优阶数的集合P〇Pt(t);对信号S_ N⑴进行自适应Pcipt(t)阶分 数阶傅立叶变换和u域窗函数滤波； 5) 对步骤4所得去除噪声后的u域上的肌电信号Semg(U)进行Pcipt(t)阶分数阶傅立叶反 变换，得到时域上的去除噪声的肌电信号S EMC(t); 6) 对步骤5得到的去除噪声后的肌电信号SEMC(t)提取特征值，判断肌肉所处的状态。2. 根据权利要求1所述的基于分数阶傅立叶变换的自适应肌电信号检测处理方法，其 特征是，所述步骤1具体为:利用传感器和信号采集电路对肌电信号进行采集，并将采集到 的数据进行放大、简单滤波和模数转化，得到肌电噪声混合信号;其中简单滤波为1.5倍增 益中心频率为50Hz的巴特沃斯带阻滤波器，主要用以滤除50Hz的工频信号干扰。3. 根据权利要求1所述的基于分数阶傅立叶变换的自适应肌电信号检测处理方法，其 特征是，所述步骤2具体为:设智能分段第i段分割时间长度为t，先预设^初值为总信号长 度L的十分之一，SPt 1 = IVlO;若求得。存在小数，则进行四舍五入求整;按时间。分割出第i 小段信号一乂^，简写为S1;并进行p = 0阶傅立叶变换，得到信号的时频关系K。4. 根据权利要求1所述的基于分数阶傅立叶变换的自适应肌电信号检测处理方法，其 特征是，所述步骤3具体为:定义瞬时信号绝对值I S11大于该小段信号绝对值的最大值I S1 Imax三分之一的肌电信号点称为该段信号的有效点，记为Slef;定义有效点数为该小段信号 Si的有效点Sief个数，记为Ni; 求31信号长度L1与有效信号点数化，若N1大于或等于整段信号长度L 1的2/3,或L1小于等 于100,则将该小段信号31的时频关系F1进行线性拟合，得到拟合斜率1^，取cmanl-lu)， 得到分数阶傅立叶变换阶数P 1 = Sa1A;反之，若该小段信号长度^为偶数则信号长度再截 半，即tFW%若0见为奇数，则信号长度加1再截半，即^=(^+1)/2;将截取后的信号记 为3 1重新进行本步骤，直至求出S1的最优分数阶傅立叶变换的旋转角度〇1￥和最优分数阶 傅立叶变换的阶数Pi 〇pt。5. 根据权利要求1所述的基于分数阶傅立叶变换的自适应肌电信号检测处理方法，其 特征是，所述步骤4具体为:设经过步骤3智能分段后将信号S EMG+N⑴一共分为η段小段信号， 分别为 Si，i = l，2,3，...，n; 定义肌电信号最优阶数集合P〇Pt(t)为所有p_t按时间分布的分段函数，具体数学表示 ^JPopt ( ti ) - Piopt ； 对SEMG+N(t)做Pcipt(t)阶傅立叶变换，对每段信号Si在U域上加窗函数fi(U)，使每小段有 效信号S1能在u域的窗函数内，在u域上得到经过滤波后的肌电信号Semc(U)。6.根据权利要求1所述的基于分数阶傅立叶变换的自适应肌电信号检测处理方法，其 特征是，所述步骤6具体为:求去除噪声后的肌电信号的中值频率fmid与信号的标准差得 出不同肌电信号的动作差异，以识别肌肉所处状态;其具体数值公式分别为：(4) (5) 其中，PSD(f)为肌电信号的功率谱密度函数，Ls为滤波后肌电信号的长度，SEMG(t)为滤 波后肌电信号的瞬时值。
【文档编号】A61B5/00GK105997067SQ201610453894
【公开日】2016年10月12日
【申请日】2016年6月21日
【发明人】郑永军, 黄强, 祝增献, 张聪, 李文军
【申请人】中国计量大学