基于Python的心电信号检测与处理
Posted 鲁棒最小二乘支持向量机
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了基于Python的心电信号检测与处理相关的知识，希望对你有一定的参考价值。
心电信号的特征提取、分析与处理
1.生物医学信号的特征提取与分析方法
2.生物医学信号的滤波方法
数据来源：
MIT-BIH数据库（可从以下数据中任选两组进行实验）
给出4组不同病例的心电信号数据，分别命名为“100-2-3”，“105-2-3”，“109-2-3”，“111-2-3”，每组数据以“.mat”形式存储。
每组数据长度N=2048，采样率fs=360Hz(硬件采集心电信号时，每秒采集360个点。注意设计FIR,IIR时可能用到该采样率。). 即2048点对应时间约为5.69s
任务内容：
（1）谱分析： 取两段心电信号数据（不同病例），对该数据进行频谱分析（幅度谱、相位谱、功率谱）；
（2）相关分析：分别计算两段心电信号的均值、方差、自相关函数与互相关函数；分析两段信号的相干函数曲线，对于相关函数进行循环相关函数与线性相关函数的对比；
（3）讨论：根据（1）、（2）得出的结果，分析并总结心电信号的特征（例如，心电信号能量主要分布范围、自相关函数与互相关函数有没有周期特性等）
（4）数字滤波器设计：取一段心电信号，添加白噪声。分别作出原始信号、加噪信号的图；作出原始信号、加噪信号的自相关曲线，分析估计心电信号的准周期；取一段心电信号，添加高频噪声（1k-2khz），如高频正弦信号，频率自己选择。结合（1）中得出的结论，即ECG的主要能量分布结果，设计数字滤波器（IIR或FIR），去除高频噪声。（也可直接设计数字滤波器去除基线漂移） 要求：对滤波器的截止频率的设置要给出说明。
（5）维纳滤波器去除工频干扰：取一段心电信号，添加频率为50Hz的高斯白噪声（工频干扰）。设计维纳滤波器，分析滤波结果。
Python代码实现
原始心电信号数据代码
#coding:UTF-8
import scipy.io as scio
import numpy as np
import matplotlib.pyplot as plt
# %matplotlib inline
#设置绘图大小
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus']=False
dataFile_100 = './100-2-3.mat'
dataFile_105 = './105-2-3.mat'
dataFile_109 = './109-2-3.mat'
dataFile_111 = './111-2-3.mat'
data_100 = scio.loadmat(dataFile_100)
data_105 = scio.loadmat(dataFile_105)
data_109 = scio.loadmat(dataFile_109)
data_111 = scio.loadmat(dataFile_111)
Y_100 = data_100['y']
Y_105 = data_105['y']
Y_109 = data_109['y']
Y_111 = data_111['y']
Fs = 360
N = 2048
frq = np.arange(N) #频率数2048个数
half_x = frq[range(int(N/2))] #取一半区间
print ("心电信号数据:",Y_100)
print ("数据长度:",len(Y_100))
# 原始数据
plt.figure(figsize=(7,5))
plt.subplot(221)
plt.title('100-2-3原始心电信号数据')
plt.plot(range(0,N), Y_100)
plt.xlabel("N")
plt.ylabel("心电信号")
plt.subplot(222)
plt.title('105-2-3原始心电信号数据')
plt.plot(range(0,N), Y_105)
plt.xlabel("N")
plt.ylabel("心电信号")
plt.subplot(223)
plt.title('109-2-3原始心电信号数据')
plt.plot(range(0,N), Y_109)
plt.xlabel("N")
plt.ylabel("心电信号")
plt.subplot(224)
plt.title('111-2-3原始心电信号数据')
plt.plot(range(0,N), Y_111)
plt.xlabel("N")
plt.ylabel("心电信号")
plt.show()
原始心电信号数据结果
部分原始心电信号数据代码
# 部分原始数据
plt.figure(figsize=(7,5))
plt.subplot(221)
plt.title('100-2-3原始心电信号部分数据')
plt.plot(frq[0:100], Y_100[0:100])
plt.xlabel("N")
plt.ylabel("心电信号")
plt.subplot(222)
plt.title('105-2-3原始心电信号部分数据')
plt.plot(frq[0:100], Y_105[0:100])
plt.xlabel("N")
plt.ylabel("心电信号")
plt.subplot(223)
plt.title('109-2-3原始心电信号部分数据')
plt.plot(frq[0:100], Y_109[0:100])
plt.xlabel("N")
plt.ylabel("心电信号")
plt.subplot(224)
plt.title('111-2-3原始心电信号部分数据')
plt.plot(frq[0:100], Y_111[0:100])
plt.xlabel("N")
plt.ylabel("心电信号")
plt.show()
部分原始心电信号数据结果
谱分析代码
# 谱分析
mf_100 = np.fft.fft2(Y_100) # 进行频谱变换（傅里叶变换）
mf_105 = np.fft.fft2(Y_105)
mag_100 = np.abs(mf_100) # 取复数的绝对值，即复数的模(双边频谱)
mag_105 = np.abs(mf_105)
gui_mag_100 = mag_100/N # 归一化处理（双边频谱）
gui_mag_105 = mag_105/N
gui_half_mag_100 = gui_mag_100[range(int(N/2))] #由于对称性，只取一半区间（单边频谱）
gui_half_mag_105 = gui_mag_105[range(int(N/2))]
f_100 = 20*np.log10(mag_100[:N//2]) # 通过 20*np.log10()将其转换为以db单位的值
f_105 = 20*np.log10(mag_105[:N//2])
双边未求绝对值的幅度谱代码
# 双边未求绝对值的幅度谱
plt.figure(figsize=(7,5))
plt.subplot(211)
plt.plot(frq,mf_100)
plt.title('100-2-3心电信号幅度谱(双边未求绝对值)')
plt.xlabel("频率(HZ)")
plt.ylabel("幅值")
plt.subplot(212)
plt.plot(frq,mf_105)
plt.title('105-2-3心电信号幅度谱(双边未求绝对值)')
plt.xlabel("频率(HZ)")
plt.ylabel("幅值")
plt.show()
双边未求绝对值的幅度谱结果
双边求绝对值的幅度谱代码
# 双边求绝对值的幅度谱
plt.figure(figsize=(7,5))
plt.subplot(211)
plt.plot(frq,mag_100)
plt.title('100-2-3心电信号幅度谱(双边求绝对值)')
plt.xlabel("频率(HZ)")
plt.ylabel("幅值")
plt.subplot(212)
plt.plot(frq,mag_105)
plt.title('105-2-3心电信号幅度谱(双边求绝对值)')
plt.xlabel("频率(HZ)")
plt.ylabel("幅值")
plt.show()
双边求绝对值的幅度谱结果
双边求绝对值归一化的幅度谱代码
# 双边求绝对值归一化的幅度谱
plt.figure(figsize=(7,5))
plt.subplot(211)
plt.plot(frq,gui_mag_100)
plt.title('100-2-3心电信号幅度谱(双边求绝对值归一化)')
plt.xlabel("频率(HZ)")
plt.ylabel("幅值")
plt.subplot(212)
plt.plot(frq,gui_mag_105)
plt.title('105-2-3心电信号幅度谱(双边求绝对值归一化)')
plt.xlabel("频率(HZ)")
plt.ylabel("幅值")
plt.show()
双边求绝对值归一化的幅度谱结果
单边求绝对值归一化的幅度谱代码
# 单边求绝对值归一化的幅度谱
plt.figure(figsize=(7,5))
plt.subplot(211)
plt.plot(half_x,gui_half_mag_100)
plt.title('100-2-3心电信号幅度谱(单边求绝对值归一化)')
plt.xlabel("频率(HZ)")
plt.ylabel("幅值")
plt.subplot(212)
plt.plot(half_x,gui_half_mag_105)
plt.title('105-2-3心电信号幅度谱(单边求绝对值归一化)')
plt.xlabel("频率(HZ)")
plt.ylabel("幅值")
plt.show()
单边求绝对值归一化的幅度谱结果
转换为以db单位的值代码
# 通过 20*np.log10()将其转换为以db单位的值
plt.figure(figsize=(7,5))
plt.subplot(211)
plt.plot(f_100,'g')
plt.title('100-2-3心电信号幅度谱')
plt.xlabel("频率(HZ)")
plt.ylabel("幅值")
plt.subplot(212)
plt.plot(f_105, 'violet')
plt.title('105-2-3心电信号幅度谱')
plt.xlabel("频率(HZ)")
plt.ylabel("幅值")
plt.show()
转换为以db单位的值结果
双边未归一化的相位谱代码
angle_mag_100 = 180*np.angle(mf_100)/np.pi #取复数的弧度,并换算成角度
angle_mag_105 = 180*np.angle(mf_105)/np.pi
# 双边未归一化的相位谱
plt.figure(figsize=(7,5))
plt.subplot(211)
plt.plot(frq,angle_mag_100,'g')
plt.title('100-2-3心电信号相位谱(双边未归一化)')
plt.xlabel("频率(HZ)")
plt.ylabel("相位")
plt.subplot(212)
plt.plot(frq,angle_mag_105,'violet')
plt.title('105-2-3心电信号相位谱(双边未归一化)')
plt.xlabel("频率(HZ)")
plt.ylabel("相位")
plt.show()
双边未归一化的相位谱结果
双边未归一化的相位谱(部分)代码
# 双边未归一化的相位谱(部分)
plt.figure(figsize=(7,5))
plt.subplot(211)
plt.plot(frq[0:100],angle_mag_100[0:100],'g')
plt.title('100-2-3心电信号相位谱(双边未归一化)')
plt.xlabel("频率(HZ)")
plt.ylabel("相位")
plt.subplot(212)
plt.plot(frq[0:100],angle_mag_105[0:100],'violet')
plt.title('105-2-3心电信号相位谱(双边未归一化)')
plt.xlabel("频率(HZ)")
plt.ylabel("相位")
plt.show()
双边未归一化的相位谱(部分)结果
相位谱代码
# 相位谱
plt.figure(figsize=(7,5))
plt.subplot(211)
plt.plot(angle_mag_100)
plt.title('100-2-3心电信号相位谱')
plt.xlabel("频率(HZ)")
plt.ylabel("相位")
plt.subplot(212)
plt.plot(angle_mag_105)
plt.title('105-2-3心电信号相位谱')
plt.xlabel("频率(HZ)")
plt.ylabel("相位")
plt.show()
相位谱结果
功率谱代码
# 功率谱
ps_100 = mag_100**2 / N
ps_105 = mag_105**2 / N
ps_100 = 20*np.log10(ps_100[:N//2])
ps_105 = 20*np.log10(ps_105[:N//2])
plt.figure(figsize=(7,5))
plt.subplot(211)
plt.plot(ps_100,'violet')
plt.title('100-2-3心电信号功率谱')
plt.xlabel("频率(HZ)")
plt.ylabel("功率(dB)")
plt.subplot(212)
plt.plot(ps_105,'g')
plt.title('105-2-3心电信号功率谱')
plt.xlabel("频率(HZ)")
plt.ylabel("功率(dB)")
plt.show()
功率谱结果
心电信号均值方差代码
#coding:UTF-8
import scipy.io as scio
import numpy as np
import matplotlib.pyplot as plt
# %matplotlib inline
#设置绘图大小
plt.rcParams['font.sans-serif']=['SimHei']
plt.rcParams['axes.unicode_minus']=False
dataFile_100 = './100-2-3.mat'
dataFile_105 = './105-2-3.mat'
dataFile_109 = './109-2-3.mat'
dataFile_111 = './111-2-3.mat'
data_100 = scio.loadmat(dataFile_100)
data_105 = scio.loadmat(dataFile_105)
data_109 = scio.loadmat(dataFile_109)
data_111 = scio.loadmat(dataFile_111)
Y_100 = data_100['y']
Y_105 = data_105['y']
Y_109 = data_109['y']
Y_111 = data_111['y']
Fs = 360
N = 2048
frq = np.arange(N) #频率数2048个数
half_x = frq[range(int(N/2))] #取一半区间
print ("心电信号数据:",Y_100)
print ("数据长度:",len(Y_100))
avr109 = np.mean(Y_109)
avr111 = np.mean(Y_111)
var109 = np.var(Y_109)
var111 = np.var(Y_111)
print("信号109的均值:",avr109)
print("信号111的均值:",avr111)
print("信号109的方差:",var109)
print("信号111的方差:",var111)
心电信号均值方差结果
心电信号自相关函数代码
Y_109 = np.array(Y_109).flatten()
RY_109 = np.correlate(Y_109, Y_109, mode='full')
plt.figure(figsize=(7,5))
plt.plot(RY_109)
plt.title('心电信号109的自相关函数')
plt.show()
Y_111 = np.array(Y_111).flatten()
RY_111 = np.correlate(Y_111, Y_111, mode='full')
plt.figure(figsize=(7,5))
plt.plot(RY_111)
plt.title('心电信号111的自相关函数')
plt.show()
心电信号自相关函数结果
心电信号互相关函数代码
Y_109 = np.array(Y_109).flatten()
Y_111 = np.array(Y_111).flatten()
RY_109_111 = np.correlate(Y_109, Y_111, mode='full')
plt.figure(figsize=(7,5))
plt.plot(RY_109_111)
plt.title('心电信号109,111的互相关函数')
plt.show()
心电信号互相关函数结果
心电信号加入高斯白噪声代码
#coding:UTF-8
import scipy.io as scio
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
# %matplotlib inline
#设置绘图大小
plt.rcParams['font.sans-serif']=['SimHei'心电信号基于matlab心率检测含Matlab源码 1993期
心电信号基于matlab GUI自适应滤波+平滑滤波+小波滤波心电信号处理含Matlab源码 1809期
心电信号基于matlab Simulink胎儿心电信号提取含Matlab源码 1550期
心电信号基于matlab NLM时间序列心电信号去噪含Matlab源码 1547期