动态脉搏信号的采集与处理

2012-01-26 07:44张爱华丑永新
中国医疗器械杂志 2012年2期
关键词:陷波波谷波峰

【作 者】张爱华,丑永新

1 兰州理工大学 电气工程与信息工程学院,甘肃,兰州,730050

2 甘肃省工业过程先进控制重点实验室,甘肃,兰州,730050

近年来,随着人们生活方式、饮食结构的不断改变,诸如高血压、冠心病等疾病成为常见病,并且发病率不断上升,特别是心血管疾病已成为威胁人类生命的主要疾病。为了避免和减轻这些疾病对人类健康的危害,如可穿戴式心电监护仪,远程医疗监护系统等应运而生。文献[1]-[6]提到的脉搏信号检测系统,将采集的人体生理信号作一定的处理,然后通过有线或无线通信的方式传递给PC机,由医生诊断并给出相应的治疗措施。这样给患者带来许多方便,提高了医生治病的灵活性和有效性。但由于心血管疾病发病的随机性和快速性的特点,提高监护系统的实时性就显得十分重要。

脉搏信号蕴含了丰富的人体生理病理信息,且用光电脉搏传感器信号检测方便[7]。本文以TMS320VC5509DSP为核心,研制包括脉搏信号采集电路在内的动态脉搏信号实时检测与分析系统。针对动态脉搏信号的特点,设计了整系数陷波和低通滤波器,提出了包络滤波法。然后,通过优化编程方法,在DSP处理器上实现动态脉搏信号的实时滤波,对滤波后的信号进行时域和频域分析,计算脉搏信号的脉率、波形高度,并对脉搏信号进行FFT变换,实现光电脉搏信号的实时检测和简单处理。

1 整系数滤波器的设计

脉搏信号是随着心脏的搏动而产生的,一般比较微弱,特别是在不影响人体正常活动的条件下,采集的信号有许多来自外界的干扰。这些干扰主要有受呼吸和皮肤接触滑动影响产生的基线漂移,人体电位变化产生的肌电干扰,50 Hz交流电及其多次谐波引起的工频干扰等[8]。这些干扰会使得脉搏采集系统的信噪比下降,甚至会淹没有用的脉搏信号,仅靠信号采集系统硬件设计上的措施并不能完全消除这些干扰。因此,采取有效的手段消除或减少基线漂移、肌电干扰和工频干扰,是识别和处理脉搏信号的前提,在脉搏检测仪器中尤为重要。如文献[9]-[13],采用的IIR和FIR数字滤波等方法可以获得相当好的滤波效果,但这些滤波器系数一般为非整数,在实际应用中,特别是实时信号处理场合,对信号的实时滤波速度有很大的影响。本文针对脉搏信号的特征,设计了整系数陷波、低通滤波器,并结合包络法滤除动态脉搏信号中的干扰信号。

1.1 陷波器的设计

为了去除脉搏信号中的工频干扰,需要设计50 Hz的陷波器,同时为保证滤波前后信号不失真和便于DSP系统的实时滤波,该滤波器还有严格的线性相位且滤波计算量不宜太大的要求。因此,采用两个具有相同传输延迟的线性相位滤波器相减得到陷波滤波,即从一个全通网络中减去一个带通网络。带通网络可以用一个整系数的梳状滤波器实现,全通网络只具有延迟作用,起到抵消梳状带通网络带来的相位延迟作用[14]。系统的采样频率为400 Hz,则50陷波滤波器的传递函数和频率响应分别为:

式(1)中第一项为全通网络,后一项为梳状带通滤波网络;l、k为正整数,决定陷波器的阻带宽度和通带波纹。通过实验,选取l=12,k=2,则式(1)变为:

显然,式(3)前后两项的纯延时为96。由MATLAB仿真可得陷波滤波器的幅频特性曲线如图1(a)所示。由式(1)可求得该陷波滤波器的输出方程为:

1.2 低通滤波器的设计

