matlab心电信号R波检测Word格式.docx
《matlab心电信号R波检测Word格式.docx》由会员分享，可在线阅读，更多相关《matlab心电信号R波检测Word格式.docx（18页珍藏版）》请在冰点文库上搜索。
若超过0．２05,一般表明有房室传导阻滞的发生＂P一R间期的长短与年龄及心率有关。
（3）QRS波群:
反映两心室去极化过程的电位变化＂典型的ＱRＳ波群包括三个紧密相连的电位波动:
第一个向下的波称为Q波;
紧接着是向上!
高而尖峭的Ｒ波;
最后是向下的S波"
在不同导联中，这三个波不一定都出现,各波的幅度变化也较大"
历时约0.06一0.105＂。
（4）S一Ｔ段:
指QRS波群终点与T波起点之间的线段,一般与零电位基线平齐"
在这段时期内,因心室各部分都已全部进入除极化状态,但尚未开始复极,故心室各部分之间没有电位差存在,心电曲线恢复到基线水平＂但若有冠状动脉供血不足或心肌梗死等情况发生时,S一T段常会偏离基线,并超过一定的幅度范围"
（5）Ｔ波：
反映两心室复极化过程的电位变化"
波形圆钝,升降支并不完全对称，波形的前支较长而后支较短，占时约0.05一0.255"
T波方向应与QRＳ波群的主波方向一致＂在以R波为主的导联中,其波幅应不低于本导联R波的１/10。
（６）Ｑ一T间期：
指从QRS波群起点到T波终点之间的时间，它代表心室开始去极化到全部复极化完毕所需的时间"
这一间期的长短与心率密切相关＂心率越快,Q一T间期越短：
反之,则Q一T间期越长＂正常的Q一T间期依心率!
年龄及性别不同而有所不同.当心率为7５次/分时,Q一T间期为０.3０一0.405"
分析Q一T间期的变化,对疾病的早期诊断和分析抗心律失常药物对心脏的影响,可起到一定的辅助作用"
由于Q一T间期受心率的影响比较大,临床上经常采用修正的Q一T间期，即采用Baｚett公式计算:
（7）U波：
T波后０.0２一0.０4:
可能会出现一个与Ｔ波方向一致的低宽U波,其成因和生理意义目前尚不十分清楚＂。
本文注重于ＱRS波的检测，而在查阅一些文献资料以后，发现QＲS波的检测主要分为基于小波变换的心电信号ORS波检测与基于EMO与Mａrｒ小波变换的心电信号ORS波检测两种。
基于小波变换的心电信号OＲS波检测
小波变换可以分为连续小波变换（CＷT）、离散栅格小波变换（ＤＷT）和离散序列的小波变换（ＤSwT）。
信号x（ｔ）的小波变换定义式是:
其中
是基本小波又称母小波函数
是母小波经过移位和伸缩所
生的一组函数，称之为小波基函数,a是尺度因子,它实现对母小波函数的伸缩变换，b是时移变量,它实现对母小波函数的移位变换,以确定对信号分析的时间中心"
在连续小波变换中,a、b、ｔ均是连续变量，而在离散小波变换中，需对它们进行离散化,常取
当
时就称之为二进离散小波变换，然而取
时,在实际信号分析中有时显得尺度跳跃跨度太大,当希望尺度ａ在ａ>
O的范围内取任意值进行分析时就需要进行连续小波变换"
下面将根据心电信号的连续小波变换模极大值线检测和定位R波峰。
心电信号的Ｒ波峰是奇异点,而且它具有较大的幅度和较高的斜率等典型特征,根据基于小波变换的信号奇异性检测理论可知,每个R波的位置都对应于小波变换的模极大值的汇聚点,所以本算法首先对心电信号作连续小波变换并对信号按照心动周期进行分段，以便分别对一个心动周期内的波形进行奇异性分析，然后分别在每一个心动周期内检测模极大值点,它们的连线就是模极大值线由此确定R波的位置,并剔除李氏指数为负对应为噪声产生的模极值线以及应用不应期策略减少噪声干扰,提高检测准确率。
具体算法实现步骤如下：
（１）对给定的心电信号作连续小波变换,小波基选用Ｈａar小波,分解尺度a=3２分解后得到的小波系数可在一个尺度一时间平面上以灰度图的形式表示。
（2）对心电信号按心动周期进行分段,分段算法是首先对尺度一时间图按尺度ａ的方向进行累加,从而得到在尺度方向上小波变换的积分值随时间变化的曲线＂对于Ｈaaｒ小波而言，该曲线在Ｒ波之前有一个波峰,Ｒ波之后有一个波谷"
再分别选其正！
负极大值的一半作为正负闭值,对积分值随时间变化的曲线进行闭值化处理,并令大于正阂值的点为+l,小于负阂值的点为一1,在两者之间的点等于0,这样在每一个R波位置的之前就有一个+l,之后有一个一１，两者之间的区域为0"
把某一个一１位置和其后出现的第一个+l位置这一段数据的中点定为心动周期的分割点,从而实现了信号的分段，每一段都包括一个心动周期，而其R波在该段的中部。
（3）对每个心动周期段信号的尺度一时间图，分别找出在每一个尺度下的正的极大值点和负的极大值点,将其连成线得到正!
负模极大值线"
由每条正！
负模极大值线的斜率求出该点对应的李氏指数,根据李氏指数判据剔除李氏指数小于O和大于１所对应的模极值线。
（4）因为信号的连续小波变换的模极值线有可能出现中断现象,所以需对每一条正!
负模极大值线进行直线拟合,以分别求出它们在尺度a二０的时间位置,若在a=0时正。
负模极大值并不收敛于同一个点,则取二者的平均值作为R波的初步位置。
（5）在初步确定为R波的位置对应10ｍs时间范围内,检测原信号的极值点，并将其最终确定为Ｒ波位置。
（6）应用不应期判据"
由于心肌细胞除极化和复极化需要一个过程,存在一个绝对不应期,所以除了室颤和室扑外一般人的心率小于３0０次/分。
一个QRS波群产生以后,其后一定时间间隔内都不会出现另一个QRS波群,我们把这个时间间隔称为不应期"
本算法中的不应期设置为Zooms"
所以检测到一个Ｒ波后将其后Zooms内的模极值都忽略，这样可以避免很多由噪声干扰所引起的误检。
基于EMO与Mａrr小波变换的心电信号ＯRS波检测
针对常规的基于EMD的QRS波检测算法在信号存在严重高频干扰的情况下会出现较多错检导致检测准确率较低的问题,本文将基于离散小波变换的QＲS波检测算法与EMＤ方法相结合，提出一种基于EMD分解与Marr小波变换的心电信号QRS波检测新算法,来克服以上算法的不足,即尝试利用ＥMD分解法将非平稳心电信号分解为一系列具有不同特征尺度的IＭF分量,然后利用Marr小波变换对相应低阶IＭＦ分量叠加得到的重构信号进行奇异性分析,从而实现对原始心电信号QＲS波的准确检测和定位。
ＥMD分解:
EＭＤ分解的低阶本征模态分量中包含原信号的骤变部分，而高阶本征模态分量中包含缓变部分。
在心电信号中,对于高瞬时幅频的QRS波群自然就被分配到低阶高频模态分量中,而且Ｒ波的局部特征在第一、二本征模函数分量中得到了明显体现。
但EＭD算法中包含局部求极值！
样条插值!
边界效应处理等步骤,其计算量相当可观,使得处理速度非常缓慢,而且目前没有快速算法，因此无法满足实时动态检测的要求"
而且每分解出一个本征模函数分量,计算量将增大一倍,所以本文根据心电信号的时频特性和检测的实时性要求,提出只对心电信号作三层经验模式分解处理,然后将分解得到的第一、二、三本征模函数分量直接相加重构得到一个新信号,通过对此新信号进行奇异性分析来实现QRＳ波的检测和定位,这样不仅可以有效抑制基线漂移,高幅P波!
T波以及伪差信号等低频干扰以及边界效应,而且还能将处理速度提高几倍。
但是由第一、二、三模函数分量相加所构成的信号中往往还会包含QRS波带宽以外的频率分量,所以直接对它进行阂值判决的Ｒ波检测算法的正确检测率必然不高,而且容易受到高频噪声的干扰,抗干扰能力较差,但是把它作为定位R波的预处理信号是不错的选择"
另外EMＤ分解中筛选过程的中止准则常用方差，但也可根据信号特点手动设定筛选次数"
研究发现，筛选次数小,QＲS波在本征模函数域对应的分量越不明显;
而筛选次数越多,中心频率越大,特别是运算量成倍增长＂通过反复实验尝试，本研究通过对心电数据进行8次筛选,以极小的分解损失换取高的计算速度，而且丝毫不影响ＱRS波的提取效果。
小波基的选取
由前面的讨论可知,在基于离散小波变换的QRS检测中,定位算法及检测效果与小波基函数的选择密切相关,Marr小波（又称Ｍexicanhａt小波）具有良好的连续性、对称性以及指数衰减性,并且还具有一阶消失矩等性质,非常适合对信号进行奇异性检测。
Marr小波的母函数是高斯函数的二阶导数与常数的乘积,表达式为:
因为它像墨西哥帽的截面，所以也常称之为墨西哥帽小波。
Ｍarｒ小波函数属于二次微分小波,在时域和频域都有很好的局部化,并且满
由于Marr小波函数具有无限光滑性以及无穷次可微,并且不对单独的噪声点敏感，再加上其独特的时域性质,能使包含信息的特征点特别突出,因此本文选用Marr小波基进行Ｒ波峰值奇异点检测,应具有良好的定位特性和分析精度"
根据Marr小波基函数,计算得到相应的小波分解低通和高通滤波器的系数l和h,如下图2-1所示:
根据人和气就可以利用Mallａt算法递归计算出信号的小波变换。
图2-1
基于Marｒ小波变换的R波峰值奇异点定位
由前面的讨论可知,信号x（t）的所有奇异点在尺度一时间平面的模极大值线上,且其小波变换
在
充分接近于零时,其模极大值点就是信号的突变点。
由于Marｒ小波是二次微分小波，而且图形是以原点左右对称的,因此原始信号的奇异点在其小波变换的各层细节信号上仍然保持为极大值,这就使得对原始心电信号R波峰值奇异点的检测可以转化为对特征尺度上细节信号的极大值点的检测"
相比之下,Marr小波能克服采用一次微分小波检测信号奇异点时存在的以下缺陷:
（l）一次微分小波检测算法需通过检测小波模极大值对的过零点来定位信号奇异点,而过零点易受到噪声干扰，使得定位精度的稳定性难以保证。
（2）一次微分小波变换算法中需借助于一对相邻的模极值点位置及两者之间的斜率间接确定Ｒ波位置,并且还要根据特征尺度进行时移修正,其计算过程相对比较复杂和繁琐。
而我们以软件为主的方法实现QRS波的检测滤波之后的信号一般经过一些变换以提高QRS波的份量，进而采用一系列阈值进行判别,这些阈值有固定阈值法，也有可变阈值法。
前者由于可能的干扰或高P、高T波的存在，若其滤波后超过其阈值便会产生假阳性（ＦP,falseposiｔiｖe）结果；
另外,当心律失常或QRS波幅度变小，阈值设置过高,会导致漏检产生假阴性（FＮ，ｆalsenegativｅ）结果。
由于固定阈值的这些缺点,有研究者提出了用可变阈值检测,以提高检测的精确率，所采用的可变阈值包括幅度阈值、斜率阈值和时间间隔阈值等。
编写的源程序:
Ｑ波和S波通常是低幅高频波，一般Q波位于Ｓ波之前,S波位于R波之后,由于他们是一般向下的波,所以他们的峰值点和极值是对应的。
因次在检测到Ｒ波向左和向右分别搜寻到极值点,对应的就是Q波和S波。
而现在我们只需要检测到Ｒ波，所以就不需要检测Q波与S波的极值点了。
具体程序如下：
ECG-R波检测：
ｃlearall;
clｃ;
z=texｔrｅaｄ（'
ECG.txt＇）;
ＥCＧ＝ｚ（:
1）;
inｐuｔ＝ＥCG（１：
2５6）;
rate=ＥＣG（1００）;
ｓig=iｎｐｕt;
ｌｅnsiｇ=ｌength（sig）；
wtsiｇ1＝ｃｗｔ（siｇ，６,'
mexh＇）;
leｎwtｓiｇ1=lｅngth（wtsig1）;
wｔsig1（1：
20）=0;
wｔsig1（leｎｗtsig1-２0:
lｅnｗtsig1）=0;
y＝wtsig1;
yabs=abs（y）；
%?
sigtｅmp=y;
sｉgleｎ＝lｅngtｈ（ｙ）;
sｉgmax＝[］;
ｆｏｒｉ=1:
sｉgｌen-２
if（y（i+1）>
y（i）&
ｙ（ｉ＋1）>
ｙ（ｉ+2））｜（y（i＋1）＜y（i）＆ｙ（i+1）<
y（ｉ+2））
ｓigmaｘ=［sigmax;
abs（ｓｉgteｍp（i＋1））,i+1]；
end；
ｅnd;
%打印原信号及变换信号
fｉgｕrｅ（１）;
sｕｂplｏｔ（2,1，１），pｌot（ｓｉg）;
tｉｔlｅ（＇ECG原信号'
）;
sｕｂpｌot（2,1,2）,pｌｏt（wｔsｉg１）;
tｉtle（＇EＣG变换信号'
）；
%取阈值,阈值为相对幅值的差的6０%
ｔhrｔemp＝sort（sigmａx）;
ｔhｒlｅｎ=leｎgth（ｓigmax）；
thr=0;
foｒi=（thrlen-7）:
tｈrｌen
thr=thr＋thrｔemp（i）;
ｅｎd；
thrmax=thr/8;
%最大幅度平均值，8个最大幅值点的平均值
zeroｔeｍｐ=soｒt（y）;
zerｏvalue=０;
forｉ=1:
1０0
zerovａｌue＝zｅrovａluｅ+ｚeｒotemｐ（i）;
end;
zeｒｏvａｌｕe＝zerovaｌｕｅ/100;
%最小幅度平均值，对消幅度，100个最小幅值点的平均值
thr=（thｒmａx-zｅrovalｕe）*0.3;
%最大、最小幅度的差值的30％为判别R波的阈值
%定位R波
rvalue=[];
foｒi=１：
tｈｒｌｅｎ
ifsｉgmaｘ（i,1）>
thr
rvalue＝[ｒｖａlue;
sigmaｘ（i,2）]；
enｄ;
ｅｎd;
ｒvalue_1=rｖaluｅ;
%排除误检,如果相邻两个极大值间距小于０．4,则去掉幅度较小的一个
lｅnvalｕｅ=ｌength（ｒvalue）；
i=２;
wｈilei<
=lenvａluｅ
if （ｒvａlue（i）-ｒvalｕe（i-１））＊raｔｅ<
0.4
ifｙabs（rvaluｅ（ｉ））>
yabs（ｒvａlｕe（i-1））
rvalｕe（i－1）＝[];
elｓe
rvａlｕe（i）=[]；
end;
lenvalue=lｅｎgth（rvaｌuｅ）；
i=i-1;
i=i＋１；
lenvaluｅ=length（rvalｕe）；
%在原信号上精确校准
fori=１：
lenｖalue
if（wtｓig1（rvalue（i））>
0）
k＝（rvalｕｅ（ｉ）－５）:
（rvaｌue（i）+5）;
［a,b］=max（sｉｇ（k））;
rｖａlue（i）=rvalｕe（i）-6+b;
else
ｋ=（rvalue（i）-５）:
（rｖalue（ｉ）+５）;
[a,b]=min（siｇ（ｋ））;
rｖａlue（i）=rvalue（ｉ）-６+b;
%打印纠正及校准前后的R波信号
fｉｇｕｒe
（2）；
subplｏt（２,1,１）,pｌｏt（1：
lenｓig,wtsｉｇ１,rvalue_1,wtｓｉg1（ｒvalue_1），'
r.＇）;
ｔｉtle（'
EＣG纠正及校准前的R波信号'
ｓubpｌoｔ（2，1,２）,plｏｔ（1:
lensig,siｇ,rvalｕe,sｉｇ（rvalue），'
r．'
tiｔｌｅ（'
ＥＣG纠正及校准后的R波信号'
NOISYECG-R波检测:
cleaｒ all;
ｃlc;
ｚ＝textread（＇NOISＹECG.tｘt'
ECG=z（:
1）；
ｉnput=ECG（１:
２５6）;
rａｔe=ECG（10０）；
sig＝inpuｔ;
ｌｅnsig=lengｔh（sｉg）；
ｗtsｉg1=cｗt（sig,6,'
mexh'
lenwｔsｉg１=lｅnｇth（ｗtsｉg1）；
ｗｔｓｉｇ１（1:
wtsig1（lenwtsig１-2０：
lenｗtｓig1）=0；
y=wtsiｇ1;
yａbs=abs（y）;
%?
sigｔemp=y;
siglｅn=leｎgth（y）;
sigmaｘ=[］;
foｒｉ=1：
ｓiｇｌen-２
if （ｙ（ｉ+1）>
y（ｉ+1）＞y（i＋2））｜（y（i+1）＜ｙ（ｉ）&
ｙ（i+１）<
ｓigmａx＝[sigmax；
ａbｓ（ｓigtｅmp（i+1））,i＋1]；
enｄ；
eｎd;
％打印原信号及变换信号
ｆiguｒｅ
（1）;
subｐlｏt（2，1,1）,plot（sig）;
tiｔle（＇NOIＳYECG原信号＇）;
sｕbｐloｔ（２，１,2）,pｌot（wtsig１）;
titｌｅ（'
NＯＩSYＥＣG变换信号'
％取阈值,阈值为相对幅值的差的６0%
thrｔemp=soｒt（ｓigmax）;
thｒlen＝length（sigmax）；
tｈｒ=0;
foｒ i=（thｒlen-7）：
thｒlｅｎ
thr=tｈr+thrtemｐ（i）;
tｈrmax=thr/８;
%最大幅度平均值，８个最大幅值点的平均值
zeｒotｅｍp=ｓort（y）;
zｅrｏｖaluｅ=０;
ｆori=１:
100
zerｏｖalue=zerovalue+ｚerｏtemp（i）；
zeｒovaｌuｅ=zerｏvalue/10０;
%最小幅度平均值,对消幅度,１０0个最小幅值点的平均值
ｔhr＝（thｒmax-zeroｖalｕe）*０．３;
%最大、最小幅度的差值的30%为判别R波的阈值
rvalue＝[];
ｆor ｉ=1：
thrlen
ｉf sｉgmaｘ（i，1）>
rvａｌue=[rvalue;
sｉgmax（i，2）];
ｅnd；
enｄ;
rvａlue_1=rvａlue;
%排除误检，如果相邻两个极大值间距小于0.4，则去掉幅度较小的一个
lｅnvalue=lengtｈ（rｖalue）；
i＝2;
whｉlei<
=lｅnvａlue
ｉf（rvalue（i）-rｖａluｅ（i-１））*raｔｅ<
if ｙabs（rvaｌue（i））>
yabs（rｖａlue（i-１））
rvaｌuｅ（ｉ-1）＝[］;
ｅｌsｅ
ｒvalｕe（i）＝［];
end;
lenvalue＝ｌeｎｇtｈ（rｖａｌue）;
ｉ＝ｉ-1;
ｅnd;
i=i+1；
ｌｅnvaluｅ=length（ｒvaluｅ）;
fori=１:
lｅnvaluｅ
if（wtsiｇ1（rvaｌue（i））>
k＝（rｖalｕe（i）-5）:
（rvａｌue（ｉ）+5）;
[a,b]=ｍaｘ（ｓｉg（ｋ））;
ｒvalue（ｉ）=rvalue（i）-6＋b;
eｌse
k=（rvalue（ｉ）-５）:
（rvalue（ｉ）＋５）；
[a,ｂ]=ｍiｎ（siｇ（k））；
rvaｌuｅ（i）=rvalue（i）－6+ｂ;
ｅｎd;
%打印纠正及校准前后的Ｒ波信号
figurｅ
（2）;
ｓuｂplot（2，1,１），plｏｔ（１：
lensig,wtｓiｇ1,rvａluｅ_1,wtsig1（rｖalue_1）,'
title（'
ＮＯISＹＥCG纠正及校准前的Ｒ波信号'
subplot（２，1,２），ploｔ（1:
lｅnsｉｇ，sig,ｒｖａlｕｅ，sig（rvalue）,'
r．＇）;
ＮOISYECG纠正及校准后的R波信号'
结论（画出要求的图形）：
ﻩ
ＥＣＧ－R波检测所获得的结果如下图３-1与图４-1所示：
图3-１
图4－１
而NOＩSＹECG-R波检测所获得的结果如下图5-１与图6－1所示:
图5－1
图6-1
总结:
实习报告分数:
指导教师:
升级会员 