Matlab心电信号QRS波检测

1、QRS波的特征

如图所示为完整的心电信号波形图，分别由P波、P—R段、P—R间期、QRS复合波、S—T段、T波和U波组成。下面将重点讲诉QRS波。并且对它进行相关分析。
QRS复波。代表两个心室兴奋传播过程的电位变化。由窦房结发生的兴奋波经传导系统首先到达室间隔的左侧面，以后按一定路线和方向，并由内层向外层依次传播。随着心室各部位先后去极化形成多个瞬间综合心电向量，在额面的导联轴上的投影，便是心电图肢体导联的QRS复合波。典型的QRS复合波包括三个相连的波动。第一个向下的波为Q波，继Q波后一个狭高向上的波为R波，与R波相连接的又一个向下的波为S波。由于这三个波紧密相连且总时间不超过0.10秒，故合称QRS复合波。QRS复合波所占时间代表心室肌兴奋传播所需时间，正常人在0.06～0.10秒之间。
2、 QRS波检测算法
QRS波群的准确检测一直是心电信号自动分析的焦点和难点。目前QRS波检测方法主要有:差分法(Derivative)、带通滤波法(Bandpass filter) 、小波变换(Wavelet Transfom)、形态学运算(Mathematical Morphology、长度和能量变换(Length and Energy Transtfoms) 等;另外还有一些新兴的研究方法，如人工神经网络(Arificial Neural Networks,ANN)、遗传算法( Genetic Algorithm，GA)、句式分析( Syntactic Methods)、隐Markov 模型( Hidden Markov Models, HMM)、匹配滤波(Matched Filter) 、Hilbert变换(Hilbert Transfom)、心电模板( Tamplate)、过零检测(Zero - Crossing detection)等,多种技术交叉融合(A lgorihms based on the fusion of several technologies)的趋势也日趋明显。
QRS波群检测的经典算法主要有以下几种:
（1）差分法
差分法是通过对信号进行一阶或二阶差分,判断其差分值是否超过特定阈值来确定QRS波的存在及其位置。差分法是一种快速算法，适合于对实时性要求较高的心电自动监护仪。
差分法的基本原理是: QRS波是心电信号波形变化最剧烈的地方,其波形的上升斜率和下降斜率与其它波形的斜率相比显著不同,通过检测心电信序列对时间的导数(即斜率的变化)来定位QRS波的位置。通常在R波的上升沿和下降沿是心电波形斜率变化最大的区域,中间出现的一阶导数过零点为R点所在的位置。
差分算法主要以一阶、二阶差分为基础，比较常用的二阶差分的公式有:
y(n)=2x(n+ 2)+x(n+1)-x(n-1)-2x(n-2)
y(n)= x(n+ 2)- 2x(n)+ x(n-2)
y(n)= 1/8[2x(n)+x(n-1)-x(n-3)-2x(n -4)]
在实际应用中，差分的方法对高频噪声非常敏感，因此，在差分处理之前，通常让ECG信号通过滤波器，以消除高频噪声于扰。一个QRS波的典型频率分量是从10Hz到25Hz,通常设计带通滤波器来消除低频和高频噪声的干扰。
经过差分运算后，R波的特征将转化为信号中的高频分量,为了强化信号中的高频成份，方便后期判别，可以对差分后的数据进行后期处理。常用的处理方法有逐点平方法、窗口求积法等非线性处理方法。
差分的方法具有直观形象、计算量小的优点，但是，也有抗干扰能力差的缺点。因此，适合于那些噪声干扰较小的情况。
（2）带通滤波法
带通滤波法实际上增强了QRS波,同时抑制了高频和低频噪声, FIR带通还能保持线性相位,但性能优良、低阶整系数FIR带通滤波器较难设计,而且频段固定，灵活不足。
带通滤波器的通带接近QRS波的主频带时，对QRS波有增强作用。FIR带通滤波器能够保持信号的线性相位,便于特征点的寻找，FIR非递归的特性也使得运算迅速。整系数的带通滤波器公式如下 ,它是FIR滤波器,能够保持线性相位。
K,L>0
如果对称且因果的FIR滤波器阶数是M，那么处理后信号延时M/2。虽然Matlab软件能够设计FIR滤波器,但是,对称、低阶、整系数并且性能优良的FIR带通滤波器的设计仍然比较困难。
（3）匹配滤波法
视QRS波段为待检信号,视非QRS段信号为噪声,匹配滤波器对信号的各频率分量起到幅度同相相加的作用,对噪声的各分量起到的功率相加的作用,综合而言,信噪比得到提高。设待检信号为s(n),则匹配滤波器的脉冲响应函数是:h(n) =c* s(k- n)。当c=1时,h(n)和s(n)关于k/2对称。即把某个QRS波的数字序列做逆序排列就得到最佳线性匹配滤波器的时域表达式h(n)。
匹配滤波器有以下特点:最大信噪比只与信号的能量和噪声强度成正比,而与信号的波形无关;对信号幅度和时延具有适应性,但对频移不具有适应性。
一般QRS波段的能量远远大于非QRS段能量,QRS波的主频段相对稳定,故可以采用匹配滤波来检测QRS波。滤波器对信号的各频率分量起到幅度同相则相加的作用,对噪声的各分量起到的实功率相加的作用,综合而言,信噪比得到提高。当QRS波的形态和宽度-致性好时,此方法是最佳线性滤波器,但QRS波形态变化显著时,检测性能下降。
（4）长度和能量变换法
对差分后的信号进行增强，较差分法稳定性提高,也适于多导心电,但对高频噪声的抵抗能力仍不尽人意。长度和能量变换的最大优点,就是不但可以用于单导联心电信号的R波检测,而且可以用于多导联R波检测。多导联心电信号的长和能量变换的定义如下:

其中n是导联数, i为时间序列号, q是分析窗长度,差分信号。
如果令n=1,则长度和能量变换就可以直接应用于单导联的R波检测。长度和能量变换与差分法相似，可以消除基线漂移的影响,但对高频干扰比较敏感。长度和能量变换检测R波的准确率要比.差分法高得多,特别在高频干扰较弱的情况下，对幅度较小的R波检测比其它方法更有效。
（5）基于小波变换的分析方法
小波变换是近年来得到迅猛发展的一种数学方法,其作为一种信号的时间尺度分析方法，具有多分辨率分析的特点，从而可以在时、频两域获得表征信号局部特征的能力。因为其可变的时间窗和频率窗,使得它对于信号具有很高的适应性。心电信号的R波、P波和T波以及各种噪声干扰都具有不同的频率内容，小波变换对心电信号作多分辨率分解后，一方面让不同的频率成分在不同的频带内分别显现,实现各个特征点以及噪声干扰的频带划分，另一方面突出它们奇异点(即突变点)的特征。
一个偶对称小波会把信号峰值变成极大值，一个奇对称小波则会把信号峰值变成一个过零点,因此小波变换能够用于信号特征点的检测。用FR滤波器组来实现的小波变换兼顾处理时间,抗千扰能力强,准确率高。但信号的干扰较大时，还有误检现象。
（6）基于数学形态学的方法
通过数学形态学运算,局部修改信号的几何特征后进行检测。基于数学形态学的算法用于心电图波形的分离与常规方法不同。心电信号是一维离散信号,用形态变换对心电信号滤波属于多值变换,基本运算包括腐蚀、膨胀、形态开和形态闭。算法无需检测QRS波，而是用数学方法对心电信号进行滤波、基线矫正等处理之后，直接去除心电图信号的QRS波群，实现波形的定性分离;检出P波和T波的起止点后便可以完成定量分离。整个过程基本通过一系列形态学算子完成的。实际方法是由两路开、闭运算组成。一路先进行开运算，再进行一次比运算;另一路则相反。然后两路运算的结果取平均值。经过这样的运算后，信号中实际被滤除成分从形态上来说，与运算使用的结构元素有关。若采用的结构元素长度为M，那么，输出的信号中将不再含有跨度小于M的波形。所以适当的调节元素结构的长度,便可以有选择的达到噪声滤除、去除QRS波群.以及基线矫正的目的。经过噪声滤除、去除QRS波群以及基线矫正后，可以得到只含有P波和T波、基线基本上为直线的图形，从定性的意义上来说实现了波形的分离。
该算法可有效消除脉冲噪声,不受基漂影响，但是修改过的信号呈“台阶”状,结构子越长则耗时越多(结构元的宽度随采样率的提高也相应增长)，时间消耗大是本方法的显著缺点。
（7）人工神经网络的方法
基于神经网络的QRS波群的检测方法是医学信号处理的-一个新兴领域。神经网络由许多并行运算的功能简单的单元组成，这些单元类似于生物神经的单元。虽然单个神经元的结构及其简单，功能有限，但大量神经元构成的网络系统所能实现的功能却是极其丰富多彩的。神经网络首先使用一些先验样本知识开始训练，这些知识以连接权重和偏置值的形式分布于神经网络的结构中。随着训练样本的增加，神经网络自动调节其结构参数,提取隐含在心电图QRS波群中的特征值，其分类能力也不断得到增强。训练结束后，既可以对训练过的心电图进行正确识别，又能对未训练过的心电图又较好的识别能力。神经网络系统具有集体运算的能力和自适应学习能力，又很强的容错性等，判别速度快，但是网络的训练时间较长。但是神经网络的方法所使用的模板受到ECG形态变异影响较大。
（8）希尔伯特变换的方法
一个时间函数x(t),它的希尔伯特变换定义为:

希尔伯特变换是信号处理中比较常见的-种处理手段,它能够对信号起到一个宽频移相的作用。希尔伯特变换器是一一个全通滤波器，只引起信号频谱的相位变化,频谱的幅度不发生变化。通过计算ECG信号的希尔伯特变换可以得到ECG信号的包络，而通过这个包络可以得到信号的带限。将得到的包络进行低通滤波后可以突出包络中的R波峰值,进而得到R波位置。使用希尔伯特变换法检测R波的流程见下图.

希尔伯特变换具有计算简单。计算量小的优点,也具有抗干扰能力差的缺点。
（9）遗传算法
遗传算法是模仿生命的进化过程和遗传变异理论的基础上的随机搜索优化技术。一般的QRS波检测器都包含两个级联的模块，一个模块用来加强QRS复波,另一个模块用来检测。在这个算法中，我们使用遗传算法来同时优化这两个模块的参数。
检测器的第一个模块是一个多项式滤波器，它用来加强QRS复波。为了实现最大效率，滤波器在非常小数目的输入采样点上进行操作，这些采样点都是由遗传算法来进行选择。第二个模块是一个简单的、自适应的最大值检测器，它可以去除伪峰值。
滤波器的参数和检测器的参数根据已经包含正确检测的ECG信号来优化。在多项式滤波器中，在时刻i的输出值从一个有N个输入的M阶多项式求得:

我们采用了一个最简单的算法，能够探测到滤波器输出的最大值。为了避免在噪声存在的情况下的错检，仅仅那些幅值比阈值Y大的而且没有落在静息期内的最大值被检测到。为了适合ECG信号的变换，采用自适应阈值。为了设计这样一个检测器，需要定义多项式滤波器的特征值、选择它的系数和最大值检测器的参数。这些变量中的一部分由人来设计，另外一部分由遗传算法来选择。
遗传算法是模仿生命的进化过程和遗传变异理论的基础上的随机搜索优化技术。遗传算法检测结果的好坏，受到ECG模板的影响较大。
3、算法设计
Q波和S波通常是低幅高频波，一般Q波位于R波之前，S波位于R波之后，由于他们是一般向下的波，所以他们的峰值点和极值是对应的。所以首先找到信号波形的波峰和波谷，波峰即是最大值，也就是R波处。然后在检测到R波向左和向右分别搜寻到极值点，对应的就是Q波和S波。
在本次设计中，我们使用基于极值的动态自适应阈值法。在利用阈值进行QRS波检测时，若固定阈值，则会造成阈值设置过高导致漏检，产生假阴性，阈值设置过低会导致多测，产生假阳性。因此本文提出基于待测信号的可变阈值，以提高检测的精确率，所采用的可变阈值包括幅度阈值和时间间隔阈值等。
基本原理:基于R波的幅值为最大的特点，根据极值的定义，筛选出所以的极大值点，这些点即可能是R波的点，然后根据心电学原理知识，确定R波的阈值，计算待测的ECG数据最大值并与阈值进行比较，若超过或达到阈值，则初步.判断已检测到一个R波，然后根据制定的规则确定R波。
Step 1：即对待测ECG的滤波处理，由于第四次实验中已经完成了对ECG信号的预处理，因此本次实验中直接调用上次实验结果。
Step2:根据QRS波波形，首先利用极值点判断，一次筛选得到所有极值点，记为sigmax:

结合心电信号原理，确定阈值：

Step3:二次筛选:对于sigmax中数据，若出现大于阈值的数据，则在该数据处和其后50ms的范围内查找，将找到的最大值作为一个R波;接着在此R波后150ms处继续进行二次筛选，直至查找完全部数据，将二次筛选得到的R波记为rvalue。
Step4:可能存在的误差:判断rvalue中前后两R波间隔是否小于0.4，若小于0.4则消去较小值，留下较大值，直至查找全部数据，得到新的R波，刷新rvalue。
Step5:寻找Q波和S波：分别在R波前面和后面寻找第一个极值点，定位到 Q波和S波。
4、部分代码

%取阈值,阈值为相对幅值的差的60%
thrtemp=sort(sigmax);
thrlen=length(sigmax);
thr=0;
for i=(thrlen-7):thrlenthr=thr+thrtemp(i);
end;
thrmax=thr/8;               %最大幅度平均值，8个最大幅值点的平均值zerotemp=sort(z);
zerovalue=0;
for i=1:100zerovalue=zerovalue+zerotemp(i);
end;
zerovalue=zerovalue/100;    %最小幅度平均值，对消幅度，100个最小幅值点的平均值
thr=(thrmax-zerovalue)*0.3; %最大、最小幅度的差值的30%为判别R波的阈值                      %定位R波
rvalue=[];
for i=1:thrlenif sigmax(i,1)>thrrvalue=[rvalue;sigmax(i,2)];end;
end;
rvalue_1=rvalue;%排除误检，如果相邻两个极大值间距小于0.4，则去掉幅度较小的一个
lenvalue=length(rvalue);
i=2;
while i<=lenvalueif (rvalue(i)-rvalue(i-1))*rate<0.4if yabs(rvalue(i))>yabs(rvalue(i-1))rvalue(i-1)=[];elservalue(i)=[];end;lenvalue=length(rvalue);i=i-1;end;i=i+1;
end;        lenvalue=length(rvalue);
% 在原信号上精确校准
for i=1:lenvalueif (sig(rvalue(i))>0)k=(rvalue(i)-5):(rvalue(i)+5);[a,b]=max(sig(k));rvalue(i)=rvalue(i)-6+b; elsek=(rvalue(i)-5):(rvalue(i)+5);[a,b]=min(sig(k));rvalue(i)=rvalue(i)-6+b; end;
end;%打印纠正及校准前后的R波信号
figure(2);
subplot(2,1,1),plot(t,sig,rvalue_1/360,sig(rvalue_1),'r.');
xlabel('时间（s）');ylabel('电压（mv）');title('定位R波信号')
subplot(2,1,2),plot(t,sig,rvalue/360,sig(rvalue),'r.');
xlabel('时间（s）');ylabel('电压（mv）');title('纠正及校准前后的R波信号') %检测Q波
wtsig2=cwt(sig,8,'mexh');
lenrvalue=length(rvalue);
qvalue=[];
for i=1:lenrvaluefor j=rvalue(i):-1:(rvalue(i)-30)if sig(rvalue(i))>0if wtsig2(j)<wtsig2(j-1)&wtsig2(j)<wtsig2(j+1)tempqvalue=j-10;                 %确定检测窗的起点break;                        %正向波，取第一个负极大值end;elseif wtsig2(j)>wtsig2(j-1)&wtsig2(j)>wtsig2(j+1)tempqvalue=j-10;                 %确定检测窗的起点break;                      %倒置R波，取第一个正极大值end;end;end;x1=tempqvalue;y1=sig(tempqvalue);x2=rvalue(i);y2=sig(rvalue(i));a0=(y2-y1)/(x2-x1);b0=-1;c0=-a0*x1+y1;                        %求直线公式参数ax+by+c=0dist=[];for k=tempqvalue:rvalue(i)tempdist=(abs(a0*k+b0*sig(k)+c0))/sqrt(a0^2+b0^2);dist=[dist;tempdist];end;                                  %求点到直线距离[a,b]=max(dist);                      %找到距离最大值，Q波就在附近tempqvalue=tempqvalue+b-1;
%     l=(tempqvalue-5):rvalue(i);
%     [c,d]=min(sig(l));
%     tempqvalue=tempqvalue-6+d;            %在最大值附近修正Q波，得到结果qvalue=[qvalue;tempqvalue];
end;%检测S波
svalue=[];
for i=1:lenrvalue-1for j=rvalue(i):1:(rvalue(i)+100)if sig(rvalue(i))>0if (wtsig2(j)<wtsig2(j-1))&(wtsig2(j)<wtsig2(j+1))tempsvalue=j+10;      %从R波开始向后寻找第一个极小值break;end;elseif (wtsig2(j)>wtsig2(j-1))&(wtsig2(j)>wtsig2(j+1))tempsvalue=j+10;     %从R波开始向后寻找第一个极大值break;end;
end;end;x1=tempsvalue;y1=sig(tempsvalue);x2=rvalue(i);y2=sig(rvalue(i));a0=(y2-y1)/(x2-x1);b0=-1;c0=-a0*x1+y1;                        %求直线公式参数ax+by+c=0dist=[];for k=rvalue(i):tempsvaluetempdist=(abs(a0*k+b0*sig(k)+c0))/sqrt(a0^2+b0^2);dist=[dist;tempdist];end;                                  %求点到直线距离[a,b]=max(dist);                   %找到距离最大值，S波就在附近tempsvalue=rvalue(i)+b-1;l=rvalue(i):(tempsvalue+10);[c,d]=min(sig(l));tempsvalue=rvalue(i)+d-1;     %在最大值附近修正S波，得到结果    svalue=[svalue;tempsvalue];
end;

Matlab心电信号QRS波检测相关推荐

检测心电信号的p波的matlab代码,matlab心电信号R波检测程序.doc
<生物医学信号处理>实习报告学生姓名: 学号:实验室名称:项目名称:心电信号的R波检测项目内容: 总结常用的QRS波检测算法: 选择一种QRS波检测算法,理解该检测算法: 编写程序,检测 ...
【代码补全】matlab心电信号R波提取
对Matlab心电信号QRS波检测中代码的补充,使其完整可用. 补充内容来自基于matlab的心电信号QRS波检测与分析. x=iswt(swa,swd,'db3'); %--心电信号去噪完成后-- ...
QRS波检测算法集锦(含源代码)
下面介绍一些我找到的一些QRS波检测算法的公开源代码 1:Pantompkin[1985年] 2:Filter-Bank[1999年] 3:Phase Space[1999年] 4:RST State ...
Matlab心电信号的PQRST模拟-实验报告
心电信号处理算法设计-实验要求 data4 是一段实际采样得到的心电数据, 采样频率为 100Hz, 波形如下图所示.设计相应的算法, 计算心率, 单位为: 次/分钟.可能会用到的知识为数字滤波器的设 ...
OSEA中QRS波检测算法
当信号经过滤波等预处理后,检测器开始检测任何一个峰值,这峰包括信号中的所有峰值.每次检测到一个峰值,它被分类为QRS波,或者噪声,或者为了后来的分类保存它.算法利用峰值的高度,相对于上一个QRS波的位 ...
OSEA中QRS波检测算法代码分析-未完待续
最近一直在搞R波检测算法,对OSEA代码主要是对注释做一个翻译,增加注释,使代码更容易理解. 一.首先看QRSDE.H /*************************************** ...
基音周期检测 matlab,语音信号基音周期检测的matlab程序
<语音信号基音周期检测的matlab程序>由会员分享,可在线阅读,更多相关<语音信号基音周期检测的matlab程序(2页珍藏版)>请在人人文库网上搜索. 1.function ...
基于压缩感知的心电信号QRS检测算法matlab仿真
目录 1.算法仿真效果 2.MATLAB源码 3.算法概述 4.部分参考文献 1.算法仿真效果 matlab2022a仿真结果如下: 2.MATLAB源码 %********************* ...
matlab求心率,心电图QRS波检测（计算心跳次数）
YURA 2021-2-18 9:54:53 谢谢!!! yangaichimantou 2021-1-6 21:01:21 谢谢分享楼主 llyy1233 2020-12-21 18:39:55 感 ...
Ecg信号QRS波峰检测：A Real-Time QRS Detection Algorithm (Pan-Tompkins法)
Pan-Tompkins法是R峰检测的经典方法,被引量高达6000+次. Pan J, Tompkins W J. A real-time QRS detection algorithm[J]. IE ...

Matlab心电信号QRS波检测

Matlab心电信号QRS波检测

Matlab心电信号QRS波检测相关推荐

最新文章

热门文章