python MIT-BIH 心拍分类实现
前言
这可能是我自认为写的最难的一篇文章,对于所有做心电信号分类的同学来说这都是一大关以下所有的代码都是自己手敲的,这整个过程很复杂,按照我一贯的思路,我会写的尽量详细一是日后方便复习,二是自己看论文也好看文章也罢,好多都是囫囵半片,看了跟没看似的,所以自己一定要写的详细,一是自己对理论的一种梳理和理解,二是给相关的研究者提供清晰的思路,三是希望能给学弟学妹们做个榜样接下来进入正题
一.本文的主要内容和研究思路
首先不管对于机器学习还是深度学习,不管你最后要做的是几分类任务,你都需要把每一个心电信号提取出将每一种信号打上标签分类放在一起,而研究心电信号最常用的数据库就是MIT-BIH数据库,python也有专门的wfdb函数库对数据集进行读取,那将单独的心电信号提取出来就是第一步,提取的时候我们根据通过函数读取的标签再将每一类分到一起,打包作为数据集,打包之后准备一个高效的读取方式能让他在将来的分类时被调用为了实现上面的步骤就需要一个庞大的过程我们从头开始，如果你有研究基础直接移步到第四步即可
二.ECG信号的读取
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.d_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数据集中有心拍种类如下表
心拍就指的是在心电信号中包含有用信息的一个数据点
| Label_story | symbol | description | 中文描述 |
|---|---|---|---|
| 0 | Not an actual annotation | ||
| 1 | N | Normal beat | 正常心搏 |
| 2 | L | Left bundle branch block beat | 左束支传导阻滞 |
| 3 | R | Right bundle branch block beat | 右束支传导阻滞 |
| 4 | a | Aberrated atrial premature beat | 异常心房早搏 |
| 5 | V | Premature ventricular contraction | 室性早搏 |
| 6 | F | Fusion of ventricular and normal beat | 心室融和心搏 |
| 7 | J | Nodal (junctional) premature beat | 交界性早搏 |
| 8 | A | Atrial premature contraction | 房性早搏 |
| 9 | S | Premature or ectopic supraventricular beat | 早搏或室上性一尾心搏 |
| 10 | E | Ventricular escape beat | 室性逸搏 |
| 11 | j | Nodal (junctional) escape beat | 交界性逸搏 |
| 12 | / | Paced beat | 起搏心拍 |
| 13 | Q | Unclassifiable beat | 无法分类的节拍 |
| 14 | ~ | Signal quality change | 信号质量变化 |
| 16 | | | Isolated QRS-like artifact | |
| 18 | s | ST change | St改变 |
| 19 | T | T-wave change | T波改变 |
| 20 | * | Systole | 心脏收缩 |
| 21 | D | Diastole | 心脏舒张 |
| 22 | " | Comment annotation | 注释注释 |
| 23 | = | Measurement annotation | 测量注释 |
| 24 | p | P-wave peak | P波峰值 |
| 25 | B | Left or right bundle branch block | 传导阻滞心搏 |
| 26 | ^ | Non-conducted pacer spike | 起搏器伪迹 |
| 27 | t | T-wave peak | T波峰值 |
| 28 | + | Rhythm change | 节律变化 |
| 29 | u | U-wave peak | U波峰值 |
| 30 | ? | Learning | |
| 31 | ! | Ventricular flutter wave | 心室颤动 |
| 32 | [ | Start of ventricular flutter/fibrillation | 开始心室扑动/颤动 |
| 33 | ] | End of ventricular flutter/fibrillation | 结束心室扑动颤动 |
| 34 | e | Atrial escape beat | 房性逸搏 |
| 35 | n | Supraventricular escape beat | 室上性逸搏 |
| 36 | @ | Link to external data (aux_note contains URL) | |
| 37 | x | Non-conducted P-wave (blocked APB) | 未下传的P波 |
| 38 | f | Fusion of paced and normal beat | 起搏融合心跳 |
| 39 | ( | Waveform onset PQ junction（begin of QRS） | QRS波开始 |
| 40 | ) | Waveform end IPT junction(J point,end of QRS) | QRS波结束 |
| 41 | r | R-on-T premature ventricular contraction | R波落在T上的室性早搏 |
wfdb函数库中有专门的函数来提取出的信号中有哪些心拍点，并读取出这些点在信号中的位置
import wfdb
signal_annotation = wfdb.rdann('mit-bih-arrhythmia-database-1.0.0/106', "atr", sampfrom=0, sampto=1000)
# 打印标注信息
print("symbol: " + str(signal_annotation.symbol))
print("sample: " + str(signal_annotation.sample))
#运行结果如下
#symbol: ['~', '+', 'N', 'N']
#sample: [ 83 229 351 724]
"atr"是标注文件的文件类型固定不需要更改，其他三个参数用法和wfdb.rdrecord()函数相同
接下来将这些点在信号中用散点图的方式打印出来
import wfdb
import matplotlib.pyplot as plt
record = wfdb.rdrecord('mit-bih-arrhythmia-database-1.0.0/106', sampfrom=0, sampto=1000, physical=True, channels=[0, ])
signal_annotation = wfdb.rdann('mit-bih-arrhythmia-database-1.0.0/106', "atr", sampfrom=0, sampto=1000)
# 打印标注信息
ECG = record.p_signal
plt.plot(ECG)
#按坐标在散点图上绘点
for index in signal_annotation.sample:
plt.scatter(index, ECG[index], marker='*',s=200)
plt.show()
打印出来的结果如下
首先我们看一下一个心电信号中最高点是R点，要提取一个心电信号我们就要找到R点的位置左右截取，可是可以看到在MIT-BIH给的所有心拍标签中有代表R点也有不代表R点的所以就要把代表R点的心拍标签的位置提取出来才好进行下一步截取操作
那那些心拍标签代表的是R点的标签呢，就是上面列表中标黄的19类
ECG_R_list = np.array(['N', 'f', 'e', '/', 'j', 'n', 'B',
'L', 'R', 'S', 'A', 'J', 'a', 'V',
'E', 'r', 'F', 'Q', '?'])
而根据AAMI标准又可以分成5类
AAMI_MIT = {'N': 'Nfe/jnBLR',#将19类信号分为五大类
'S': 'SAJa',
'V': 'VEr',
'F': 'F',
'Q': 'Q?'}
所有的患者编号文件名如下
file_name =['100', '101', '102', '103', '104', '105', '106', '107',
'108', '109', '111', '112', '113', '114', '115', '116',
'117', '118', '119', '121', '122', '123', '124', '200',
'201', '202', '203', '205', '207', '208', '209', '210',
'212', '213', '214', '215', '217', '219', '220', '221',
'222', '223', '228', '230', '231', '232', '233', '234']
四.分类读取
先介绍点函数吧最后的程序里用到的
numpy.isin()函数
numpy.isin()会生成一个布尔型的数组
import numpy as np
A = np.array(['N', 'N', 'F', 'C', 'N', 'A', 'N'])
B = np.array(['N', 'C', 'Q'])
C = np.isin(A, B)
print(f'C={C}')
D = A[C]
print(f'D={D}')
#C=[ True True False True True False True]
#D=['N' 'N' 'C' 'N' 'N']
例如这里生成了一个和A长度一样的布尔型数组，每一位的值是True还是False取决于A对应位置的元素在不在B里，并且这个布尔型数组可以作为索引扔到原数组里这样就可以把对应的位置的元素提取出来组成一个新的集合
set（）函数
set() 函数创建一个无序不重复元素集，可进行关系测试，删除重复数据，还可以计算交集、差集、并集等
import numpy as np
A = np.array(['N', 'N', 'F', 'C', 'N', 'A', 'N'])
B = np.array(['N', 'C', 'Q'])
C = np.isin(A, B)
D = A[C]
E = set(D)
print(f'C={C}')
print(f'D={D}')
print(f'E={E}')
#C=[ True True False True True False True]
#D=['N' 'N' 'C' 'N' 'N']
#E={'C', 'N'}
这样我们就知道在A中有哪几种元素在B里，后面的分类需要用到
enumerate（）
enumerate() 函数用于将一个可遍历的数据对象(如列表、元组或字符串)组合为一个索引序列，同时列出数据和数据下标，一般用在 for 循环当中
seasons = ['Spring', 'Summer', 'Fall', 'Winter']
print(list(enumerate(seasons)))
#[(0, 'Spring'), (1, 'Summer'), (2, 'Fall'), (3, 'Winter')]
链接: 菜鸟教程enumerate
Python List extend()
extend() 函数用于在列表末尾一次性追加另一个序列中的多个值（用新列表扩展原来的列表）
aList = [123, 'xyz', 'zara', 'abc', 123];
bList = [2009, 'manni'];
aList.extend(bList)
print ("Extended List : ", aList);
#Extended List : [123, 'xyz', 'zara', 'abc', 123, 2009, 'manni']
链接: 菜鸟教程extend
完整分类代码
import wfdb
import numpy as np
import csv
AAMI_MIT = {'N': 'Nfe/jnBLR',#将19类信号分为五大类
'S': 'SAJa',
'V': 'VEr',
'F': 'F',
'Q': 'Q?'}
# 5分类存储字典
AAMI = {'N': [],'S': [],'V': [],'F': [], 'Q': []}
file_name =['100', '101', '102', '103', '104', '105', '106', '107',
'108', '109', '111', '112', '113', '114', '115', '116',
'117', '118', '119', '121', '122', '123', '124', '200',
'201', '202', '203', '205', '207', '208', '209', '210',
'212', '213', '214', '215', '217', '219', '220', '221',
'222', '223', '228', '230', '231', '232', '233', '234']
ECG = {'N': [], 'f': [], 'e': [], '/': [], 'j': [], 'n': [], 'B': [],
'L': [], 'R': [], 'S': [], 'A': [], 'J': [], 'a': [], 'V': [],
'E': [], 'r': [], 'F': [], 'Q': [], '?': []}
for F_name in file_name:
signal_annotation = wfdb.rdann(f'mit-bih-arrhythmia-database-1.0.0/{F_name}', "atr", sampfrom=0, sampto=650000)
record = wfdb.rdrecord(f'mit-bih-arrhythmia-database-1.0.0/{F_name}', sampfrom=0, sampto=650000, physical=True, channels=[0, ])
# 删除非R点的标签
ECG_R_list = np.array(['N', 'f', 'e', '/', 'j', 'n', 'B',
'L', 'R', 'S', 'A', 'J', 'a', 'V',
'E', 'r', 'F', 'Q', '?'])
# 获取表示R点的心拍类型的索引
Index = np.isin(signal_annotation.symbol, ECG_R_list)
# 将标签从list列表类型转化为数组类型为了用下面的索引值提取
Label = np.array(signal_annotation.symbol)
# 提取表示为R点的心拍标签
Label = Label[Index]
# 提取表示为R点的坐标
Sample = signal_annotation.sample[Index]
# 获取心拍种类
Label_kind = list(set(Label))
# 读取每一种R点在信号中的位置
for k in Label_kind:
index = [i for i, x in enumerate(Label) if x == k]
Signal_index = Sample[index]
length = len(record.p_signal)
#print(Signal_index)
#截取
for site in Signal_index:
if 130 < site < length-130:
ECG_signal = record.p_signal.flatten().tolist()[site-130:site+130]
#print(ECG_signal)
ECG[str(k)].append(ECG_signal)
# 打印种类
for key, value in ECG.items():
print(f'{key} = {len(value)}')
for ECG_key, ECG_value in ECG.items():
for AAMI_MIT_key, AAMI_MIT_value in AAMI_MIT.items():
if ECG_key in AAMI_MIT_value:
AAMI[AAMI_MIT_key].extend(ECG_value)
# 5分类并保存到CSV格式的文件里
for key, value in AAMI.items():
with open(f'{key}.csv', 'w',newline='\n') as f:
writer = csv.writer(f)
# 将列表的每条数据依次写入csv文件， 并以逗号分隔
# 传入的数据为列表中嵌套列表或元组，每一个列表或元组为每一行的数据
writer.writerows(value)
运行的结果如下
运行时间因电脑而异但是要基上运行10多分钟毕竟要处理40个通道的40×650000个点
首先你会在串口处打印出每一种类型心拍点的数量
可以看到有四种是0所以有的论文里也会说有15类
N = 75021
f = 982
e = 16
/ = 7025
j = 229
n = 0
B = 0
L = 8071
R = 7256
S = 2
A = 2546
J = 83
a = 150
V = 7129
E = 106
r = 0
F = 802
Q = 33
? = 0
其次你会看见你程序文件夹的目录下生成了5个CSV文件
点开一个感受一下看一下每一行就是一个心电信号的数据
五.信号读取
读取之后在AAMI字典生成五个键值对,每一个键值对中的值是list of list(列表的列表)
包含着每一种心拍类型的心电信号数据的列表
读取只需要10秒左右就可以
map（）函数
map() 会根据提供的函数对指定序列做映射。
第一个参数 function 以参数序列中的每一个元素调用 function 函数，返回包含每次 function 函数返回值的新列表
def square(x): # 计算平方数
return x ** 2
print(list(map(square, [1,2,3,4,5])))# 计算列表各个元素的平方
#[1, 4, 9, 16, 25]
链接: 菜鸟教程map
好了最后再上读取的代码实际使用中用他
import csv
AAMI_key = {'N','S','V','F','Q'}
AAMI = dict()
def Tolist(x):
return list(map(float,x))
for key in AAMI_key:
with open(f'{key}.csv', 'r') as f:
reader = csv.reader(f)
AAMI[key] = list(map(Tolist,list(reader)))
最后转换成数组看一眼形状
import csv
import numpy as np
AAMI_key = {'N','S','V','F','Q'}
AAMI = dict()
def Tolist(x):
return list(map(float,x))
for key in AAMI_key:
with open(f'{key}.csv', 'r') as f:
reader = csv.reader(f)
AAMI[key] = list(map(Tolist,list(reader)))
for i in AAMI_key:
print(f'{i}={np.array(AAMI[i]).shape}')
结果
V=(7235, 260)
Q=(33, 260)
S=(2781, 260)
N=(98600, 260)
F=(802, 260)
结束语
这里只是实现了分类但是没有滤波,可以分类之前滤也可以分类之后滤，这个是我下一步要做的，想看懂所有代码确实还是需要一定的python基础,代码有看不懂的地方去搜一下函数的功能就好了,基本我都贴上了，也可以作为黑盒直接使用不关心内部实现只调用最后的字典即可