CN101548885A - 去除电生理信号中工频干扰信号的方法 - Google Patents

Abstract

一种去除电生理信号中工频干扰信号的方法，步骤包括：1)对于采集到的电生理信号，把工频干扰信号看作正弦信号，重构该正弦信号，该正弦信号即为类似于工频干扰的正弦信号；2)从采集到的电生理信号中减去重构的正弦信号，即消除干扰，得到纯净的电生理信号。本方法用简单的减法应用就消除了干扰，而且效果很好。

Description

去除电生理信号中工频干扰信号的方法

技术领域：

本发明涉及去除工频干扰信号的方法，尤其是电生理信号中工频干扰信号的方法，本发明方法属于自适应信号处理范畴的降噪算法，是将锁定放大器软件化后的结果。

背景技术：

针对电生理信号中的工频干扰，传统的方法是使用50Hz的陷波器。然而，由于电网、用电设备的影响，因此传统的陷波器难以达到令人满意的结果。在IEEE Trans这些领先的杂志上，每年都有针对去除工频干扰信号的文章发表，但是，这些文章大多从随机信号处理的角度出发，虽然处理效果不错，但是往往存在算法复杂度高，难以理解等缺点。因此，在对电生理信号进行预处理时，大家都迫切地需要一种易于实现，易于理解，同时效果又相当不错的去除工频干扰的算法。

例如，在心电图检测中，现有技术中，在硬件上，尽管通过引进更好的仪表放大器以及右腿驱动电路，共模干扰能够被极大削弱，然而，硬件改良并非在所有场合下都适用。在软件上，传统的做法是采用50Hz陷波器。然而，工频干扰的频率相位并非一成不变，因此单纯地设计50Hz附近一个很窄的频段作为陷波器效果不佳。但是，武断地增大陷波器带宽的做法则会引起两个问题：一是50Hz附近的有用频率也会被滤除，造成信号失真，二是会在信号的时域图上产生振铃现象。

发明内容

为了解决现有技术中存在的上述问题，本发明提供了一种新的去除电生理信号中工频干扰信号的方法，具体技术方案如下：

一种去除电生理信号中工频干扰信号的方法，步骤包括：

1)对于采集到的电生理信号，把工频干扰信号看作正弦信号，重构该正弦信号(该正弦信号即为类似于工频干扰的正弦信号)；

2)从采集到的电生理信号中减去重构的正弦信号，即消除干扰，得到纯净的心电图信号。

所述步骤1)中，探测正弦信号的频率、相位和幅度参量，重构正弦信号，步骤是，把工频干扰信号看作有用信号，而将电生理信号看作干扰信号，通过抑制电生理信号，精确地检测工频干扰的幅度、频率、相位，然后重构干扰信号，具体来说，步骤包括：

1.1)设定一个初始参考信号，幅度、频率和相位都已知；

1.2)然后把含有干扰的电生理输入信号Input X(t)和初始参考信号Reference作比较，不断观察比较得到的Delta Frequency，Delta Phase和两个信号幅度的乘积；

1.3)如果Delta Frequency或Delta Phase大于相应的阀值，就根据DeltaFrequency和Delta Phase不断修正参考信号；所述Delta Frequency是两个信号频率的差值，Delta Phase是两个信号相位的差值；

1.4)重复步骤1.1～1.3，输出的Delta Frequency和Delta Phase越来越小，最终小于阈值，此时，参考信号的频率、相位非常逼近干扰信号的频率相位；

1.5)应用步骤1.4修正最终得到的参考信号的Frequency、Phase作为干扰信号的频率和相位，参考信号和输入信号的幅度乘积，作为干扰信号的幅度，重构仿真得到干扰信号。

所述步骤2)中，用输入信号Input X(t)减去重构仿真的干扰信号，去除输入信号中的干扰信号，得到纯净的电生理信号。

所述步骤1)是在锁定放大器中完成的，所述锁定放大器包括两个乘法器和两个低通滤波器LPF；每个乘法器对应设在一个低通滤波器LPF的前端，所述心电图输入信号Input X(t)分别传入一个乘法器；所述初始参考信号Reference有两组，分别是一对正交的信号Vref和Vref90°；所述Vref和Vref90°分别传入一个乘法器中；测量两个低通滤波器LPF输出信号，运算后即得干扰信号的幅度、相位和频率。

本发明还将所述锁定放大器的硬件电路软件化，使本发明更具有灵活性和稳定性。

所述Delta Frequency的阈值y＝0.005Hz，Delta Phase的阈值是0.005 Radian。

