基于Python的心电信号检测与处理

基于Python的心电信号的特征提取、分析与处理：谱分析，相关分析，滤波器设计等，配套数据、完整代码、分析报告

鲁棒最小二乘支持向量机

10732人浏览 · 2021-11-07 10:13:35

鲁棒最小二乘支持向量机 · 2021-11-07 10:13:35 发布

心电信号的特征提取、分析与处理

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)