超入门级-基于中值滤波处理ECG信号的基线漂移-Python-MIT-BIH数据集

中值滤波处理心电信号的基线漂移

距离上次发东西已经8个月,我已经本科毕业成为了一名研究生,但是我已经暂时弃硬从软,暂时开始做深度学习方向了,这篇文章就算一个我研究生学习的第一次笔记分享,我也会争取写的非常详细,ECG信号的处理是我第一个课题中的一部分,因为现实中的信号在采集的过程中都会存在噪声和干扰,所以如果想用ECG(Electrocardiogram)心电图中的特征作为深度学习的输入,那就要尽可能还原这个信号去除干扰之后原本的样子,这些干扰在专业上又被称为噪声,ECG信号有三大噪声-基线漂移，工频干扰，肌电干扰,而中值滤波就被用于处理基线漂移干扰,为了得到有效的ECG信号这类噪声是必须消除的

一.ECG中的基线漂移介绍

基线漂移通俗的讲就是由于测量电极和人体自身呼吸的原因,会导致测出来的信号会有上下的抖动,这类抖动将难以定位ECG信号的特征点,使测出来的数据可用性降低,接下来放两张图片
1.原始的心电信号

2.去除基线漂移之后的心电信号

可以看到第一张的原始信号中的心电波形"错落有致",这就好像,我们要给一个班级的同学测量身高数据,结果这个班的同学到操场以后有的站在砖头上有的站在土坑里,而我们是以每个同学头顶到地心的距离来表示每个人的身高的,就好像以地心为原点建立了一个坐标轴,而这个时候由于有土坑和砖头的影响,测得的数据不能很好的反应每个同学真实的身高关系(这块应该不难理解,可以稍微想一下)如果以此时的数据作为代表学生身高的依据的话,无论是用这个数据排队列还是排集体照都不会得到一个我们想要的效果,而如果我们要测得更具有参考价值的数据,那我们可以采用填坑垫砖头的方法,让所有的同学都能站到站到一个水平面上,这个时候测得的数据就会变的更具有参考价值,而每一个心电信号就像一个同学,填坑垫砖头的过程就是去除基线漂移的过程,去除基线漂移之后如上面的第二个图,也会让信号看起来’板正’了不少

二.中值滤波

什么叫中值滤波,怎么个滤法呢,简单点解释,就是让每一个数,被包括他周围在内的数中的中位数代替

对列表里的每一个元素进行中值滤波计算,我们使用的是Python中的scipy库的signal.medfilt(滤波数据, 窗大小)

这里有一个简单的例子

import scipy                                                         # 导入scipy模块
window_size = 3
x=[70, 80, 30, 20, 10, 90,  0, 60, 40, 50]             # 定一个列表数组
print(scipy.signal.medfilt(x,window_size))                          #一维中值滤波
# [70 70 30 20 20 10 60 40 50 40]                       #打印滤波后的输出值

接下来我们关注对应的输出值是如何实现的

2.1对于中间数据

输入为[70, 80, 30, 20, 10, 90, 0, 60, 40, 50]
输出为[70, 70, 30, 20, 20, 10, 60, 40, 50, 40]
我们先用数组中的第二个数据举例子原来为80滤波后为值变为70
首先我们关注到scipy.signal.medfilt(x,window_size)
这里的参数窗的大小window_size设置为了3在一维数据中指的是以对应数为中心左右各取(window_size-1)/2个数以这些数作为滤波的基础数据

(window_size-1)/2为了让这个数除尽所以window_size值必须为奇数不然会报错

求这三个数的中位数赋给对应的值
我们再把输入的数据拿出来
x=[70, 80, 30, 20, 10, 90, 0, 60, 40, 50]
对这三个数据进行排序 70,80,30 → 30,70,80
我们再把输出拿出来[70 70 30 20 20 10 60 40 50 40]
对应80这个位置数据变为了70

2.2对应两侧数据

那两侧的数据在原序列中取不到3个数据怎么办呢
例如第一个数据70左边没有数据这不就不够3个数了吗
对应这种情况我们运算的时候会在左边偷偷补一个0
0,[70, 80, 30, 20, 10, 90, 0, 60, 40, 50]
算0,70,80这三个数的中位数
所以对应输出[70 70 30 20 20 10 60 40 50 40]这个数还是70
同理最右边的50
[70, 80, 30, 20, 10, 90, 0, 60, 40, 50],0
对这三个数据进行排序 40,50,0 → 0,40,50
所以对应输出[70 70 30 20 20 10 60 40 50 40]这个数变成了40

三.ECG信号的读取