本发明把工频干扰(就是市电50Hz干扰)看成正弦波(实际情况是类正弦波而不是绝对的正弦波)，然后通过一定的方法，将这个正弦波重构出来，最后再将这个重构的正弦波减去，达到了消除干扰的目的。本发明用简单的减法应用就消除了干扰。而在信号处理领域中，几乎所有的去除干扰的行为，就可以概括为卷积滤波，业内人士也对这种方法非常熟悉，甚至已经产生了审美疲劳或是思维惯性。但是，本发明与传统不同，就是简单的相减，而且效果很好。

本发明在探测频率、相位、幅度时，在宏观上构成了一个反馈环，通过多次反馈，逐渐逼近了那三个量的精确值。由于是基于反馈原理，因此最终的精度非常的高，可以达到小数点后6位。这种利用反馈的思想，将探测精度提高了不止3个数量级。

将锁定放大器的硬件电路软件化，由于锁定放大器的强烈的抑制噪声的本领，本算法工作时，不需要避开QRS波群。

附图说明

图1是本方法原理示意图

图2是锁定放大器原理框图

图3是测试一中原始心电信号

图4是测试一中含有工频干扰的心电信号

图5是测试一中还原的信号(分段长度为128个点)

图6是测试一中干扰信号(分段长度为128个点)

图7是测试二中实际采集的ECG信号

图8是测试二中还原的信号。

具体实施方式

一种去除电生理信号中工频干扰信号的方法，步骤包括：

1)对于采集到的电生理信号，把工频干扰信号看作正弦信号，重构该正弦信号(该正弦信号即为类似于工频干扰的正弦信号)；

2)从采集到的电生理信号中减去重构的正弦信号，即消除干扰，得到纯净的心电图信号。

所述步骤1)中，探测正弦信号的频率、相位和幅度参量，重构正弦信号，步骤是，把工频干扰信号看作有用信号，而将电生理信号看作干扰信号，通过抑制电生理信号，精确地检测工频干扰的幅度、频率、相位，然后重构干扰信号，具体来说，步骤包括：

1.1)设定一个初始参考信号，幅度、频率和相位都已知；

1.2)然后把含有干扰的电生理输入信号Input X(t)和初始参考信号Reference作比较，不断观察比较得到的Delta Frequency，Delta Phase和两个信号幅度的乘积；

1.3)如果Delta Frequency或Delta Phase大于相应的阀值，就根据DeltaFrequency和Delta Phase不断修正参考信号；所述Delta Frequency是两个信号频率的差值，Delta Phase是两个信号相位的差值；

1.4)重复步骤1.1～1.3，输出的Delta Frequency和Delta Phase越来越小，最终小于阈值，此时，参考信号的频率、相位非常逼近干扰信号的频率相位；

1.5)应用步骤1.4修正最终得到的参考信号的Frequency、Phase作为干扰信号的频率和相位，参考信号和输入信号的幅度乘积，作为干扰信号的幅度，重构仿真得到干扰信号。

所述步骤2)中，用输入信号Input X(t)减去重构仿真的干扰信号，去除输入信号中的干扰信号，得到纯净的电生理信号。

所述步骤1)是在锁定放大器中完成的，所述锁定放大器包括两个乘法器和两个低通滤波器LPF；每个乘法器对应设在一个低通滤波器LPF的前端，所述心电图输入信号Input X(t)分别传入一个乘法器；所述初始参考信号Reference有两组，分别是一对正交的信号Vref和Vref90°；所述Vref和Vref90°分别传入一个乘法器中；测量两个低通滤波器LPF输出信号，运算后即得干扰信号的幅度、相位和频率。

本发明还将所述锁定放大器的硬件电路软件化，使本发明更具有灵活性和稳定性。

所述Delta Frequency的阈值y＝0.005Hz，Delta Phase的阈值是0.005 Radian。

下面以本方法在去除心电图中工频干扰信号中的应用，结合附图对本发明作进一步说明。

参考图1，输入信号包括两部分：ECG信号和工频干扰。此时ECG被看作干扰信号，而工频干扰则被看成有用的信号。假设存在一个装置，他能够抑制ECG信号，比较干扰信号和参考信号Vref，从而得到输入信号中干扰信号和参考信号二者的频率差，相位差和幅度的乘积。