由于经过陷波滤波器处理以后,脉搏信号中仍然存在较为严重肌电干扰,且其频率范围广泛,不能从原始信号中根本的消除,只能通过低通滤波器减少其对脉搏信号的干扰。为此,采用零极点对消法设计简单的整系数低通滤波器。该滤波器的传递函数及频率响应如下:

由式(7)可以看出,此低通滤波器有M个零点和极点。其中,M个零点均匀的分布在单位圆上,在z=1处有一个极点,用来抵消第一个零点,其他的极点都集中在原点处,因而具有低通特性。整个单位圆圆周角取值范围是0-2π,实际采样频率为fs,其对应的圆周频率为0~2 πfs。为了使得低通截止频率相对较小,假设幅频特性曲线第一次降为零的频率为2πf(0)/fs,所有零点的个数为M=fs/f(0)。由于DSP系统采样频率为400 Hz,且必须保证M为整数,取f(0)=66.67 Hz,得到M=6,则传递函数为:

由MATLAB仿真可得低通滤波器的幅频特性曲线如图1(b)所示。其第一次降为零点时的截至频率为66.7 Hz,在通带范围内,滤波器具有线性相位,满足滤波要求。由式(9)可求得该低通滤波器的输出差分方程为:

图1 滤波器的幅频特性曲线Fig.1 Magnitude response of fi lters

1.3 包络法去基线漂移

基线漂移属超低频信号,可以通过高通滤波器从原始的脉搏信号中消除。但在高频滤波过程中,或多或少会使脉搏信号中的有效成份丢失。为此,提出一种新颖的滤波方法—包络法,去除脉搏信号中的基线漂移。

如图2所示,上、下两包络线是脉搏信号的峰值和谷值经过多项式曲线拟合法拟合的曲线,中线为峰值和谷值包络曲线的均值曲线,可以近似反映脉搏信号的基线漂移。用原始的脉搏信号减去均值曲线,即可去除脉搏信号中的基线漂移。

图2 脉搏信号包络示意图Fig.2 Envelope diagram of pulse signal

图2所示的脉搏包络曲线是在PC机上仿真得到,在检测出脉搏信号的峰值和谷值之后,需要用多项式曲线拟合法拟合出脉搏信号的峰值和谷值包络。但由于DSP系统运行速度的限制,很难用多项式曲线拟合法实时地拟合出脉搏信号的包络曲线。我们用相邻两峰值(谷值)之间的连线近似的代替拟合曲线,作为脉搏信号的包络曲线,这样极大地降低了算法的计算量,可以在DSP系统内两采样点之间实时地实现包络曲线的估算和基线漂移的去除。

2 脉搏信号的时域和频域分析

2.1 脉搏信号的时域分析

时域分析主要是指在时域提取和分析脉搏特征,如幅值特征、时间特征等。本文主要分析脉搏信号的波高和脉率。

2.1.1 波高检测

波高是指脉搏波波峰和相邻波谷幅度的差值,要计算波高,首先要检测脉搏波波峰和波谷。因为DSP系统内设置的数据缓冲区宽度为1024点,只能存储1024点的脉搏数据,我们只需要在这1024点数据中求出脉搏波的波峰(波谷)即可。由于脉搏波波峰(波谷)为其邻域内最大值(最小值),同时为了避免潮波和重搏波的影响,通过对大量脉搏波进行反复分析,可以得到此邻域的最佳宽度为100点。于是,在数据缓冲区内,取100点宽度为一个固定窗,设为W[j],然后将此窗固定在数据缓冲区内,让数据流过此窗,求出此窗内数据的最大值(最小值),再根据脉搏信号波峰(波谷)判决条件判断此最大值(最小值)是否为波峰(波谷)。设脉搏信号数据序列为X[i],其中i∈[0,1024],由上面分析可得到在波峰邻域内波峰X[ip]的判决条件为:

式中max{W[j]}为固定窗内数据最大值,ip为波峰位置。由于最大值邻域宽度的限制,ip的取值不能小于50,不能大于1153,由此决定j的取值上限和下限。j的取值决定着固定窗在数据缓冲区内的位置,一旦 j值取定,就可以在数据缓冲区内确定一个长度为100点的固定窗。而k的取值就是为了限制j的取值,使其满足要求,k可以为其定义域内的任意整数。与波峰判决条件类似,可得波谷X[il]的判决条件为:

由上面方法求得的相邻波峰与波谷幅度差便可获得波形高度,即一个周期的波形高度是由此周期内的波峰幅值减去其前面最邻近的波谷幅值。

2.1.2 脉率计算

脉率是指每分钟脉搏的次数,正常情况下与心率一致,健康成年人在安静状态下脉率为每分钟60-100次。且脉率受诸如年龄、性别、体型等因素的影响,在一定的范围内波动。而脉率的不同,其对应的脉象也不同,反映着人体不同的生理和病理特征,故在临床诊断中有必要对脉率进行实时检测。

检测脉率的思路为:先检测相邻两个脉搏波的波峰(波谷)值;然后求出两个脉搏波峰(波谷)之间的数据点数(设为n);由采样率求出这n点所对应的时间,即为脉搏信号的周期(设为T);最后用60/T便可计算出实时脉率(设为N次/分)。具体检测脉率的算法如下:

1) 脉搏信号峰值检测并确定峰值位置。用式(11-1)(式(11-2))检测出两个连续周期波峰(波谷)及其位置,并求出波峰(波谷)之间的点数n。

2) 确定脉搏信号的周期T,计算脉率值。采样频率fs =400 Hz,则采样周期为1/fs,故脉搏信号的周期为T=n/fs。根据定义,可求得脉率N=60/T=(60×fs)/n。

2.2 脉搏信号的频域分析

虽然用时域法研究脉搏信息比较简单直观,但用其在做脉象和脉势等方面研究已经不能完全满足要求,有一定的局限性。脉搏信号有着自身独特的频域特征,这些频域特征蕴含着丰富的脉搏信息,求出脉搏信号的幅度谱十分必要。对信号进行离散傅里叶变换(DFT),其基本变换公式为:

式(13)中,每项求和分别是原始序列的N/2个偶数点和N/2个奇数点的DFT。相比DFT,计算量极大的减少。可应用式(13)对DSP系统数据缓冲区内的128点脉搏数据进行FFT变换,求出其频谱图。

3 基于DSP的动态脉搏信号采集与处理

为了实时采集和处理动态脉搏信号,以TI公司生产的TMS320VC5509 DSP芯片为核心,设计了动态脉搏信号的实时检测系统,如图3所示。此系统主要由三个模块组成:1)光电脉搏传感器;2)脉搏信号调理电路;3)DSP实时检测与处理系统。光电脉搏传感器采集人体的生理脉搏信号,经信号调理电路,被送入DSP系统进一步做信号滤波、时域和频域分析及结果显示,实现脉搏信号的实时检测和分析。

图3 脉搏信号检测系统结构图Fig.3 Architecture diagram of pulse signal detection system

3.1光电脉搏传感器及脉搏信号的调理电路

光电脉搏传感器主要完成生理脉搏信号向电信号的转换。传感器内置一个近红光(940 nm)发光二极管,用来提供光源;一个光敏二极管,用来感知光信号。由于动脉血液对光的吸收量随动脉搏动而变化,光敏二极管可将变化的光信号转化为变化的电流信号,实现将指尖脉搏的波动转化为电流信号输出。

脉搏信号的调理电路实现动态脉搏信号的初步调理,先通过I/V转换电路,将光电脉搏传感器拾取的电流信号转换为电压信号;然后,通过一个截止频率为0.16 Hz二阶高通模拟滤波器,去除基线漂移;再由一个截止频率为80 Hz的二阶低通模拟滤波器滤波,降低肌电干扰。在滤波的同时,通过两级放大器,将微弱的脉搏信号放大500倍,转换为可供DSP内ADC转换器采样的电压信号。

3.2 DSP实时检测与处理系统