好了我们现在知道对于一维信号来说中值滤波是什么样子的了
而我们在进行处理时每一个心电信号其实就是被记录为一串数字的数据，每一个数据代表的是当前时刻信号的幅值，每两个数据的采样间隔时间为频率的倒数(这里要是懵住了可以稍微想一下)我们只要把这串数字扔进去,算好每一个数滤波后对应的中值所组成的一维数组,而这个新的一维数组的数据整体表现出来的就是总体的偏移的情况，也就叫做是基线，这块不好理解的话到后面有滤波前后的基线对比帮助大家理解，到这总体什么个思路大致了解了，那心电信号那初始的一维数据哪来呢
这里我们用的是MIT-BIH数据集,在这个数据库里就有我们要分析的心电信号数据

MIT-BIH 是由美国麻省理工学院提供的研究心律失常的数据库

链接: MIT-BHI官方下载地址.
点进去之后下翻点击Download the ZIP file就可以下载MIT-BHI数据集

下载成功之后你将会得到一个叫这个名字的压缩文件

我们建立一个工程将该文件夹解压在工程目录下

之后我们先来简单读取一下心电信号打出来个图看看

# 导入心电信号处理函数库
import wfdb
# 导入python的类matlab绘图函数库
import matplotlib.pyplot as plt
# 读取本地的100号记录，从0到25000，读取模拟信号,通道0
record = wfdb.rdrecord('mit-bih-arrhythmia-database-1.0.0/100', sampfrom=0, sampto=25000, physical=True, channels=[0, ])
# 读取，从第145个数据到第756个数据
ventricular_signal = record.p_signal[144:756]
# 打印标题
plt.title("ventricular signal")
# 打印信号
plt.plot(ventricular_signal)
# 显示图像
plt.show()

wfdb 和 matplotlib 两个库需要自己在电脑的cmd命令提示符操作框内使用pip命令安装安装

命令分别如下

pip install wfdb
pip install matplotlib

其中这里是读取心电信号最关键的一句

record = wfdb.rdrecord('mit-bih-arrhythmia-database-1.0.0/100', sampfrom=0, sampto=25000, physical=False, channels=[0, ])

解释一下这里用到的各个参数
首先第一个参数’mit-bih-arrhythmia-database-1.0.0/100’
这个参数要表达的意思是要读取mit-bih-arrhythmia-database-1.0.0文件夹下的第100号心电数据的参数这个文件夹的名字我们毛都没变直接给粘过来了

sampfrom=0, sampto=25000 这两句的意思是提取对应心电信号中第0到第25000个样本点
physical=False 这里设置成False 则提取出的信号为d_signal数字信号也就是没有转换成心电信号电压幅度之前的原始信号的值可以代表最大的精度而设置成Ture
则会读取出每一个点的模拟信号值也就是对应采样点的mv电压值p_signal
上图对比一下,对比的代码我会放在最后,可以看到两个一模一样的波只不过纵坐标轴的尺度不一样

channels=[0, ]每一次读取的数据都包含两个通道的信号可以设置为channels=[0, ]选择第一个信号,channels=[, 1]选择第二个信号,channels=[0, 1]选择两个信号
好的我们解释清楚我们的每一个函数了直接运行pycharm运行效果如下,我们终于显示出了一个心电信号波形

到这就完成心电信号显示的hello-world

四.MIT-BIH 数据集中值滤波

好了读取完心电信号数据，我们就要把他扔到滤波器里，为了让大家更好的理解里面的细节我们分部分讲解

4.1准备部分代码

import wfdb
import matplotlib.pyplot as plt
from scipy.signal import medfilt
fliter = int(0.8*360)

上面三个导入库函数可以理解，这里的fliter = int(0.8*360)这个东西看起来比较奇怪，这是个什么东西，其实这个东西就说我们要输入到中值滤波函数里的窗的大小，那0.8和360这俩数又代表什么呢，这里要有点耐心，我们一个一个来，先说360，360这个数是因为MIT-BIH数据集每1秒采集360个信号也就是说如果滤波器如果一次对360个数据进行处理那就是作用在了1s的信息上对吧，那为啥要成乘以0.8呢显而易见那就是处理了0.8秒的信息那为啥是0.8秒呢，因为人每次心跳差不多0.8秒一次，每个ECG信号（每一次心跳产生的信号）形状都差不多，也就是说一个周期内代表一个心电信号的所包含的样值点里，取中值应该差距不大，所以，以一个心电信号的窗的大小在整个连续的心电信号图中不断平移对信号进行处理取中值处理的话，那如果信号没有上下抖动(基线漂移)，那对由好几个单独的ECG信号组成的一个连续的信号中的每一个点，加窗，滤波，取得的中值所构成的线（基线，代表信号抖动程度都线），应该很平稳对不对，那如果差距大了说明这个信号抖了,有漂移了对不对，那我们只需要记录下对应点的漂移情况减去这个漂移的量，这不就把信号抻平了嘛