具体来说，首先设定一个初始参考信号(幅度，频率和相位都已知)，然后把含有干扰的输入信号Input X(t)和初始参考信号Reference交给这所述装置，然后不断观察输出的Delta Frequency，Delta Phase和两个信号幅度的乘积。如果Delta Frequency或Delta Phase过大(意味着干扰信号和参考信号reference的差别很大)，就根据Delta Frequency和Delta Phase不断修正参考信号，然后再把修正过的参考信号reference和Input送给所述装置，重复以上的过程。这样经过若干次反复以后，输出的Delta Frequency和Delta Phase将越来越小，最终小于阈值(本发明中设为阈值Delta Frequency＝0.005Hz，阈值Delta Phase＝0.005Radian)，而此时就可以认为参考信号的频率、相位非常逼近干扰信号的频率相位。而又可以通过输出的两个信号的幅度乘积，得到干扰信号的幅度。最后，可以使用已被修正好的参考信号的Frequency，Phase和干扰信号的幅度这三个参数，重构一个仿真的干扰信号，然后只需要用输入信号减去仿真的干扰信号，就可以去除输入信号中的干扰信号。

所述装置选用锁定放大器，如图2，Vref和Vref90°是一对正交的信号，Vref就是所述初始参考信号，LPF是低通滤波器。

参考图2，对本算法作进一步说明。

先设X(t)为单频信号(即输入的信号是单纯的工频干扰信号)，并且可以写成

Vref和Vref90°为参考信号，所有初始参数均由人为设定，设初始参考信号的形式是，

那么，经过乘法器以后，A和B处的信号为：

式中

和

为高频信号，在经过低通滤波器后就会被滤除，所以，C和D处的信号分别为：

最后一步的运算是：

V1，V2可看成两条曲线，从这两条曲线中，可以得到Delta Frequency(频率差)和Delta Phase(相位差)的信息。通过测量V1，便可以从(1)式中得到输入信号X(t)的幅度Ei， $\mathrm{Ei}=\frac{2V1}{\mathrm{Er}}.$

对于2)式，在t比较小的时候，可以近似看成：

式中，DeltaP即为DeltaPhase，Delta F即为DeltaFrequency。

当t＝0时，即在V2曲线的初始点，V2＝Delta P等于干扰信号和参考信号的相位差。而 $\frac{d\left(v2\right)}{\mathrm{dt}}=2\mathrm{\π}\×\mathrm{DeltaF},$ 即V2曲线的斜率，就等于干扰信号和参考信号的频率差。

这样，通过测量V1和V2，便可以得到干扰信号的幅度、干扰信号和参考信号的相位差，频率差。

但是，实际情况中，输入信号将不是一个单频信号，而是ECG信号和干扰信号的混合体。为了将X(t)普适化，令X(t)为一个单频信号和多个频率信号的混杂体，这样更接近实际采集到的信号，设

其中

相当于噪声(在实际中相当于ECG波形)参考信号不变，设为：

那么，A和B处对应的信号就变为：

同样的道理，经过低通滤波器之后，斜体字部分的高频信号将被滤除。而对于频率为f_n1-f₂，f_n2-f₂，…，f_nx-f₂的信号来说，除非f_n1-f₂，f_n2-f₂，…，f_nx-f₂非常接近直流，否则仍然会被滤除，无法对后续产生干扰，对于心电信号来说，能量主要集中在30Hz以下，50Hz附近的能量非常小，很难对探测结果产生影响，因此可以忽略这种计算中的干扰。

在LPF的截止频率很低的情况下，

算法改进：

上面的部分已经将算法数学化，清晰地表明了理论推导部分和可能产生误差的部分，包括式(3)、(4)和(5)。以上的理论推导出的信号，在实际中会遭到很强的干扰，以致无法直接应用上述算法。因此，在实际算法中还需消除误差和抑制噪声，精确地检测到干扰信号。

本方法提出了acquisition(大致获得)和tracking(跟踪)的方案：

在Acquisition的过程中，精度是次要的，而获得Delta Frequency和Delta Phase的大致值则是其主要目的，因此算法需要鲁棒性(稳定性)。在Acquisition的过程中，为了得到Delta Phase的大致值，首先令t＝0，求出Delta Phase＝V2。为了求Delta Frequency的大致值，只分析V2序列的前20个点(t比较小，保证式7)中约等号不致引起过大误差，保证算法的稳定性)，LPF的截止频率也被设定为5Hz，这意味至少在理论上，算法可以追踪从45Hz-55Hz的干扰信号。然后，立即使用Acquisition过程中得到的Delta Frequency和Delta Phase修正参考信号。

由于在Acquisition的阶段中，Delta Frequency和Delta Phase的大致值已经确定，因此，经过修正以后的参考信号与要追踪的干扰信号大大接近了，即此时的Delta Frequency和Delta Phase已经比Acquisition以前大大缩小了。故可以进入Tracking过程。

