关键词：
这是前段时间做的一个课程设计，做的比较简单，没有考虑到太细，只是初步地达到了想要的效果。这次设计主要是对心电信号进行预处理，将其信号中包含的一些干扰滤除或者抑制掉。
一、心电信号
（1）心电信号的特性
①工频干扰②基线漂移③肌电干扰
心电信号具有以下几个特点：
·信号极其微弱，一般只有0.05～4mV，典型值为1mV；
·频率范围较低，频率范围为0.1～35Hz，主要集中在5～20Hz；
·存在不稳定性。人体内部各器官问的相互影响以及各人的心脏位置、呼吸、年龄、是否经常锻炼等因素，都会使心电信号发生相应变化；
·干扰噪声很强。对心电信号进行测量时，必然要与外界联系，但由于其自身的信号非常微弱，因此，各种干扰噪声非常容易影响测量。
其噪声可能来自工频(50Hz)干扰、电极接触噪点、运动伪迹、肌电噪声、呼吸引起的基线漂移和心电幅度变化以及其他电子设备的机器噪声等诸多方面。
（2）心电信号的选择
本次实验所采用的心电信号来自MIT-BIH库（心律失常数据库，关于这个库的介绍可以参考：点击打开链接），库中有48组失常的心电信号，要在其中找出符合实验要求的心电信号（即含有肌电干扰、工频干扰和基线漂移）。
（3）正常心电信号波形
图1 正常心电信号波形
图1是正常心电信号在一个周期内的波形，由P波、QRS波群和T波组成。
P波是由心房的去极化产生的，其波形比较小，形状有些圆，幅度约为0.25mV，持续时间为0.08~0.11s。窦房结去极化发生在心房肌细胞去极化之前，因而在时间上要先于P波，只是窦房结处于心脏内部，其电活动在体表难以采集。
P-R间期是指P波起点和QRS波群起点所跨越的时间，是窦房结产生的兴奋，经过右心房、左心房、房室交接区、房室束、左右束支之后，传到到心室所需要的时间。在正常的体表心电图中，P-R间期的值为0.12~0.2s，其中大部分时间是兴奋在房室交界区内传导所需要的时间。P-R间期也称为房室传导时间。
P-R段是指P波终点和QRS波群起点之间所跨越的时间。在正常的体表心电图中，P-R段的心电信号电位值都是接近基线水平的很小点位。在P-R段期间，左右心房同时兴奋，因而两者产生的综合电场对体表心电图的影响较小。另外，此时的兴奋还处于房室交界区和房室束特殊传导系统中，没有到达心室，因而没有产生较大波动的体表心电图信号。
QRS波群是左右心室肌细胞一次发生去极化所产生的膜外负电位在体表的反应。QRS波群的持续时间为0.06~0.1s。由于心室肌细胞在兴奋过程中的综合电场向量多次发生改变，因而形成了体表心电图中大小和方向多次发生变化的心电信号，其中QRS波群中第一个向下的波为Q波，第一个向上的波为R波，R波后面的为S波。
S-T段是指QRS波群终点和T波起点之间所跨越的时间。S-T段期间，左右心室的肌细胞都处于兴奋期间，因而两者形成的综合电场向量在体表心电图中的贡献非常小，导致S-T段心电信号处于大约基线的水平。
T波由心室肌细胞的复极化产生，其幅度为0.1~0.8mV，持续时间为0.05~0.25s。由于复极化差异的存在，T波的方向和QRS波群主波的方向一致。在R波向上的情况下，T波的幅度一般都超过R波幅度的1/10。Q-T间期是指QRS波群起点和T波终点所跨越的时间段，代表心室肌细胞开始去极化到结束复极化所需要的时间，与心率呈负相关。
二、滤波器的选择
1.肌电干扰的滤除—低通滤波器
通常来说，肌电信号的频率为20~5000HZ，其主要成分的频率与肌肉的类型有关，一般在30~300HZ，而心电信号的频率主要集中在5~20HZ，所以选择低通滤波器来滤除肌电干扰。
巴特沃斯滤波器的特点是通频带内的频率响应曲线最为平坦，没有起伏，而在阻频带则逐渐下降为零。巴特沃斯滤波器的振幅对角频率单调下降，并且滤波器的阶数越高，在阻频带幅度衰减速度越快，其他滤波器高阶的振幅对角频率图和低阶数的振幅对角频率有不同的形状。
2.工频干扰的抑制—带陷滤波器
工频干由于供电网络无所不在，因此50Hz的工频干扰是最普遍的，也是心电信号的主要干扰来源。50HZ陷波器的软件设计方法多种多样，常见方法有小波变换滤波、自适应滤波、模板匹配滤波等，但都需要手工计算获得滤波器的参数，运算比较复杂。
滤波器设计中，使用IIR滤波器，可使阶数降低，运算量减少，但破坏了相位特性；使用FIR滤波器既能得到很好的滤波效果，是波形失真达到最下，而且，FIR滤波器可以做成线性相位特性，这正好是心电信号滤波所需要的。
利用MATLAB设计FIR滤波器的方法有窗函数法、频率抽样法和切比雪夫逼近法等，本次课设采用窗函数法设计50HZ陷波滤波器。窗函数方法的基本思想是：首先根据要求选择一个适当的理想低通滤波器，因为其脉冲响应是非因果且无限长的，用最优化窗结构窗函数来截取它的脉冲响应，从而得到线性相位和因果的FIR滤波器。Kaiser窗是接近最优化窗结构的窗函数，它可以根据不同的参数调整滤波器的各项指标，因此采用Kaiser窗函数进行滤波器设计扰的抑制—带陷滤波器
3.基线漂移的纠正—零相移滤波器
零相移滤波器是指一个信号序列经过该滤波器滤波后相位不发生变化，即该滤波器系统函数的相位响应为零。显然，对于因果系统来说是不可能实现零相移的，在事先无法知道信号相位谱的情况下，实现零相移是不可能的。零相移只能是对非因果系统来说的。具体而言，零相移滤波器使用了当前信号点前面和后面的信号点所包含的信息，从本质上说就是使用了“未来的信息”来消除相位失真。
三、程序及结果
1.心电信号读取
因为对MIT-BIH库不是很熟悉，在官网上看过之后，还是不懂（全英文，而且是医学方面的。。。）。所以，此处的心电信号的读取程序是来自网上的 rddata.m ，详细地可以参考点击打开链接。如果自己要用的话，在选取好要处理的心电信号后，把路径更改，并选取合适的样本数，就可以了。我选取的是MIT-BIH中的 109，样本数为1500，下图为心电信号读取后的图形：
图2 MIT-BIH库109心电信号图（双通道）
从图2红色曲线可以看到，波形上存在许多“毛刺”，并且其相位在发生变化（以波峰为例，各波峰大致不在一条水平线上，即所说的“基线漂移”），部分波形收到的干扰比较严重，比较符合对信号处理的要求。
2.心电信号的预处理
（1）肌电信号的滤除
clc;
%------------------------------低通滤波器滤除肌电信号------------------------------
Fs=1500; %采样频率
fp=80;fs=100; %通带截止频率，阻带截止频率
rp=1.4;rs=1.6; %通带、阻带衰减
wp=2*pi*fp;ws=2*pi*fs;
[n,wn]=buttord(wp,ws,rp,rs,'s'); %'s'是确定巴特沃斯模拟滤波器阶次和3dB
截止模拟频率
[z,P,k]=buttap(n); %设计归一化巴特沃斯模拟低通滤波器，z为极点，p为零点和k为增益
[bp,ap]=zp2tf(z,P,k) %转换为Ha(p),bp为分子系数，ap为分母系数
[bs,as]=lp2lp(bp,ap,wp) %Ha(p)转换为低通Ha(s)并去归一化，bs为分子系数，as为分母系数
[hs,ws]=freqs(bs,as); %模拟滤波器的幅频响应
[bz,az]=bilinear(bs,as,Fs); %对模拟滤波器双线性变换
[h1,w1]=freqz(bz,az); %数字滤波器的幅频响应
m=filter(bz,az,M(:,1));
figure
freqz(bz,az);title('巴特沃斯低通滤波器幅频曲线');
figure
subplot(2,1,1);
plot(TIME,M(:,1));
xlabel('t(s)');ylabel('mv');title('原始心电信号波形');grid;
subplot(2,1,2);
plot(TIME,m);
xlabel('t(s)');ylabel('mv');title('低通滤波后的时域图形');grid;
N=512
n=0:N-1;
mf=fft(M(:,1),N); %进行频谱变换（傅里叶变换）
mag=abs(mf);
f=(0:length(mf)-1)*Fs/length(mf); %进行频率变换
figure
subplot(2,1,1)
plot(f,mag);axis([0,1500,1,50]);grid; %画出频谱图
xlabel('频率(HZ)');ylabel('幅值');title('心电信号频谱图');
mfa=fft(m,N); %进行频谱变换（傅里叶变换）
maga=abs(mfa);
fa=(0:length(mfa)-1)*Fs/length(mfa); %进行频率变换
subplot(2,1,2)
plot(fa,maga);axis([0,1500,1,50]);grid; %画出频谱图
xlabel('频率(HZ)');ylabel('幅值');title('低通滤波后心电信号频谱图');
wn=M(:,1);
P=10*log10(abs(fft(wn).^2)/N);
f=(0:length(P)-1)/length(P);
figure
plot(f,P);grid
xlabel('归一化频率');ylabel('功率(dB)');title('心电信号的功率谱'); 
图3 图4
图5 图6
图3是所设计的巴特沃斯数字低通滤波器的幅频响应曲线，图3是在时域滤波前后心电信号的波形图，图5是在频域滤波前后心电信号的频谱图，图6是心电信号的功率谱图
（2）工频干扰的抑制
%-----------------带陷滤波器抑制工频干扰-------------------
%50Hz陷波器：由一个低通滤波器加上一个高通滤波器组成
%而高通滤波器由一个全通滤波器减去一个低通滤波器构成
Me=100; %滤波器阶数
L=100; %窗口长度
beta=100; %衰减系数
Fs=1500;
wc1=49/Fs*pi; %wc1为高通滤波器截止频率，对应51Hz
wc2=51/Fs*pi ;%wc2为低通滤波器截止频率，对应49Hz
h=ideal_lp(0.132*pi,Me)-ideal_lp(wc1,Me)+ideal_lp(wc2,Me); %h为陷波器
冲击响应
w=kaiser(L,beta);
y=h.*rot90(w); %y为50Hz陷波器冲击响应序列
m2=filter(y,1,m);
figure
subplot(2,1,1);plot(abs(h));axis([0 100 0 0.2]);
xlabel('频率(Hz)');ylabel('幅度(mv)');title('陷波器幅度谱');grid;
N=512;
P=10*log10(abs(fft(y).^2)/N);
f=(0:length(P)-1);
subplot(2,1,2);plot(f,P);
xlabel('频率(Hz)');ylabel('功率(dB)');title('陷波器功率谱');grid;
figure
subplot (2,1,1); plot(TIME,m);
xlabel('t(s)');ylabel('幅值');title('原始信号');grid;
subplot(2,1,2);plot(TIME,m2);
xlabel('t(s)');ylabel('幅值');title('带阻滤波后信号');grid;
figure
N=512
subplot(2,1,1);plot(abs(fft(m))*2/N);axis([0 100 0 1]);
xlabel('t(s)');ylabel('幅值');title('原始信号频谱');grid;
subplot(2,1,2);plot(abs(fft(m2))*2/N);axis([0 100 0 1]);
xlabel('t(s)');ylabel('幅值');title('带阻滤波后信号频谱');grid;
其中，ideal_lp（）函数在另一个M文件中，具体如下：
%理想低通滤波器
%截止角频率wc，阶数Me
function hd=ideal_lp(wc,Me)
alpha=(Me-1)/2;
n=[0:Me-1];
p=n-alpha+eps; %eps为很小的数，避免被0除
hd=sin(wc*p)./(pi*p); %用Sin函数产生冲击响应 
图7 图8
图9
图7是带陷滤波器的幅度谱和功率谱，从图中可以看到在50Hz处，滤波器的幅度很大，而且功率在-150以下，说明带陷性能较好。图8是在时域滤波前后的心电信号图，可以看出，滤波后波形有了略微的改善。图19是在频域滤波前后的心电信号频谱图。
（3）基线漂移的纠正
%------------------IIR零相移数字滤波器纠正基线漂移-------------------
Wp=1.4*2/Fs; %通带截止频率
Ws=0.6*2/Fs; %阻带截止频率
devel=0.005; %通带纹波
Rp=20*log10((1+devel)/(1-devel)); %通带纹波系数
Rs=20; %阻带衰减
[N Wn]=ellipord(Wp,Ws,Rp,Rs,'s'); %求椭圆滤波器的阶次
[b a]=ellip(N,Rp,Rs,Wn,'high'); %求椭圆滤波器的系数
[hw,w]=freqz(b,a,512);
result =filter(b,a,m2);
figure
freqz(b,a);
figure
subplot(211); plot(TIME,m2);
xlabel('t(s)');ylabel('幅值');title('原始信号');grid
subplot(212); plot(TIME,result);
xlabel('t(s)');ylabel('幅值');title('线性滤波后信号');grid
figure
N=512
subplot(2,1,1);plot(abs(fft(m2))*2/N);
xlabel('频率(Hz)');ylabel('幅值');title('原始信号频谱');grid;
subplot(2,1,2);plot(abs(fft(result))*2/N);
xlabel('频率(Hz)');ylabel('幅值');title('线性滤波后');grid;s
ubplot(2,1,2);plot(abs(fft(result))*2/N);
xlabel('线性滤波后信号频谱');ylabel('幅值');grid;
图10 图11
图10是在时域滤波前后的心电信号图，可以看出，滤波后基线漂移得到了改善，图11是在频域滤波前后的心电信号频谱图。
心电信号基于matlabgui自适应滤波+平滑滤波+小波滤波心电信号处理含matlab源码1809期(代码片段)
一、心电信号处理简介1引言ECG是一种基本的人体生理信号,具有重要的临床诊断价值。其特点是信号微弱,信噪比小,一般正常人的心电信号频率在0.05～100Hz范围内,幅度为10μV(胎儿)～5mV(成人)。在检测心电信号时,易受到仪... 查看详情
心电信号基于matlab瞬时抑制心电信号iir滤波含matlab源码1533期(代码片段)
...方式1：完整代码已上传我的资源：【心电信号】基于matlab心电信号PTT+HRV+PRV【含Matlab源码1551期】获取代码方式2：通过订阅紫极神光博客付费专栏，凭支付凭证，私信博主，可获得此代码。备注：... 查看详情
心电信号基于matlab心电信号ptt+hrv+prv含matlab源码1551期(代码片段)
...方式1：完整代码已上传我的资源：【心电信号】基于matlab心电信号PTT+HRV+PRV【含Matlab源码1551期】获取代码方式2：通过订阅紫极神光博客付费专栏，凭支付凭证，私信博主，可获得此代码。备注：... 查看详情
心电信号基于matlabsimulink胎儿心电信号提取含matlab源码1550期(代码片段)
一、心电信号简介0引言心电信号是人类最早研究的生物信号之一,相比其他生物信号更易于检测,且具有直观的规律。心电图的准确分析对心脏病的及早治疗有重大的意义。人体是一个复杂精密的系统,有许多不可抗的外界因素,得... 查看详情
心电信号基于matlab心电图峰值检测含matlab源码1548期(代码片段)
...方式1：完整代码已上传我的资源：【心电信号】基于matlab心电图峰值检测【含Matlab源码1548期】获取代码方式2：通过订阅紫极神光博客付费专栏，凭支付凭证，私信博主，可获得此代码。备注：订阅紫... 查看详情
心电信号基于matlab心电图pqrst检测含matlab源码1549期(代码片段)
...方式1：完整代码已上传我的资源：【心电信号】基于matlab心电图PQRST检测【含Matlab源码1549期】获取代码方式2：通过订阅紫极神光博客付费专栏，凭支付凭证，私信博主，可获得此代码。备注：订阅紫... 查看详情
心电信号基于matlabnlm时间序列心电信号去噪含matlab源码1547期(代码片段)
一、心电信号简介0引言心电信号是人类最早研究的生物信号之一,相比其他生物信号更易于检测,且具有直观的规律。心电图的准确分析对心脏病的及早治疗有重大的意义。人体是一个复杂精密的系统,有许多不可抗的外界因素,得... 查看详情
心电信号基于matlab小波阙值心电信号去噪含matlab源码2188期
一、小波阈值法去噪概述电能质量扰动信号的噪声大多以高斯白噪声的形式存在，利用小波变换对信号进行多分辨率分解，由于小波变换具有去除数据相关性的特点，故可以将有用信号与噪声的能量分离开来。信号中... 查看详情
基于python的心电信号检测与处理(代码片段)
心电信号的特征提取、分析与处理1.生物医学信号的特征提取与分析方法2.生物医学信号的滤波方法数据来源：MIT-BIH数据库（可从以下数据中任选两组进行实验）给出4组不同病例的心电信号数据，分别命名为“100-... 查看详情
基于stm32的多功能心电信号监测系统设计
...、运动姿态以及体温生理信息，本文设计实现了一种基于STM32系列MCU的多功能心电信号监测系统，系统包含心电信号采集模块、体表温度采集模块、运动信息测量模块、无线数据传输模块、系统控制模块。实现对人体心电... 查看详情
心电信号基于matlab逻辑算法rst波检测含matlab源码2386期
...)和噪声信号ε(t)经过小波变换多尺度分解后的细节系数。基于小波阈值去噪的流程图如图1所示。图1小波去噪流程小波阈值去噪的具体步骤如下:(1)多尺度分解:根据含噪信号的特点选择适合的小波基和分解层数，经过离散小波... 查看详情
数字信号去噪基于matlab粒子滤波器与fbsmap平滑方法数字信号去噪含matlab源码2179期
一、平滑分解简介根据奈奎斯特定理,采样频率必须大于等于有用信号最高频率的2倍。假设对心电信号的采样满足奈奎斯特采样定理,则实测信号采样频率的1/2为有用心电信号的最高频率。通过三点平滑滤波,可以将频率大于1/2采... 查看详情
基于stm32的无线蓝牙心电监护仪系统设计-毕设课设资料
本次设计的基于STM32的心电信号采集系统由四个模块组成：STM32F103VET6主控模块、OLED显示屏模块，蓝牙模块，AD8232模拟前端模块。在本设计中首先通过三路导联获取人体的心电信号，三路导联线分别与人体的左胸、... 查看详情
毕设题目：matlab数字信号处理
...成案例（代码+参考文献）1.【数字信号处理】基于matlabGUIIIR低通+FIR高通信号时域+频谱分析【含Matlab源码1029期】2.【数字信号处理】基于matlabGUI频谱分析仪【含Matlab源码932期】3.【数字信号处理】基于matlabGUI手机... 查看详情
全套完结数字信号处理----全套matlab实验报告建议保存
...实验一、熟悉MATLAB环境和基本信号的产生与运算实验二、基于MATLAB的离散时间系统的响应实验三、基于MATLAB的离散时间信号的频域分析实验四、离散傅立叶变换的应用实验五、基于MATLAB的模拟滤波器设计实验六、基于MATLAB的IIR数... 查看详情
深度学习heartpy心电图分析(代码片段)
目录1heartpy介绍2使用Pands读取数据2.1数据说明2.2心电图的绘制3心电滤波4心电特性数据读取处理5ECG信号处理相关的开源Python库5.1NeuroKit25.2hrv1heartpy介绍该库提供了处理以下几种信号的方法：来自智能手表和智能手环的常规PPG信... 查看详情
基于matlab的语音信号处理(代码片段)
一、图形界面设计1.新建GUI界面2.新建空白页3.命名为"yydsp"，打开界面 4.拖放控件 5.按预定功能修改界面6.填写Callback函数未填写前的代码：functionvarargout=yydsp(varargin)%YYDSPMATLABcodeforyydsp.fig%YYDSP,byitself,createsanewYYDSPorraisestheexis 查看详情
将 MIT-BIH 心律失常心电图数据库加载到 MATLAB
】将MIT-BIH心律失常心电图数据库加载到MATLAB【英文标题】：LoadMIT-BIHArrhythmiaECGdatabaseontoMATLAB【发布时间】：2011-09-1103:40:46【问题描述】：我正在使用涉及模式识别的神经网络进行心电图信号处理。由于我需要从Matlab收集所有数... 查看详情