4.2滤波部分代码

record = wfdb.rdrecord('mit-bih-arrhythmia-database-1.0.0/100', sampfrom=0, sampto=25000, physical=True, channels=[0, ])
Original_ECG = record.p_signal[0:2000].flatten()
ECG_baseline = medfilt(Original_ECG, fliter+1)
Filtered_ECG = Original_ECG - ECG_baseline

第一句已经非常详细的解释过了，我们从第二句开始说起
按理说Original_ECG = record.p_signal[0:2000]这不就把一部分心电信号读取出来了吗直接扔进去不就完事了吗，为啥还要后面加一个flatten()函数呢
先简单科普一下flatten()

import numpy as np
a = np.array([[3,7],[9,1]])   # 建立一个2×2的数组
print(a)
#   [[3 7][9 1]]
print(a.flatten())              # 将a展开后输出
# [3 7 9 1]

显而易见这个函数的功能是将多维的数组水平展开成一列成为一个一维数组
我曾经天真的以为Original_ECG = record.p_signal[0:2000]这读出来不就是一个一维的包含2000个数的数组吗直接扔理呗
扔里之后运行程序你会得到这样一个警告

  warnings.warn('kernel_size exceeds volume extent: the volume will be '

然后你的基线就会是一条直线,这就意味着你的基线没有输出
为啥会这样呢这难道不是1维的吗,我要打印看一下然后我就打印了他的形状和数据

好~嘛！真是错怪人家了，一维？，人家可是正经八百的2000维数据，一行就一个数，你一维滤波的窗大小都快到300了能好使就怪了，原来MIT-BIH每个数据读出来的时候单独的数就是是一维
那我们展开一下再看看他的数据的和数据的形状，好的这回是一维两千个数据，舒服了！

好了我们往下到这句了这句里就一个地方要说一下

ECG_baseline = medfilt(Original_ECG, fliter+1)

这个fliter+1为啥要加一还记得窗的大小得是一个奇数吧而我们一开始定义的fliter = int(0.8*360)这必然是个偶数啊，所以务必得加上一个1让它变成奇数才可以
最后一句我们之前讲过了ECG_baseline也就是我们的基线记录了我们信号的便宜情况滤波后的信号Filtered_ECG只需要用Original_ECG(原信号)减去ECG_baseline(基线)就可以得到了

Filtered_ECG = Original_ECG - ECG_baseline

4.3 显示原始信号基线和过滤后的波

plt.figure(figsize=(100, 8))
plt.subplot(3, 1, 1)
plt.ylabel("mv")
plt.title("Original ECG signal")
plt.plot(Original_ECG)#输出原图像
plt.subplot(3, 1, 2)
plt.title("ECG_baseline")
plt.plot(ECG_baseline)#输出基线轮廓
plt.subplot(3, 1, 3)
plt.ylabel("mv")
plt.title("Filtered ECG signal")
plt.plot(Filtered_ECG)#基线滤波结果
plt.show()

最后我们将3条线打印出来

结果如下从上到下依次为原信号,基线,滤波后的信号

好的我们这里就简单实现的对于MIT-BIH数据集的中值滤波

4.3.1处理信号两边的失真

而这个基线两边嫩么高是因为什么呢,是因为是因为我们的滤波器比较宽,在对两边的点做中值处理时,补0比较多取中值时就会有很大的误差,其实有一种很简单的方法可以解决这个问题就是把窗的尺寸变小减小补0取中值的影响,但是结果显而易见抑制到是抑制了但这失真也太严重了吧几乎把信号大部分的特征点给抹平了很显然这是不可取的

fliter = int(0.04*360)

这里我推荐用另一种方法放弃两边需要补0计算中值的特征点只要中间内一部分,这里我们思考一下哪些点在中值滤波的过程中是不需要补0计算的呢,本例中fliter = 0.8*360 = 288 , fliter/2 = 144 也就是在滤波完成之后我们舍弃掉前144和后144个点的值认为这些信号为两边的滤波中产生的噪声直接去除掉

import wfdb
import matplotlib.pyplot as plt
from scipy.signal import medfilt
fliter = int(0.8*360)
Give_up_size = int(fliter/2)
record = wfdb.rdrecord('mit-bih-arrhythmia-database-1.0.0/100', sampfrom=0, sampto=25000, physical=True, channels=[0, ])
Original_ECG = record.p_signal[0:2000].flatten()
print(Original_ECG.shape)
print(Original_ECG)
ECG_baseline = medfilt(Original_ECG, fliter+1)
Filtered_ECG = Original_ECG - ECG_baseline
plt.figure(figsize=(100, 8))
plt.subplot(3, 1, 1)
plt.ylabel("Original ECG signal")
plt.plot(Original_ECG[Give_up_size:-Give_up_size])#输出原图像
plt.subplot(3, 1, 2)
plt.ylabel("ECG_baseline")
plt.plot(ECG_baseline[Give_up_size:-Give_up_size])#输出基线轮廓
plt.subplot(3, 1, 3)
plt.ylabel("Filtered ECG signal")
plt.plot(Filtered_ECG[Give_up_size:-Give_up_size])#基线滤波结果
plt.show()

这里可以看到信号完全没有失真抻平的同时比起最开始的波形两边的波形没有产生失真整理波形的特征得到了很好的保留

4.3.2 处理减去基线之后的幅度值总体偏移

好像已经很不错了,等等好像还有一点瑕疵,仔细观察一下,原波形和,第三个波形,滤波后的波形怎么比原来的波形高了一点,在我认真观察了之后,发现在对心电信号中的模拟信号进行中值滤波的例子里,这里的基线值都是负值,也就是波形减去基线的过程实际上是每一个点都加了一个正数,有点像把每一个波都往上垫,然后垫到同一个高度,虽然信号平了但是总体的幅值高了,那怎么把幅值降下来又不破坏我们好不容易抻平的信号形态呢,我忽然灵机一动,我发现每一个点减去一下总体提升幅度值的平均值不就完事了吗,这不就好像通过几次简单的变换不但抑制了基线漂移保留了信号特征,而且完全没有引入新的噪声

import wfdb
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import medfilt
fliter = int(0.8*360)
Give_up_size = int(fliter/2)
record = wfdb.rdrecord('mit-bih-arrhythmia-database-1.0.0/100', sampfrom=0, sampto=25000, physical=True, channels=[0, ])
Original_ECG = record.p_signal[0:2000].flatten()
print(Original_ECG.shape)
print(Original_ECG)
ECG_baseline = medfilt(Original_ECG, fliter+1)
Filtered_ECG = Original_ECG - ECG_baseline
plt.figure(figsize=(100, 8))
plt.subplot(3, 1, 1)
plt.ylabel("Original ECG signal")
plt.plot(Original_ECG[Give_up_size:-Give_up_size])#输出原图像
plt.subplot(3, 1, 2)
plt.ylabel("ECG_baseline")
plt.plot(ECG_baseline[Give_up_size:-Give_up_size])#输出基线轮廓
plt.subplot(3, 1, 3)
plt.ylabel("Filtered ECG signal")
plt.plot(Filtered_ECG[Give_up_size:-Give_up_size]+np.sum(ECG_baseline[Give_up_size:-Give_up_size])/(2000-fliter))#基线滤波结果
plt.show()

对比上一张观察两次第三张图的Y轴这里明显感觉好多了
接下来我们进行一个简单的优化并打印出被滤波后ECG的基线图
我们用两个变量来代替截取点的上下限

Initial_intercept_point = 0
Final_intercept_point = 2000
Original_ECG = record.p_signal[Initial_intercept_point:Final_intercept_point].flatten()
#代替
Original_ECG = record.p_signal[0:2000].flatten()

目的是可以将总体偏移误差单独提取出来

Totality_Bias = np.sum(ECG_baseline[Give_up_size:-Give_up_size])/(Final_intercept_point-Initial_intercept_point-fliter)

最后打印去除总体偏移之后的信号就换成了

Final_Filtered_ECG = Filtered_ECG[Give_up_size:-Give_up_size]-Totality_Bias

经过上面几步我们就可以只改变Initial_intercept_point,和Final_intercept_point的值后面的就跟着一起变了
然后我们用下面这三行代码对Filtered_ECG_signal求一下基线打出来看看

Filtered_ECG_signal_baseline = medfilt(Filtered_ECG, fliter+1)
plt.plot(Filtered_ECG_signal_baseline[Give_up_size:-Give_up_size])#基线滤波结果
plt.show()

整体代码如下

import wfdb
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import medfilt
fliter = int(0.8*360)
Initial_intercept_point = 0
Final_intercept_point = 2000
Give_up_size = int(fliter/2)
record = wfdb.rdrecord('mit-bih-arrhythmia-database-1.0.0/100', sampfrom=0, sampto=25000, physical=True, channels=[0, ])
Original_ECG = record.p_signal[Initial_intercept_point:Final_intercept_point].flatten()
ECG_baseline = medfilt(Original_ECG, fliter+1)
Totality_Bias = np.sum(ECG_baseline[Give_up_size:-Give_up_size])/(Final_intercept_point-Initial_intercept_point-fliter)
Filtered_ECG = Original_ECG - ECG_baseline
Final_Filtered_ECG = Filtered_ECG[Give_up_size:-Give_up_size]-Totality_Bias
plt.figure(figsize=(100, 8))
plt.subplot(4, 1, 1)
plt.ylabel("Original ECG signal")
plt.plot(Original_ECG[Give_up_size:-Give_up_size])#输出原图像
plt.subplot(4, 1, 2)
plt.ylabel("ECG baseline")
plt.plot(ECG_baseline[Give_up_size:-Give_up_size])#输出基线轮廓
plt.subplot(4, 1, 3)
plt.ylabel("Filtered ECG signal")
plt.plot(Final_Filtered_ECG)#基线滤波结果
plt.subplot(4, 1, 4)
plt.ylabel("Filtered ECG signal baseline")
Filtered_ECG_signal_baseline = medfilt(Filtered_ECG, fliter+1)
plt.plot(Filtered_ECG_signal_baseline[Give_up_size:-Give_up_size])#基线滤波结果
plt.show()

结果如下看纵坐标,可以看到滤波之后的偏移量在0.03mv以内不会为0是因为每种代码都有计算精度上限还是会有一点微小的误差但是总体效果还是不错的

总结

伟大的马克思主义哲学中的唯物辩证法中的否定之否定规律告诉所有的实物是循环往复发展的，都是从是自己变得不是自己又变回更好的自己所以，我最后要做的就是给最第三部分中最开始显示的第一个波形滤波后输出，这里我们将整个ECG的中值滤波过程封装成一个函数直接进行调用

import wfdb
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import medfilt
record = wfdb.rdrecord('mit-bih-arrhythmia-database-1.0.0/100', sampfrom=0, sampto=25000, channels=[0, ])
Original_ECG = record.p_signal[0:1000].flatten()
fliter = int(0.8*360)def medfilt_ecg(original_ecg, window_size):# 确保窗大小变为奇数 window_size = window_size+1 if window_size % 2 == 0 else window_sizegive_up_size = int(fliter / 2)ecg_baseline = medfilt(Original_ECG, window_size)totality_bias = np.sum(ecg_baseline[give_up_size:-give_up_size])/(len(original_ecg)-2*give_up_size)filtered_ecg = original_ecg - ecg_baselinefinal_filtered_ecg = filtered_ecg[give_up_size:-give_up_size]-totality_biasreturn final_filtered_ecgplt.title("ventricular signal")
plt.ylabel('mv')
plt.plot(medfilt_ecg(Original_ECG, fliter))
plt.show()

physical=True 为默认值直接省略
用len(original_ecg)表示数据长度

结果如下

附件模拟信号数字信号对比代码

import wfdb
import matplotlib.pyplot as plt
# 读取数字信号
record = wfdb.rdrecord('mit-bih-arrhythmia-database-1.0.0/100', sampfrom=0, sampto=25000, physical=False, channels=[0, ])
ventricular_d_signal = record.d_signal[0:1000]
plt.figure(figsize=(20, 4))
plt.subplot(1, 2, 1)
plt.title("ventricular d_signal")
plt.plot(ventricular_d_signal)
# 读取模拟信号
record = wfdb.rdrecord('mit-bih-arrhythmia-database-1.0.0/100', sampfrom=0, sampto=25000, physical=True, channels=[0, ])
ventricular_p_signal = record.p_signal[0:1000]
plt.subplot(1, 2, 2)
plt.title("ventricular p_signal")
plt.ylabel("mv")
plt.plot(ventricular_p_signal)
plt.show()

后记

第一次写信号处理类文章诚惶诚恐,这一篇文章其实就浓缩了我对这个知识点的学习经历,中值滤波窗大小的设置这里也是加入了我自己的理解,甚至4.3.1和4.3.2这两个方法是我写到一半边写边思考忽然想到的,然后去实验发现非常好使就写上了,想到这两个东西的时候真的没有看任何的文献,所以,虽然简单,但是也感觉非常开心,学习路漫漫,如果有什么不对的地方恳请大家批评指正,最后感谢我同学用催更表达对我的认可哈哈哈哈