在Tracking的过程中，精度非常重要。随着Delta Frequency和Delta Phase的缩小，在中(Delta P即为Delta Phase)，即使t取得比较大，约等号的精度仍然可以保证，这就允许分析V2中更多的数据。并且，随着Tracking过程的不断迭代，DeltaFrequency会越来越小，因此我们就可以分析V2中越来越多的数据。Acquisition中的仅取V2中前20个数据是为了保证算法的稳定性，而Tracking中随着DeltaFrequency的降低，分析V2中越来越多的数据则可以平滑信号中的噪声，保证了算法的精度。

另外，由于在Tracking过程中，Delta Frequency不断缩小，因此低通滤波器的截止频率就可以随着Delta Frequency的缩小而缩小，以便取得更强的抑制干扰的能力。

在算法中，Acquisition和Tracking过程相当于一个反馈环：Acquisition首先大致定位，然后不断重复Tracking的过程。在Tracking中，随着Delta Frequency的减小，V2中越来越多的点纳入计算之中，LPF的截止频率也随Delta Frequency的减小而降低。这两样措施都会进一步提高精度，抑制误差，反过来又减小了Delta Frequency的值。这个反馈和模拟电路中的反馈机制有类似的作用，可以不断消除误差，使参考信号越来越逼近干扰信号，在模拟电路中反馈环付出的代价是使电路放大倍数减小，而Tracking过程付出的代价则是更多的计算量。

测试：

测试一(参考图3～6)，取MIT-BIH数据库的一段信号(采样频率256Hz，长度为6s)，这段信号受到了人为添加的幅度、频率、相位不同的噪声的污染。

表1：

表2：

表3：

从图4中，还原信号的第1025点到1152点的噪声比其他段大，也比较明显。这是因为，对原始的第1025点到1152点添加的噪声幅度达到5V(所有噪声中最大)，因此与之对应的还原信号中，噪声自然也相应较大。不过，从表格中可以看出，算法的效果(SNR改进比)，在第5段没有明显下降，约为35dB。

测试二(参考图7、8)：

再取一段实际采样得到的心电信号.采样频率为1000Hz，信号持续时间为1560ms。应用本方法后，SNR(信噪比)改善程度：30.26dB

总结：

测试证明，本算法继承了锁定放大器出色的抑制噪声的能力。即使在工频干扰幅度很大的情况下，(如测试一中，信号的最大值仅为1V而噪声的幅度达1.5V甚至5V，测试二中ECG信号几乎全部淹没在工频干扰的噪声中)，算法仍能恢复出ECG信号。由于采用简单抵消的方法，而非卷积，相比传统的陷波器，我们的算不存在相移和振铃现象。本方法反应速度更快，几乎不需要学习的时间。很多去除50Hz工频干扰的算法都需要做QRS检测，并要小心翼翼地避开QRS波以防止干扰，然而，由于LIA(锁定放大器)出色的抗噪声能力，本算法并不需要做QRS检测，大大地简化了算法流程。

由于参考信号的频率可以任意调节，因此从理论上来说，本算法不但可以消除50Hz的工频干扰，亦可以消除60Hz的干扰和100Hz的倍频干扰。这种算法的效果相当于一个效果明显的陷波器，却没有滤波器的相移作用，也不会像陷波器那样对50Hz附近频率的信号也有影响，而是仅仅对工频干扰有去除作用。由于仅仅涉及了初等数学运算，这种算法无疑要比很多复杂的算法容易理解，容易实现得多。

Claims (9)

1、一种去除电生理信号中工频干扰信号的方法，其特征是步骤包括：

1)对于采集到的电生理信号，把工频干扰信号看作正弦信号，重构该正弦信号，该正弦信号即为类似于工频干扰的正弦信号；

2)从采集到的电生理信号中减去重构的正弦信号，即消除干扰，得到纯净的电生理信号。

2、根据权利要求1所述的去除电生理信号中工频干扰信号的方法，其特征是所述步骤1)中，探测正弦信号的频率、相位和幅度参量，重构正弦信号，步骤是，把工频干扰信号看作有用信号，而将电生理信号看作干扰信号，通过抑制电生理信号，精确地检测工频干扰的幅度、频率、相位，然后重构干扰信号。

3、根据权利要求2所述的去除电生理信号中工频干扰信号的方法，其特征是所述步骤1)中，检测工频干扰的幅度、频率、相位的的步骤包括：

1.1)设定一个初始参考信号，幅度、频率和相位都已知；

1.2)然后把含有干扰的电生理输入信号Input X(t)和初始参考信号Reference作比较，不断观察比较得到的Delta Frequency，Delta Phase和两个信号幅度的乘积；

