王静,张大慧,张宏
(浙江大学生物医学工程与仪器科学学院生物医学工程系,浙江杭州310027)
心血管系统疾病是当今危害人类健康的主要疾病之一。心电信号的实时分析是诊断心血管疾病的重要方法。QRS复合波的检测是心电自动分析的关键环节。与其他波形分量相比,QRS波具有较高的幅值和特定的频率范围,典型QRS波群的频率分量约在10 Hz~25 Hz之间[1]。心电信号是一种不平稳信号,并受到基线漂移(baseline wander)、肌电噪声(muscle artifact)、电极移动噪声(electrode motion artifact)等多种因素的影响,因此,QRS波检测必须克服各种干扰的影响,才能正确提取出波形特征。
本设计采用离散小波变换算法,其准确率高,但算法计算复杂,对处理器的性能有较高要求。FPGA的应用解决了这一难题。FPGA具有全软件和全硬件的用户可定制性和可重配置性,可以根据需要灵活改变FPGA中构成DSP系统的硬件结构来改变系统的功能。
系统采用离散小波变换(DWT)算法,将心电信号分解到不同的尺度,其实质是多尺度的带通滤波。本文采用二次样条小波(Quadratic Spline Wavelet)作为母小波,基于二次样条小波的小波变换离散的低通(H(ω))和高通(G(ω))滤波器分别是[2]:
可以看出,二次样条小波对应的系数简单,可以用2n(n∈整数)的多项式表示。在FPGA中,这类运算通过移位和加法实现,不必使用资源耗费量更大的乘法器[3]。
人体心电信号通过合适的电极经过右腿驱动电路进入系统,经放大、滤波等环节之后,通过ADC转换成数字信号进入FPGA,进行数字信号处理(DSP)。DSP部分由FIFO数据存储、DWT变换、阈值检测、单尺度过零点(Zero Crossing)检测、检测合成、心率及R波检测、结果显示等7个模块构成,整个DSP部分由状态机(State Machine)实现系统控制。信号处理算法通过FPGA开发工具DSP Builder实现。FPGA片上数字信号处理结构如图1所示。
FIFO用来缓存由ADC采集的原始心电数据。FIFO的读写时钟独立设置,写时钟频率与ADC采样率相同,为360 Hz;读时钟与DSP处理速度相适应,为50 MHz。设置不同的FIFO读写时钟,使得AD采样部分和后续的处理部分可以按照各自时钟独立进行,不受干扰,保证了系统的实时性。
心电信号频率范围通常在0.05 Hz~80 Hz[4]。根据ECG信号、噪声和伪信号的能量谱分布和表1可知,QRS复合波分量主要分布于第1~4尺度上,高于尺度4的分量中,QRS波的分量显著减少,同时运动伪迹和噪声显著增加。
表1 不同尺度对应滤波器组的3分贝带宽
设计采用尺度1~4进行数据处理。原始信号经过4组高、低通滤波器组进行滤波,得到信号的1~4尺度的频率分量,如图2所示。
如图3所示,原始波形是一组类QRS的波形,WF1~WF4分别是对原始信号进行1~4尺度小波变换得到的波形[5]。可见,原始波形的突变和回落会在4个尺度上产生一对正极大值-负极小值对(简称极值对)。利用极值对可以成功地将QRS波和肌电、工频、基漂等干扰区别开,尤其是对于由信号饱和产生的伪迹,由于其只能产生孤立的模极大值列,因而能有效地将其剔除。由于患者在同一尺度的不同时刻和不同尺度上的极大值和极小值都不一定相同,所以为检测所设定的阈值也是动态的,J时刻的最大值、最小值阈值Tmax[J]、Tmin[J]为含当前最值M[J]和前一次最值M[J-1]的多项式:
由图3看出,R波发生的位置恰好为4个尺度的极大、极小值之间的过零点位置,因此,用1标示出这个过零点位置,其他位置均为0。过零点会出现在两种情况中:(1)一个极大值紧跟一个极小值;(2)一个极小值紧跟一个极大值。两种情况出现其一即认为该位置存在过零点。当然,如果信号被高频噪声干扰,则在低尺度上噪声和无关信号量影响较大,这时可以先从频率较低的较大尺度(第四尺度)上找到极值对,然后在其附近逐级向较小尺度上搜索极值对,检测过零点并最终获得4个由0、1组成的序列。1~4尺度DWT变换及各自的过零点检测结果如图4所示。
有研究者认为,只在尺度2或尺度3上进行过零点检测即能代表QRS波的位置[6],但是为了尽可能减少因噪声干扰而产生的假阳性点,本设计中考虑在尺度1~4中,至少在3个尺度上均存在过零点,才表示这个位置可能有一个QRS波出现,其他情况均为噪声影响。因此,需要对3个尺度的过零点进行与(and)操作。从图3中可以看出,随着尺度数增加,过零点的位置有所延迟,且高尺度的信号被延迟更多。理论上,如果原始波形是对称的,则尺度j处过零点相对于原始波形延迟2j-1-1个周期(处理时钟)[2],不对称波形(如图3中原始波形的第二个峰值处)带来的延迟比对称波形多。本设计中根据理论值和实验情况,得出各尺度的延迟如表2所示。其中理论延迟为原始波形对称的情况下产生的延迟,系统延迟包括波形不对称造成的延迟时间增加值,总延迟时钟周期数为理论延迟和系统延迟的和。
表2 小波变换尺度1~4的延迟周期数
由于延迟的存在,要对4个尺度的0、1序列进行补偿。补偿过后的每3个序列进行与操作,所得结果再进行或(or)操作,得到一组0、1序列,1的位置即为信号中R波的位置,该序列配合心电采样率即可获得实验者的心率。为提高算法的准确率,此处进行了序列优化。考虑到人的心跳频率不会超过300次/min,即在0.2 s时间窗内不会有2个“1”出现,因此如果检测到在0.2 s时间窗内出现了多于1个“1”,则保留第一次出现的“1”,其余位置均为0,借此去除大部分假阳性点。
DSP Builder状态机根据FIFO中存储数据的个数来实现系统的控制。具体控制方法为:初始状态为NULL;当FIFO数据个数大于某个设定值N(由公式5可知,N至少为3)时,进入READ状态,读取FIFO中的信号并使能后续模块,进行信号处理;当FIFO中的数据个数等于0时,进入READY状态,禁止后续模块以节省能耗,同时等待接收数据。当FIFO中数据个数大于0时再次进入READ状态,进行下一轮处理。由于显示模块需要随时显示结果,因此不被状态机控制,也不会进入节能状态。
采用MIT-BIH数据库中第一组23个信号中的第一通道信号对算法进行检测,每组信号截取90 000个采样点约250 s信号进行实验。表3列出了检测出假阳性、假阴性点的8个检测结果,其余15个文件中检测准确率均为100%,表中未列出。表3中文件104有较多的肌电干扰,导致假阳性点较多;文件108患者有窦性心律失常现象,基线漂移和噪声也较严重。对MIT-BIH数据库的23组信号进行检测,准确率高达99.5%,验证了本方法的准确性和可用性。
表3 MIT-BIH数据库检测结果
图5记录了6组特殊原始信号的检测结果。除图5(f)之外,其他信号都是采用MLⅡ(modified limb leadⅡ)导联。从图5中可以看出,本系统所应用的算法对含有基线漂移、噪声、某种特定心脏疾病以及MLⅡ导联和modified lead V5导联的心电信号,都能够正确检测出其QRS复合波的位置,具有较高的鲁棒性。
图5(a)~图5(f)的6组图中,每组2个波形中,上面的波形表示原始数据,下面幅值为1的脉冲表明QRS复合波的位置。图5(a)信号取自文件121,有较严重的基线漂移;图5(b)信号取自文件119,采集者有室性三联律现象;图5(c)信号取自文件114,信号被取反;图5(d)信号取自文件111,信号中的噪声较严重;图5(e)信号取自文件108,有尖锐的P波;图5(f)信号取自文件104,采用了与其他信号不同的校正的V5导联。
本文所设计的基于FPGA的实时QRS波检测系统,延时T主要由3个部分组成:基于采样频率的采样延时Ts、基于FPGA时钟的信号处理延时Tp和显示延时Td。如式(5)所示:
由以上分析可知,采样延时Ts是延时的主要因素。如果要进一步提高系统的实时性,可以适当提高信号的采样频率(各尺度的3 db带宽会随采样率变化)。但从理论上讲,5.56 ms的延时人眼基本无法分辨,360 Hz的采样频率对于QRS波检测的实时系统已经足够了。
长期以来,小波变换算法的复杂性所带来的大运算量和高成本限制了其在实时数字信号处理领域的广泛应用,FPGA的应用解决了这一难题。FPGA构成的DSP模块采用硬件逻辑实现,节省了软件系统中处理器执行指令流所消耗的时间。一个采样点的信号处理在几个处理时钟周期中完成,而采样时钟和处理时钟独立的结构,将信号采样和信号处理过程区分开来,使得信号处理过程的延时能够完全控制在一个采样周期之内,为心电信号的实时检测提供了保证。此外,利用FPGA内部丰富的硬件资源和硬件可定制性,将整个信号处理系统在一个FPGA芯片上实现,不但减小了电路的复杂度,降低了成本,方便了系统设计和调试,且提高了系统的整体性能。利用FPGA的可重构特性,通过修改FPGA中的程序,能快速、方便地实现功能的升级、更新,能适应不同应用对处理性能和结构的不同要求,为复杂数字信号处理提供了一个新的解决方案。
[1]KOHLER B U,HENNIG C,ORGLMEISTER R.The principles of software QRS detection-reviewing and comparing algorithms for detecting this important ECG waveform[J].IEEE Engineering in medicine and biology,2002,21(1):42-57.
[2]LI Cui Wei,ZHENG Chong Xun,TAI Chang Feng.Detection of ECG characteristic points using wavelet transforms[J].IEEE Transactions on biomedical engineering,1995,(41):21-28.
[3]Al-haj A M.An FPGA-based parallel distributed arithmetic implementation of the 1-D discrete wavelet transform[J].Informatica,2005,(29):241-247.
[4]王保华.生物医学测量与仪器[M].上海:复旦大学出版社,2003.
[5]陈永利.动态心电自动分析中QRS复合波检测算法研究[D].杭州:浙江大学,2006:21-36.
[6]熊沁,方祖祥,宋海浪,等.基于小波变换的QRS波实时检测方法与实现[J].中国医疗器械,2007,31(4):242-244.