脉搏动态信号的检测与处理系统其核心为TMS320VC5509DSP处理芯片,具有高处理速度(400MIPS)的定点内核DSP芯片。自带RAM(128 K×16 bit), 片上FLASH(4 M×16 bit),2路10bit A/D转换器,最大采样率为12.5 KHz,采集电压范围-3V-+3V。该系统主要完成脉搏信号的AD转换、整系数滤波、波高和脉率计算、频谱计算,日期、时钟和脉率显示,并对脉率进行监控。当脉率异常时,声光报警。

图4 系统程序流程图Fig.4 Program fl ow chart

系统程序流程图如图4所示。设计思路为:采用DSP处理器自带的定时器通过定时中断的方式控制AD转换器采集数据,每采集一个数据,就在定时中断子程序内完成一次数据处理。图4(a)为程序主流程图,主程序主要完成系统初始化、计算脉搏信号的频谱、显示数据处理结果。图4(b)为AD转换定时中断子程序流程图。本设计采取定时中断方式控制AD转换,采样率为400 Hz。每次定时结束,定时器0便会向系统申请中断,系统响应中断以后,开始运行中断子程序。中断子程序主要完成工作为:首先启动AD转换器开始转换数据,读取并存储数据;然后对数据缓冲区的脉搏数据进行陷波、低通和包络滤波,检测脉搏信号的峰值和谷值,计算波高和脉率,并判断脉率是否异常,若脉率异常,发出声光报警;最后返回主程序。

3.2 滤波算法在DSP上的编程简化

在实际应用中,影响系统程序运行速度除硬件设备的因素外,另一主要的因素就是算法内乘法的运算次数。虽然DSP处理器自带硬件乘法器,但是采取有效的措施简化程序,减少运算过程中乘法次数和循环次数,对于实现信号的实时处理和充分利用DSP内部的硬件资源,有十分重要的意义。

如图4(b)所示,在AD转换的子流程图中,要在两个采样点之间绝对实时地实现信号的陷波滤波、低通滤波、波高计算和脉率计算等,计算量非常大。为了解决这个问题,我们对编程方法作如下改进:

1)将乘法和除法运算改为移位运算。如式(6)所示,运算中的乘法全为×2运算,可优化为左移运算。

2)采用单点滤波。在DSP处理系统中,每采集一个数据点,就要对数据缓冲区内的所有数据进行移位,进行数据的更新,然后再对数据进行滤波。若采用全数据滤波,由式(3)和式(9)可知,每次滤波都需要对系统的群时延(206个数据点)进行初始化。为了解决这个问题,每次在滤波结束时都要保留206个数据,用来进行下次滤波的初始化,于是每次参与滤波的数据为1024+206=1230个。这不仅增加了乘法和循环的次数,而且还占用了多余的存储空间。若采用单点滤波,即每采集一个数据点,只对此数据点进行滤波运算,如此所有的滤波过程只需要初始化群时延一次,同时降低了计算量和存储空间的占有量。

4 结果

图5(b)为经过陷波和低通滤波以后的脉搏信号图,对比图5(a)原始脉搏信号可以看出,经过陷波器和低通滤波器后,脉搏信号中的高频干扰得到了有效的抑制。图5(c)为经包络滤波后的脉搏信号,其幅值由400 mV~750 mV降为-158 mV~158 mV,基线漂移得到了有效地去除。 经过滤波处理后的脉搏信号,基本上能满足实时信号处理要求。图6为将脉搏信号进行FFT变换后的脉搏信号频谱图,由图看出脉搏信号的低频和高频干扰得到了有效地抑制,方便脉搏信号的频域特征分析。

图5 脉搏信号预处理Fig.5 Pulse signal preprocessing

5 结论

系统运行结果表明,所设计的整系数滤波和包络滤波算法在保存原始有用信息的前提下,有效地抑制了动态脉搏信号的噪声干扰。在此基础上,分别从时域和频域对所采集的脉搏信号进行了实时处理,计算出了脉搏信号的波形高度和脉率,以及FFT变换求得其频谱图。通过特征分析,实时显示人体状态和参数,并可实现异常状态报警。该系统在健康监护领域具有较强的实用性和应用前景。

