基于Python的心电信号检测与处理_python 检测r波位置-程序员宅基地
技术标签： python 笔记 一起学Python 傅立叶分析 iir滤波器 数字信号处理
心电信号的特征提取、分析与处理
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程序较多，大家按需参考学习使用。
本文希望对大家有帮助，当然上文若有不妥之处，欢迎指正。
分享决定高度，学习拉开差距
智能推荐
while循环&CPU占用率高问题深入分析与解决方案_main函数使用while(1)循环cpu占用99-程序员宅基地
文章浏览阅读3.8k次，点赞9次，收藏28次。直接上一个工作中碰到的问题，另外一个系统开启多线程调用我这边的接口，然后我这边会开启多线程批量查询第三方接口并且返回给调用方。使用的是两三年前别人遗留下来的方法，放到线上后发现确实是可以正常取到结果，但是一旦调用，CPU占用就直接100%(部署环境是win server服务器)。因此查看了下相关的老代码并使用JProfiler查看发现是在某个while循环的时候有问题。具体项目代码就不贴了，类似于下面这段代码。while(flag) {//your code;}这里的flag._main函数使用while(1)循环cpu占用99
【无标题】jetbrains idea shift f6不生效_idea shift +f6快捷键不生效-程序员宅基地
文章浏览阅读347次。idea shift f6 快捷键无效_idea shift +f6快捷键不生效
node.js学习笔记之Node中的核心模块_node模块中有很多核心模块,以下不属于核心模块,使用时需下载的是-程序员宅基地
文章浏览阅读135次。Ecmacript 中没有DOM 和 BOM核心模块Node为JavaScript提供了很多服务器级别，这些API绝大多数都被包装到了一个具名和核心模块中了，例如文件操作的 fs 核心模块 ，http服务构建的http 模块 path 路径操作模块 os 操作系统信息模块// 用来获取机器信息的var os = require('os')// 用来操作路径的var path = require('path')// 获取当前机器的 CPU 信息console.log(os.cpus._node模块中有很多核心模块,以下不属于核心模块,使用时需下载的是
数学建模【SPSS 下载-安装、方差分析与回归分析的SPSS实现(软件概述、方差分析、回归分析)】_化工数学模型数据回归软件-程序员宅基地
文章浏览阅读10w+次，点赞435次，收藏3.4k次。SPSS 22 下载安装过程7.6 方差分析与回归分析的SPSS实现7.6.1 SPSS软件概述1 SPSS版本与安装2 SPSS界面3 SPSS特点4 SPSS数据7.6.2 SPSS与方差分析1 单因素方差分析2 双因素方差分析7.6.3 SPSS与回归分析SPSS回归分析过程牙膏价格问题的回归分析_化工数学模型数据回归软件
利用hutool实现邮件发送功能_hutool发送邮件-程序员宅基地
文章浏览阅读7.5k次。如何利用hutool工具包实现邮件发送功能呢？1、首先引入hutool依赖<dependency> <groupId>cn.hutool</groupId> <artifactId>hutool-all</artifactId> <version>5.7.19</version></dependency>2、编写邮件发送工具类package com.pc.c..._hutool发送邮件
docker安装elasticsearch，elasticsearch-head，kibana，ik分词器_docker安装kibana连接elasticsearch并且elasticsearch有密码-程序员宅基地
文章浏览阅读867次，点赞2次，收藏2次。docker安装elasticsearch，elasticsearch-head，kibana，ik分词器安装方式基本有两种，一种是pull的方式，一种是Dockerfile的方式，由于pull的方式pull下来后还需配置许多东西且不便于复用，个人比较喜欢使用Dockerfile的方式所有docker支持的镜像基本都在https://hub.docker.com/docker的官网上能找到合..._docker安装kibana连接elasticsearch并且elasticsearch有密码
随便推点
Python 攻克移动开发失败！_beeware-程序员宅基地
文章浏览阅读1.3w次，点赞57次，收藏92次。整理 | 郑丽媛出品 | CSDN（ID：CSDNnews）近年来，随着机器学习的兴起，有一门编程语言逐渐变得火热——Python。得益于其针对机器学习提供了大量开源框架和第三方模块，内置..._beeware
Swift4.0_Timer 的基本使用_swift timer 暂停-程序员宅基地
文章浏览阅读7.9k次。//// ViewController.swift// Day_10_Timer//// Created by dongqiangfei on 2018/10/15.// Copyright 2018年 飞飞. All rights reserved.//import UIKitclass ViewController: UIViewController { ..._swift timer 暂停
元素三大等待-程序员宅基地
文章浏览阅读986次，点赞2次，收藏2次。1.硬性等待让当前线程暂停执行，应用场景：代码执行速度太快了，但是UI元素没有立马加载出来，造成两者不同步，这时候就可以让代码等待一下，再去执行找元素的动作线程休眠，强制等待 Thread.sleep(long mills)package com.example.demo;import org.junit.jupiter.api.Test;import org.openqa.selenium.By;import org.openqa.selenium.firefox.Firefox.._元素三大等待
Java软件工程师职位分析_java岗位分析-程序员宅基地
文章浏览阅读3k次，点赞4次，收藏14次。Java软件工程师职位分析_java岗位分析
Java：Unreachable code的解决方法_java unreachable code-程序员宅基地
文章浏览阅读2k次。Java：Unreachable code的解决方法_java unreachable code
标签data-*自定义属性值和根据data属性值查找对应标签_如何根据data-*属性获取对应的标签对象-程序员宅基地
文章浏览阅读1w次。1、html中设置标签data-*的值 标题 11111 222222、点击获取当前标签的data-url的值$('dd').on('click', function() { var urlVal = $(this).data('ur_如何根据data-*属性获取对应的标签对象