1.3)如果Delta Frequency或Delta Phase大于相应的阀值，就根据DeltaFrequency和Delta Phase不断修正参考信号；所述Delta Frequency是两个信号频率的差值，Delta Phase是两个信号相位的差值；

1.4)重复步骤1.1～1.3，输出的Delta Frequency和Delta Phase越来越小，最终小于阈值，此时，参考信号的频率、相位非常逼近干扰信号的频率相位；

1.5)应用步骤1.4修正最终得到的参考信号的Frequency、Phase作为干扰信号的频率和相位，参考信号和输入信号的幅度乘积，作为干扰信号的幅度，重构仿真得到干扰信号。

4、根据权利要求1或2或3所述的去除电生理信号中工频干扰信号的方法，其特征是所述步骤2)中，用输入信号Input X(t)减去重构仿真的干扰信号，去除输入信号中的干扰信号，得到纯净的电生理信号。

5、根据权利要求3所述的去除电生理信号中工频干扰信号的方法，其特征是所述Delta Frequency的阈值y＝0.005Hz，Delta Phase的阈值是0.005Radian。

6、根据权利要求2或3所述的去除电生理信号中工频干扰信号的方法，其特征是所述步骤1)是在锁定放大器中完成的，所述锁定放大器包括两个乘法器和两个低通滤波器LPF；每个乘法器对应设在一个低通滤波器LPF的前端，所述心电图输入信号Input X(t)分别传入一个乘法器；所述初始参考信号Reference有两组，分别是一对正交的信号Vref和Vref90°；所述Vref和Vref90°分别传入一个乘法器中；测量两个低通滤波器LPF输出信号，运算后即得干扰信号的幅度、相位和频率。

7、根据权利要求6所述的去除电生理信号中工频干扰信号的方法，其特征是将所述锁定放大器的硬件电路软件化，算法的步骤包括：

X(t)为电生理信号和工频信号的混合体，即单频信号和多个频率信号的混杂体，设

其中

相当于噪声，即为电生理信号波形；

设参考信号为：

经过乘法器以后，Vref和Vref90°分别对应的信号为：

和

；

经低通滤波后，高频信号被滤除，得到Vref和Vref90°分别对应的信号为：

和

；

对于频率为f_n1-f₂，f_n2-f₂，…，f_nx-f₂的信号来说，除非该信号非常接近直流，否则仍然会被滤除，无法对后续产生干扰，对于电生理信号来说，能量主要集中在30Hz以下，50Hz附近的能量非常小，很难对探测结果产生影响，因此可以忽略这种计算中的干扰；

                                $V1=\sqrt{C^2+D^2}\≈\frac{\mathrm{EiEr}}{2}$

在LPF的截止频率很低的情况下，       ；

                               

V1，V2可看成两条曲线，通过这两条曲线中，可以得到信号X(t)和参考信号的频率差Delta Frequency和相位差Delta Phase的信息；

对于V1，信号X(t)的幅度 $\mathrm{Ei}=\frac{2V1}{\mathrm{Er}};$

对于V2，在t比较小的时候，可以近似看成：

式中Delta F即为Delta Frequency，Delta P即为Delta Phase；

当t＝0时，即在V2曲线的初始点，V2＝Delta Phase等于干扰信号和参考信号的相位差。而 $\frac{d\left(v2\right)}{\mathrm{dt}}=2\mathrm{\π}\×\mathrm{DeltaF}$ 即V2曲线的斜率，即为X(t)信号和参考信号的频率差；

通过测量V1和V2，即得工频干扰信号的幅度、相位和频率。

8、根据权利要求7所述的去除电生理信号中工频干扰信号的方法，其特征是所述算法中还包括acquisition大致获得和tracking跟踪过程；Acquisition和Tracking过程相当于一个反馈环：Acquisition首先大致定位，然后不断重复Tracking的过程，步骤包括：

a)先进行Acquisition的过程：

求Delta Phase的大致值：

首先，令t＝0，V2＝Delta Phase，只分析V2序列的前20个点，LPF的截止频率也被设定为5Hz；然后，立即使用Acquisition过程中得到的Delta Frequency和Delta Phase修正参考信号；

b)再进入Tracking的过程，并不断重复Tracking的过程：

随着Delta Frequency和Delta Phase的缩小，在

中，即使t取得比较大，精度仍然可以保证，则分析V2中更多的数据；并且，随着Tracking过程的不断迭代，Delta Frequency会越来越小，因此可分析V2中越来越多的数据。

9、根据权利要求1或2所述的去除电生理信号中工频干扰信号的方法，其特征是所述算法在不避开QRS波群和情况下进行。

Images