CNN心电识别笔记
参考CSDN文献:https://blog.csdn.net/qq_43750573/article/details/105930152
一、心电信号的简介
-
心电信号的特征
1. 微弱性:幅值在10uV-5mⅤ,属于低幅值信号 2. 不稳定性:易受干扰,大量噪声(卷积神经网络要解决降噪问题) 3. 低频性:频率范围在0.05-100Hz,主要能量集中分布在0.5-40Hz
-
心电波形(一个心拍)
-
降噪手段
1. 经典数字滤波器。根据频率范围降噪 对于基线漂移使用高通滤波器去噪 对于肌电干扰使用低通滤波器去噪 对于工频干扰使用带通滤波器去噪 2. 基于小波变换的阈值去噪。 根据信号和噪声的频率在不同尺度上的分布,先对信号进行小波变换,再根据阈值对各层小波系数进行处理,最后重构信号进行去噪。
二、心率失常数据库
-
MIT-BIH心律失常数据库
下载地址:https://www.physionet.org/content/mitdb/1.0.0
-
数据读取
MIT-BIH 使用了自定义的 format212 格式,需要进行解码
import wfdb """ 读取心电信号的数据 data_path:路径。不必加拓展名,自动补齐 sampfrom:读取心电信号的起始位置 sampto:读取心电信号的结束位置 channel_names:列表值。读取导联的名字["MLII"] channels:设置读取第几个心电信号 return:返回的是一个元组 (信号,记录描述符) 信号维度 [650000,1],如果不指定导联,得到 [65000,2] 的结果 """ # 仅仅读取“MLII”导联的信号 record = wfdb.rdrecord('data_path', channel_names=['MLII']) """ 读取注解文件 extension:文件拓展名。string “atr” 其他参数同上 return:一个注释对象 """ annotation = wfdb.rdann("data_path",extensiom = "atr") # 返回每一个心拍的 R 波的尖锋位置的信号点,与心电信号对应。返回值是一个 ndarray 数组 annotation.sample # 返回每一个心拍的类型 N,L,R 等等。返回值是一个数组 annotation.symbol # 返回标注的数量 annotation.ann_len # 返回标注的文件名 annotation.record_name # 查看所有心拍的类型 wfdb.show_ann_labels()
三、心电信号去噪
-
噪声分类
1. 工频干扰 电磁干扰,幅值低,噪声频率为50Hz左右,其波形很像一个正弦信号,会影响 P 波和 T 波的检测。 2. 肌电干扰 人体运动肌肉不自主颤抖造成,这种干扰无规律可言,波形形态会急速变化,频率很高,并且分布很广,范围在0-2000Hz内,能量集中在30-300Hz内,持续时间一般为50ms 3. 基线漂移 属于低频干扰,频率分布在0.15-0.3Hz内,幅度和频率都会时刻变化着。心电信号中的PR波段和ST波段非常容易受到影响产生失真。
-
心电信号预处理-小波变换
-
步骤:
- 选择一个小波基函数,对心电信号进行尺度分解,得到各尺度上的小波系数
- 对各尺度上的小波系数进行阈值处理。幅值较大的小波系数为信号,较小的为噪声。噪声小波系数调用阈值函数处理
- 重构信号
-
阈值函数
-
硬阈值函数。系数绝对值大于阈值,保持不变,小于阈值令其为零
-
软阈值函数。系数绝对值大于阈值,减去 λ ,小于阈值令其为零
-
λ 计算公式
median为中位数,N为信号长度
-
-
-
代码实现
""" 时频分析,小波分解 wavelet:小波基 level:尺度(要变换多少层) return:一个列表值,[低频系数,高频系数] 低频系数是一个 array 数组 高频系数是一个元组(水平高频系数,垂直高频系数,对角线高频系数),有几层小波分解就有几个元组 """ coeffs = pywt.wavedec(data,wavelet= "db5",level = 9) cA9, cD9, cD8, cD7, cD6, cD5, cD4, cD3, cD2, cD1 = coeffs """ 阈值处理,噪声滤波 pywt.threshold(data,value,mode,subtitute) value:阈值 mode:模式,有四种。soft软阈值函数;hard硬阈值函数;greater大于等于阈值的值不变化其余的置为subtitute;less将大于阈值的值替换为subtitute,其余的值不变 subtitute:替换值 """ # 计算阈值 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")
四、构建卷积神经网络
-
构建数据集
1. 构建心拍数据 通过 annotation.sample 获取 R 波尖峰位置 指定一个 R 波尖峰位置的前99个信号点与后200个信号点,一共300个信号点为一个心拍 2. 构建标签数据 通过 annotation.symbol 获取每一个心拍的类型 仅选取 ['N', 'A', 'V', 'L', 'R'] 五种心拍类型,舍去不属于这五种心拍类型的心拍数据 将标签数据集转化为 one-hot 编码 3. 将数据集和标签集拼接,打乱 4. 划分出训练数据集与测试数据集
-
构建卷积神经网络
样本的数据维度为 [None,300] 标签数据维度为 [None,5] 卷积 激活 relu函数激活 池化 max_pooling,mean_pooling
-
全连接层
计算交叉熵损失
梯度下降降低损失
求得准确率
-
收集变量用于可视化操作
五、整体代码实现
import wfdb
import pywt
import numpy as np
from matplotlib import pyplot as plt
import tensorflow as tf
FLAGS = tf.app.flags.FLAGS
tf.app.flags.DEFINE_string("data_path", "./mit-bih-arrhythmia-database-1.0.0/", "心电数据的路径")
def de_noise(data):
"""
小波去噪函数
:param data:原始数据
:return: rdata 去噪后的数据
"""
# 1、时频分析,小波去噪
coeffs = pywt.wavedec(data,wavelet = "db5",level = 9)
cA9, cD9, cD8, cD7, cD6, cD5, cD4, cD3, cD2, cD1 = coeffs
# 2、阈值处理,噪声滤波
threshold = (np.median(np.abs(cD1)) / 0.6745) * (np.sqrt(2 * np.log(len(cD1)))) # 计算阈值
cD1.fill(0)
cD2.fill(0) # 高频信号置零
# 将其他中低频信号软阈值滤波
for i in range(1,len(coeffs) -2 ):
coeffs[i] = pywt.threshold(coeffs[i],threshold,mode="soft")
# 3、信号进行小波重构
rdata = pywt.waverec(coeffs=coeffs,wavelet="db5")
return rdata
def process_data(num,X_data,Y_data):
"""
获取数据并进行小波变换,截取心拍
:param num: 第几条数据
:param X_data: 样本数据集
:param Y_data: 标签数据集
:return:None
"""
# 1、获取原始信号数据.signal_data为 [650000,1] , 需要展开
rec = wfdb.rdrecord(FLAGS.data_path + num, channel_names=["MLII"])
signal_data = rec.p_signal
transpose_data = signal_data.flatten() #展开
# 2、原始数据进行小波分解
rdata = de_noise(transpose_data)
# 3、获取心电信号中 R 波尖峰位置以及对应的标签数据
annotation = wfdb.rdann(FLAGS.data_path + str(num), extension= "atr")
R_location = annotation.sample # 返回每个心拍 R 波尖峰的位置,numpy 数组
R_class = annotation.symbol # 返回每个心拍的类型,list 列表
# 4、去除前后不稳定的数据
start = 10
end = 5
i = start # 从第10个心拍开始,到最后5个心拍结束
j = len(R_class) - end
# 5、截取符合条件的心拍信号,标签数据
class_list = ['N', 'A', 'V', 'L', 'R'] # 仅对此五种心拍类型进行训练
while i < j:
try:
# 获取当前心拍在心拍类型中的位置,即属于哪种心拍。如果索引失败即此心拍不属于五种心拍,弹出错误,跳过此心拍
label = class_list.index(R_class[i])
local = R_location[i] # 当前心拍的 R 波尖峰位置
train_data = rdata[local - 99 : local + 201] # 截取出一个完整的心拍,R 波尖峰前99个到后200个
X_data.append(train_data) # 将一个完整的心拍样本加入到数据集中
Y_data.append(label) # 将心拍的标签加入到标签集中
i += 1
except:
i += 1
return None
def shuffle_data(X_data,Y_data):
"""
将数据集与标签机一一对应后合并,合并后打乱
:param X_data:数据集
:param Y_data:标签集
:return: data 打乱后的数据集
"""
dataSet = np.array(X_data).reshape(-1,300)
labelSet = np.array(Y_data).reshape(-1,1)
# 数组拼接
stitch_data = np.hstack((dataSet,labelSet))
# 打乱数组
np.random.shuffle(stitch_data)
return stitch_data
def weight_variables(shape):
"""
初始化权重
:param shape:
:return: w 指定形状的权重变量
"""
w = tf.Variable(tf.random_normal(shape=shape,mean=0.0,stddev=1.0))
return w
def bias_variables(shape):
"""
初始化偏置
:param shape:
:return: 指定形状的偏置变量
"""
b = tf.Variable(tf.constant(0.0,shape= shape))
return b
def CNN_model():
"""
自定义的卷积神经网络
:return: x,y_true,y_predict 变量
"""
# 1、准备占位符
with tf.variable_scope("data"):
x = tf.placeholder(tf.float32,[None,300])
y_true = tf.placeholder(tf.int32,[None,5])
# 2、一层卷积.filter 1*5*1 32个 strides = 1
with tf.variable_scope("conv1"):
# 随机初始化权重,偏置
w_conv1 = weight_variables([1,5,1,32])
b_conv1 = bias_variables([32])
# 对 x 进行形状变化 [None,300] --> [None,1,300,1]
x_reshape = tf.reshape(x,[-1,1,300,1])
# filter卷积 [None,1,300,1] --> [None,1,300,32]
convolution1 = tf.nn.conv2d(x_reshape,w_conv1,strides=[1,1,1,1],padding="SAME")
# relu 函数激活
x_relu1 = tf.nn.relu(convolution1 + b_conv1)
# 池化 [None,1,300,32] --> [None,1,150,32]
x_pool1 = tf.nn.max_pool(x_relu1,ksize=[1,1,2,1],strides=[1,1,2,1],padding="SAME")
# 3、二层卷积 1*5*32 64个
with tf.variable_scope("conv2"):
# 初始化权重偏置
w_conv2 = weight_variables([1,5,32,64])
b_conv2 = weight_variables([64])
# filter卷积 [None,1,150,32] --> [None,1,150,64]
convolution2 = tf.nn.conv2d(x_pool1,w_conv2,strides=[1,1,1,1],padding="SAME")
# 激活
x_relu2 = tf.nn.relu(convolution2 + b_conv2)
# 池化 1*2 [None,1,150,64] --> [None,1,75,64]
x_pool2 = tf.nn.avg_pool(x_relu2,ksize=[1,1,2,1],strides=[1,1,2,1],padding="SAME")
# 4、三层卷积 1*5*64 32个
with tf.variable_scope("conv3"):
w_conv3 = weight_variables([1,5,64,32])
b_conv3 = bias_variables([32])
# 卷积、激活、池化 [None,1,75,64] --> [None,1,75,32]
convolution3 = tf.nn.conv2d(x_pool2,w_conv3,strides=[1,1,1,1],padding="SAME")
x_relu3 = tf.nn.relu(convolution3 + b_conv3)
# 1*3 [None,1,75,32] --> [None,1,25,32]
x_pool3 = tf.nn.max_pool(x_relu3,ksize=[1,1,3,1],strides=[1,1,3,1],padding="SAME")
# 5、全连接层 [None,1,25,32] --> [None,1*25*32] * [1*25*32,5] + [5] = [None,5]
with tf.variable_scope("full_connect"):
# 初始化权重和偏置
w_fc = weight_variables([1*25*32,5])
b_fc = bias_variables([5])
# 修改形状
x_fc_reshape = tf.reshape(x_pool3,[-1,1*25*32])
# 矩阵运算得出预测结果
y_predict = tf.matmul(x_fc_reshape,w_fc) + b_fc
# 高维变量收集
tf.summary.histogram("weights", w_fc)
tf.summary.histogram("biases",b_fc)
return x,y_true,y_predict
def main():
"""
定义主函数进行卷积分类
:return:
"""
# 1、加载数据集,将心拍数据以及标签数据加载到列表中
data_list = []
label_list = []
number_list = ['100', '101', '103', '105', '106', '107', '108', '109', '111', '112', '113', '114', '115',
'116', '117', '119', '121', '122', '123', '124', '200', '201', '202', '203', '205', '208',
'210', '212', '213', '214', '215', '217', '219', '220', '221', '222', '223', '228', '230',
'231', '232', '233', '234']
for num in number_list:
process_data(num,data_list,label_list)
# 2、打乱数据集,分离实际的数据集、标签集
sf_data = shuffle_data(data_list,label_list)
sf_data_list = sf_data[:,:300] # 前 300 列为数据
sf_label_list = sf_data[:,300] # 最后一列为标签
one_hot_label = tf.one_hot(sf_label_list,depth=5) # 转换为独热编码
# 3、定义模型得出输出
x,y_true,y_predict = CNN_model()
# 4、交叉熵计算
with tf.variable_scope("soft_cross"):
# 求出平均交叉熵损失
loss = tf.reduce_mean(tf.nn.softmax_cross_entropy_with_logits(labels=y_true,logits=y_predict))
# 5、梯队下降降低损失
with tf.variable_scope("optimizer"):
train_op = tf.train.GradientDescentOptimizer(0.001).minimize(loss)
# 6、计算准确率
with tf.variable_scope("acc"):
equal_list = tf.equal(tf.argmax(y_true,1),tf.argmax(y_predict,1))
accuracy = tf.reduce_mean(tf.cast(equal_list,tf.float32))
# 7、收集变量
# 单个数字值收集
tf.summary.scalar("losses",loss)
tf.summary.scalar("acc",accuracy)
# 定义一个初始化的变量 op
init_op = tf.global_variables_initializer()
# 定义一个合并变量的op
merged = tf.summary.merge_all()
# 创建一个 saver
saver = tf.train.Saver()
# 开启会话运行
with tf.Session() as sess:
# 初始化变量 op
sess.run(init_op)
# 建立 events 文件并写入
filewriter = tf.summary.FileWriter("./tmp/summary/test/", graph=sess.graph)
one_hot_label = sess.run(one_hot_label)
# 循环训练
for i in range(1000):
# 获得真实的特征值,标签值,每次投喂90条数据
true_x = sf_data_list[i * 90 : (i + 1) * 90,:]
true_y = one_hot_label[i * 90 : (i + 1) * 90,:]
# 运行 train_op 进行训练
sess.run(train_op,feed_dict={ x : true_x , y_true : true_y})
# 写入每一步的训练值
summary = sess.run(merged,feed_dict={ x : true_x , y_true : true_y})
filewriter.add_summary(summary,i)
print("第 %d 步训练,准确率为:%f " % (i,sess.run(accuracy,feed_dict={x : true_x , y_true : true_y})))
# 保存模型
saver.save(sess,"./tmp/ckpt/fc_model")
return None
if __name__ == '__main__':
main()