基于Python的心电信号检测与处理
基于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']
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))
def awgn(x, snr, seed=7):
'''
加入高斯白噪声 Additive White Gaussian Noise
:param x: 原始信号
:param snr: 信噪比
:return: 加入噪声后的信号
'''
np.random.seed(seed) # 设置随机种子
snr = 10 ** (snr / 10.0)
xpower = np.sum(x ** 2) / len(x)
npower = xpower / snr
noise = np.random.randn(len(x)) * np.sqrt(npower)
return x + noise
SNR=10*np.log(100/8)
SNRY_109 = awgn(Y_109,SNR)
plt.figure(figsize=(7,5))
plt.subplot(211)
plt.plot(Y_109)
plt.title('原始109-2-3心电信号')
plt.xlabel("频率(HZ)")
plt.ylabel("幅值")
plt.subplot(212)
plt.plot(SNRY_109)
plt.title('加入高斯白噪声')
plt.xlabel("频率(HZ)")
plt.ylabel("幅值")
plt.show()
心电信号加入高斯白噪声结果
心电信号加入高频噪声1000HZ代码
dt = 1/1023.5
t = np.array(2)
Y_109_high=np.sin(2*np.pi*1000*t)/10
Y_109_out = Y_109 + Y_109_high
plt.plot(Y_109_out)
plt.title('加入高频噪声1000HZ')
plt.xlabel("频率(HZ)")
plt.ylabel("幅值")
plt.show()
心电信号加入高频噪声1000HZ结果
IIR零相移数字滤波器纠正基线漂移代码
# IIR零相移数字滤波器纠正基线漂移
Wp = 1.4*2/Fs # 通带截止频率
Ws = 0.6*2/Fs # 阻带截止频率
devel = 0.005 # 通带纹波
Rp = 20*np.log10((1+devel)/(1-devel)) # 通带纹波系数
Rs=20 # 阻带衰减
N, Wn = signal.ellipord(Wp,Ws,Rp,Rs,True) # 求椭圆滤波器的阶次
b, a = signal.ellip(N,Rp,Rs,Wn,'high', True) # 求椭圆滤波器的系数
w, h = signal.freqs(b, a, 512)
result = signal.lfilter(b,a,Y_109)
plt.figure(figsize=(7,5))
plt.subplot(211)
plt.plot(Y_109)
plt.title('原始109-2-3心电信号')
plt.xlabel("频率(HZ)")
plt.ylabel("幅值")
plt.subplot(212)
plt.plot(result)
plt.title('线性滤波后信号')
plt.xlabel("频率(HZ)")
plt.ylabel("幅值")
plt.show()
N = 512
plt.figure(figsize=(7,5))
plt.subplot(211)
plt.plot(np.abs(np.fft.fft2(Y_109))*2/N)
plt.title('原始109-2-3心电信号频谱')
plt.xlabel("频率(HZ)")
plt.ylabel("幅值")
plt.subplot(212)
plt.plot(np.abs(np.fft.fft2(result))*2/N)
plt.title('线性滤波去掉基线漂移频谱')
plt.xlabel("频率(HZ)")
plt.ylabel("幅值")
plt.show()
IIR零相移数字滤波器纠正基线漂移结果
其他数字滤波器尝试代码
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
# 巴特沃斯滤波器
b, a = signal.butter(4, 100, 'low', analog=True) # 4阶低通临界频率为100Hz
w, h = signal.freqs(b, a) # h为频率响应,w为频率
result = signal.lfilter(b,a,Y_109)
# 巴特沃斯带通滤波器
N, Wn = signal.buttord([20, 50], [14, 60], 3, 40, True)
b, a = signal.butter(N, Wn, 'band', True)
w, h = signal.freqs(b, a, np.logspace(1, 2, 500))
result = signal.lfilter(b,a,Y_109)
# Chebyshev I型滤波器
b, a = signal.cheby1(4, 5, 100, 'low', analog=True)
w, h = signal.freqs(b, a)
result = signal.lfilter(b,a,Y_109)
#
N, Wn = signal.cheb1ord(0.2, 0.3, 3, 40)
b, a = signal.cheby1(N, 3, Wn, 'low')
w, h = signal.freqz(b, a)
result = signal.lfilter(b,a,Y_109)
# Chebyshev II型滤波器
b, a = signal.cheby2(4, 40, 100, 'low', analog=True)
w, h = signal.freqs(b, a)
result = signal.lfilter(b,a,Y_109)
#
N, Wn = signal.cheb2ord([0.1, 0.6], [0.2, 0.5], 3, 60)
b, a = signal.cheby2(N, 60, Wn, 'stop')
w, h = signal.freqz(b, a)
result = signal.lfilter(b,a,Y_109)
结果整理和报告分析
4组不同病例的心电信号数据：link
心电信号的特征提取、分析与处理的报告分析及完整代码：link
本人能力有限，解释尚不清楚明了，如遇任何问题，大家可留言或私信。开源Matlab程序较多，大家按需参考学习使用。
本文希望对大家有帮助，当然上文若有不妥之处，欢迎指正。
分享决定高度，学习拉开差距
华为开发者空间，是为全球开发者打造的专属开发空间，汇聚了华为优质开发资源及工具，致力于让每一位开发者拥有一台云主机，基于华为根生态开发、创新。
更多推荐
23
0- 0
已为社区贡献3条内容
所有评论(0)