赞
踩
EGG信号具有微弱、低幅值、低频、随杋性的特点,很容易被噪声干扰,而噪声可能来自生物体内,如呼吸、肌肉颤抖,也可能因为接触不良而引起体外干扰。是ECG信号主要的三种噪声为工频干扰、肌电干扰和基线漂移3,也是在滤波过程中急需被抑制去除的噪声干扰。
小波变换(Wavelet Transform, WT)可以进行时频变换,是对信号进行时域以及频域分析的最为理想工具。本文对含噪心电信号采用基于小波变换的去噪处理方法,分为以下3个步骤:
由于噪声和信号混杂在一起,首先选择一个小波基函数,由于噪声和信号混杂在一起,所以要用小波变换对含噪心电信号进行某尺度分解得到各尺度上的小波系数。
心电信号经过小波变换尺度分解后,幅值比较大的小波系数就是有用的信号,幅值比较小的小波系数就是噪声,根据心电信号和夹杂噪声的频率分布,对各尺度上的小波系数进行阈值处理,把小于阈值的小波系数置零或用阈值函数处理。
分别处理完小波尺度分解后的低频系数和高频系数,再重构信号。
4尺度小波分解所得各尺度系数示意图如下,9尺度小波分解可以类比之:
小波系数处理的阈值函数有硬阈值和软阈值之分。
其中w为原始小波系数,W为处理后的小波系数,λ为给定的阈值,N为信号长度。λ的计算公式为
λ
=
m
e
d
i
a
n
∣
w
∣
∗
2
∗
l
n
N
0.6745
\lambda=\frac{median|w|*\sqrt{2*lnN}}{0.6745}
λ=0.6745median∣w∣∗2∗lnN
median为中位数
由上文分析可知,软阈值去噪法的小波系数在连续性上优于硬阈值法,故本文采取了软阈值法结合小波基对信号进行仿真实验。
关于小波变换的介绍可以参考这篇文章:https://blog.csdn.net/tbl1234567/article/details/52662985
去噪理论大致如此,下面来实战编写代码。以编号为100的心电数据记录为例,上一篇文章已经通过wfdb将心电数据读取到了record变量中,record为一个大小为(650000,1)的numpy数组。要对其进行小波分解,首先要通过np.flatten()将其转换为shape=(650000,)的数组。
record = wfdb.rdrecord('ecg_data/' + 100, channel_names=['MLII'])
data = record.p_signal.flatten()
Python中使用pywt包进行信号的时频分析,其中包括傅里叶变换、小波变换等工具。使用pywt.wavedec()将时域信号进行9尺度变换到频域,小波基选用db5,返回值即为各尺度系数。
# 用db5作为小波基,对心电数据进行9尺度小波变换
coeffs = pywt.wavedec(data=data, wavelet='db5', level=9)
cA9, cD9, cD8, cD7, cD6, cD5, cD4, cD3, cD2, cD1 = coeffs
由于MIT-BIH心电信号的采样频率为360Hz,根据奈奎斯特采样定理,原始心电信号的最大频率在180Hz以下,因此信号分解的D1层最大频率为180Hz。在对原始信号进行分解之后,我们可以知道,1-2层的细节分量的能量与原始信号的高频干扰保持一致。表明1-2层是高频噪声(工频干扰以及肌电噪声)集中的主要地方。因此我们需要滤除D1层和D2层的细节分量,通过将其置0来达到去除的目的。然后将信号分解得到的3~9层小波系数通过软阈值公式对信号的阈值进行处理。pywt.threshold()函数提供了阈值滤波功能,并且默认是mode='soft’的软阈值滤波。
threshold = (np.median(np.abs(cD1)) / 0.6745) * (np.sqrt(2 * np.log(len(cD1))))
# 将高频信号cD1、cD2置零
cD1.fill(0)
cD2.fill(0)
# 将其他中低频信号按软阈值公式滤波
for i in range(1, len(coeffs) - 2):
coeffs[i] = pywt.threshold(coeffs[i], threshold)
最后对小波系数进行反变换,获得去噪后的信号。
rdata = pywt.waverec(coeffs=coeffs, wavelet='db5')
对预处理前后的100号信号截取前1500个信号点进行对比如下,可以看出降噪效果还是不错的:
plt.figure(figsize=(20, 4))
plt.plot(data)
plt.show()
plt.plot(rdata)
plt.show()
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。