技术领域
本发明涉及计算机辅助心电信号诊断，更具体地说，涉及一种基于极值法的快速P波检测方法。
背景技术
心电图(ECG)检测因其对人体无害、非侵入性、能较好地表示心脏的电活动等优势，广泛用于临床的心脏功能的分析。P波作为心电图波组中出现的第一个波，它代表着心房的去极化，因此P波的准确检测可以最早发现是否存在心房病变或者是否为异位心搏。然而，由于P波幅度相较于QRS波群和T波更低，能量集中在0.5-10Hz，信噪比低且易被噪声破坏，从而影响检测结果，主要体现在：(1)检测算法可能无法区分因噪声而产生的波动和真实P波，而错误判断P波的位置，影响心脏起搏位置的判断；(2)噪声叠加在P波上使得检测到的P波振幅不在正常范围内，被认为是假性P波；(3)因为噪声的干扰，而误认为当前P波存在，从而忽略了心房颤动、室性早搏和高血钾症等心脏疾病的可能性，影响了病人的及时就医。对上述P波的特点和存在的问题，相关学者提出了小波变换法、波形拟合法、面积法、QRS-T消除法、自适应阈值法等，这些方法很多都已经有90％以上的准确率，但是方法往往都过于复杂，不利于临床应用。现有的很多基于极值法的P波检测算法都是在一定阈值空间内直接搜索极值点，或者固定搜索窗口内的原信号上搜索最大值作为P波峰值点。虽然满足了方法快速、简单的要求，但前者容易忽略阈值空间之外的潜在P波，后者可能会将搜索窗口的边界值或者因基线漂移而抬高的假性P波，误作为P波峰点。在P波存在性的判断上，现有的方法大都使用P波振幅、时限或者结合P波形态的神经网络方法，前者可能因为P波的低信噪比(SNR)而使得P波的振幅不在标准阈值内，从而错误地把正确的P波认为是假性P波，而后者使用神经网络太耗时。
发明内容
本发明要解决的技术问题在于，针对现有技术的不足，提供一种基于改进极值法的快速P波检测方法，并使用伪TP间隔作为验证P波存在的依据。该方法在噪声干扰下更具有鲁棒性，且快速、简单，能更好地适用于临床。
本发明通过下述技术方案实现：一种基于极值法的快速P波检测方法，包括如下步骤：
A)对采集得到的心电信号s(t)直接进行R峰检测，以分割单个心拍；
B)对所述采集得到的心电信号s(t)进行信号预处理，以提高信号的信噪比和后续检测的准确性；
C)利用所述R峰的位置确定P、T波的搜索范围；
D)在所述的P、T波搜索范围内搜索所有的极大值点，并根据得到的极大值点的个数和极大值点所在处的振幅确定P、T波峰的位置；
E)利用三角形法和拟合最小值法确定P、T波的起始点；
F)计算伪TP间隔，并根据伪TP间隔是否小于某个阈值来验证其是否为真实P波。
更进一步地，所述步骤B)中，进一步包括如下步骤：
B1)使用加窗FIR带通滤波器对采样信号进行初步滤波，并将其应用于前向、后向级联滤波器，得到一个没有相移的纯净信号；
B2)对所述纯净信号连续进行两次中值滤波以消除QRS波群，得到剩余信号；
B3)对所述剩余信号减去由移动滤波器得到的近似基线，以消除基线漂移。
附图说明
图1是本发明一种基于极值法的快速P波检测方法流程图
图2是本发明中各预处理子步骤处理后的信号结果图
图3是本发明中基于改进极值法确定P波的原理图
图4是本发明中基于改进极值法确定T波的原理图
图5是本发明中基于三角形法确定P、T波起终点的原理图
图6是本发明在QT公开数据库中的sel100记录的检测结果图
具体实施方式
下面将结合附图对本发明实施例作进一步说明，但本发明的实施方式不限于此。
如图1所示，本发明基于极值法的快速P波检测方法，包括如下步骤：
步骤A)：对采集得到的心电信号是S(t)，对其直接进行R峰检测，以分割单个心拍：基于心电信号的计算机分析往往需要逐拍分析，因此R峰的准确检测是心电信号分析的第一步。本发明采用了已经相对成熟和广泛应用的Pon&Pomkings法：
所述Pon&Pomkings法先将ECG信号转换为ECG斜率信号，以增强R峰与其他部分采样增幅的差异；
然后将所述ECG斜率信号通过平方运算将所有数据值转换为正数，从而使ECG信号中的高频成分得到增强；
最后使用移动窗口平均滤波器来逐个积分，并将阈值设置为积分信号输出的平均值和最大值的乘积来检测R波峰。
尽管R峰的检测已经很准确，但P、T峰的检测窗口非常依赖R峰的位置。因此，为了使二者的检测更为准确，本发明在检测到R峰后对其进行了修正。考虑到R峰在II导联通常是高尖向上的，我们在初步检测到的R峰附近搜索最大值来修正R峰位置，搜索区域的大小为：
[rloc-0.12*fs:rloc+0.06*fs]
rloc为检测到的初始R峰位置，fs为心电信号记录的采样频率。
步骤B)：对采集得到的心电信号s(t)，进行预处理：心电信号作为一种非线性非平稳的弱电信号，主要幅值范围在10uV～4mv之间，频率范围在0.05-100Hz之间。在心电信号的采集过程中常伴有工频干扰(50Hz)、肌电干扰(30-300Hz)及基线漂移(0.15-0.3Hz)等噪声，与心电信号的频带有所重叠，故在对ECG信号进行检测和分类之前需要对其进行滤波，更进一步的：
步骤B1)：本发明首先使用截止频率选为[1，100]Hz，滤波器阶数为32的加窗FIR带通滤波器(选用’hamming’窗)初步滤除所述采样信号s(t)中的上述信号噪声，并将其应用于前向、后向级联滤波器，得到一个没有相移的纯净信号S1(t)；
步骤B2)：对所述的纯净信号S1(t)连续进行两次中值滤波(窗口L＝2*0.03*fs+1)以消除QRS波群，突出P、T波，得到剩余信号S2(t)；
步骤B3)对所述剩余信号减去由移动滤波器(窗口长度L＝2*0.05*fs)得到的近似基线，以消除基线漂移，得到滤波信号S3(t)。
从图2中我们可以看到，整个预处理流程结束后，QRS基本消除且P、T波的形态更加明显，极值点个数也明显减少，有利于后续的极值点检测。
步骤C)：利用所述R峰的位置确定P、T波的搜索范围：在得到滤波信号S3(t)和已知R峰的前提下，根据P、T波相对R峰的位置和当前心跳的RR间隔设计了P、T波的搜索范围如下：
用于P波检测的搜索窗口：已知正常的PR间期在0.12-0.2s之间，PR段持续时间在0.02-0.12s之间，QRS的持续时间<0.12s。故设置P波搜索窗口为：
[rloc(i)-0.35*RR(i):rloc(i)-0.1*RR(i)]
rloc(i)为当前R峰值位置，RR(i)为当前RR间隔大小，i＝1，2，…，N，N为当前心电记录的总心拍个数。
同理，用于T波检测的搜索窗口：QT间期是指QRS复合波的起点至T波终点之间的持续时间，其正常值在0.35-0.44s之间，QT间期的异常可引起室性心律失常，因此T波搜索窗口的确定与QT间期相关。故设置T波搜索窗口为：
[rloc(i)+0.15*RR(i):rloc(i)+0.57*RR(i)]
步骤D)：在所述的P、T波搜索范围内搜索所有的极大值点，并根据得到的极大值点的个数和极值点所在处的振幅确定P、T波峰的位置；由所述P、T搜索窗口得到对应的信号片段tecg_piece，但是所述信号片段tecg_piece仍然可能存在基线漂移。因此在搜索极值前需要对其再做一次去均值得到ecg_bay。
然后，通过平方运算将信号片段tecg_piece的所有数值转换为正值，从而在经分化后的ECG信号中的高频成分得到进一步的增强，得到ecg_piece。继而计算搜索窗口内的所有极大值，针对不同的极大值个数进行不同的处理：
P波检测：若极值点个数＝0，则我们把当前RR间隔的0.2作为P峰位置(图3(a))；
如果，极值点个数>1，且最后一个极值点位置处对应的ecg_piece处的振幅小于设置的阈值，则以倒数第二个极值点为P峰(图3(b))，定义如下：
ecg_bay[loc]-min(ecg_bay)＜THR
THR＝0.1*[max(ecg_bay)-min(ecg_bay)]
loc为极值点位置，ecg_bay[loc]为极值点位置处的振幅，min(ecg_bay)和max(ecg_bay)分别为ecg_bay的最小值和最大值。
否则，以最后一个极值点为P峰(图3(c))。
T波检测：若极值点个数＝0，则我们把当前RR间隔的0.15作为T峰位置(图4(a))；
如果，极值点个数>2，且第二个极值点位置处对应的ecg_bay处的振幅小于设置的阈值，则以第三个极值点为T峰(图4(b))；
否则，若极大值个数>1并且第一个极值点位置处对应的ecg_bay处的振幅过小，则以第二个极值点为T峰(图4(c))；
除此之外的其它情况，以第一个极值点为T峰(图4(d))。
然而，在得到P、T波峰后，考虑到先前滤波的相移问题，可能会导致原P、T波峰的位置发生偏移，造成检测误差，因此本文针对检测到的P、T波峰进行了局部优化处理。当波峰位置大于0.03倍的当前RR间隔时，在原信号中在波峰前后0.03倍的RR间隔寻找最大值，否则，在波峰到波峰后的0.03倍的RR间隔寻找最大值作为最终的P、T峰位置。
步骤E)利用三角形法和拟合最小值法确定P、T波的起始点，得到伪TP间隔：在检测到P、T波峰后，为了后续的P波存在性检测，就需要在波峰前后的信号片段内搜寻P、T波的起终点。不同于最开始的FIR滤波器截止频率的选择，在进行波形的起终点检测前，我们选择截止频率选为[10，60]Hz的加窗FIR带通滤波对采样信号进行滤波。根据P、T波的形态和相对位置，P波起点和T波的起终点，我们都用三角形法来检测，特别的，对P波终点，我们使用拟合最小值法来检测。具体如下：
P波起点和T波的起终点(三角形法)：如图5所示，以T波起点K点的检测为例，我们先固定T波峰点(固定点po)，然后大致在T波峰前定义J点(临时固定点pt)，从当前R峰值之后的当前RT间隔的四个选项中选择J点位置：{0.3,0.35,0.4,0.45}。连接J、T两点，则K点在JT波段上(移动顶点px)，当三角形JK’T的面积达到最大值时，对应的顶点K’即为K点。三角形的面积我们用海伦公式求得：
pipj表示三角形的边长，s为三角形的半周长：
类比于T波起点的检测，P峰起点pstart为P峰到P峰之前的0.2*RR(i)区间内最大三角形处对应的点；T波终点为T峰到T峰后0.2*RR(i)区间内最大三角形处对应的点。
P波终点(拟合最小值法)：从当前P峰值之后的当前PR间隔的四个选项中选择Pend’点位置：{0.3,0.35,0.4,0.45}，作为终点搜索区间。在所述终点搜索区间内做线性拟合，使用时间间隔的最小二乘近似曲线，求最小值即为终点。
步骤F)计算伪TP间隔，并根据伪TP间隔是否小于某个阈值来验证其是否为真实P波：针对以上的P、T波检测，都是在假定搜索窗口内一定存在P波或T波，但实际心电信号并不是每个心拍都存在P波，若判断错误，则会将异常心拍误认为是正常心拍，不利于后续的诊断，因此，我们会对P波的存在性进行进一步的判段，本发明提出了使用TP间隔来进行判别。
根据假设，每个心跳的P波与前一个心跳的T波的之间一定存在一个TP间隔，我们将TP间隔定义为T波终点与P波起点之间的持续时间。通过T、P波之间的相对位置和P波消失时TP间隔相对较小，结合心率HR，最终将用于判定P波存在的阈值设置为：
THR1＝0.16，若HR<105
THR1＝HR*0.55/fs，其它
THR2＝0.1*RR(i)
并且规定：如果TP间隔小于min{THR1,THR2}，则我们认为不存在P波，由此完成P波检测。
本发明所提供的技术方案为本发明的一个较佳的实现方式，但本发明的实施方式并不局限于此，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化、均应为等效的置换方式，都包含在本发明的保护范围之内。