本课题组对脉搏信号中的高频信号有了一定的研究。为了保存脉搏信号中的有用高频信号, DSP处理系统采用400 Hz作为系统的采样率,实时完成了脉搏信号的检测和处理。本设计所提出的滤波算法和脉搏信号的处理算法,同样可以在低采样率的系统中实现动态脉搏信号的实时检测和处理。

图6 脉搏信号频谱图Fig.6 Spectrum of pulse signal

[1] 于洋, 刘静. 手机无线心电监测技术系统的实现及性能测评[J].中国医疗器械杂志, 2010, 32(5): 391-395.

[2] 孟淑婷,王磊. 基于DSP的家庭健康监护仪的设计[J]. 国际生物医学工程杂志, 2010, 22(9): 280-282.

[3] Woosik Shin, Yong Dae Cha, Gilwon Yoon. ECG/PPG Integer signal processing for a ubiquitous health monitoring system[J].Journal of Medical Systems, 2009, 34: 891-898.

[4] Cheng Wen, Ming-Feng Yeh, Kuang-Chiung Chang, et al. Realtime ECG telemonitoring system design with mobile phone platform[J]. Measurement, 2007, 41(2008): 891-898.

[5] Sera fi m Tabakov, Ivo Iliev, Vessela Krasteva. Online digital fi lter and qrs detector applicable in low resource ECG Monitoring Systems[J]. Annals of Biomedical Engineering, 2008, 36(11):1805-1815.

[6] E Pinheiro, O Postolache. Non-Intrusive device for real-time circulatory system assessment with advanced signal processing capabilities [J]. Measurement Science Review, 2010, 10(5): 166-175.

[7] 张珣. 光电脉搏传感器的设计与改进[J]. 中国医疗器械杂志,2009, 33(5): 344-346.

[8] 杨颖飞. 强背景噪声下的脉搏信号处理算法研究[D]. 西安:西安电子科技大学, 2009.

[9] Jinseok Lee, Ki H Chon. An autoregressive model-based particle fi ltering algorithms for extraction of respiratory rates as high as 90 breaths per minute from pulse oximeter[J]. IEEE Transactions on Biomedical Engineering, 2010, 57(9): 2158-2167.

[10] Woosik Shin, Yong Dae Cha, Gilwon Yoon. ECG/PPG Integer Signal Processing for a Ubiquitous Health Monitoring System[J].Journal of Medical Systems, 2009, 34: 891-898.

[11] Kristjan Pilt, Kalju Meigas, Rain Ferenets, et al.Photoplethysmographic signal processing using adaptive sum comb filter for pulse delay measurement[J]. Estonian Journal of Engineering, 2010, 16(1): 78-94

[12] K Shafqat, D P Jones, R m langford, et al. Filtering techniques for the removal of ventilator artifact in oesophageal pulse oximetry[J].Med Bio Eng Comput, 2006, 44:729-737.

[13] S Poornachandra, N Kumaravel. A novel method for the elimination of power line frequency in ECG signal using hyper shrinkage function [J]. Digital Signal Processing, 2007, 18(2008): 116-126.

[14] 兰瑞芬, 胡广书. 高采样率下简单整系数工频陷波器的设计[J].航天医学与医学工程, 2008, 21(2): 152-156.

[15] 赵洪亮, 卜凡亮, 黄鹤松, 等. TMS320C5X应用系统设计[M]. 北京: 北京航空航天大学出版社, 2008

猜你喜欢
陷波波谷波峰
炮制工程骗钱的“甲方”
板厚与波高对波纹钢管涵受力性能影响分析
梅缘稻
波峰焊接技术现状及绿色化设计方向
中国、英国、美国、日本规范关于直墙波谷力计算方法的对比
作用于直立堤墙与桩柱的波峰高度分析计算
基于数字递归陷波的多通道瞬变电磁法周期噪声去除研究
频域陷波对直接序列扩频信号接收性能影响分析
中空玻璃胶接结构界面脱粘缺陷的超声与X射线检测研究
双陷波补偿算法在火箭推力矢量控制中的扩展应用