CN105832323A - 一种基于rc电路的光电血管容积图信号处理方法 - Google Patents

Abstract

本发明公开了一种基于RC电路的光电血管容积图信号处理方法，包括高通滤波、低通滤波、信号平滑化和信号求导四个部分。本发明基于RC电路理论，利用RC电路硬件特性，提出了一种低复杂度且可靠的PPG信号处理方法来分析和处理脉搏波，并计算心率，不仅能准确地计算出心率，同时也大大节约了能源，即便使用纽扣电池，使用本发明的脉搏血氧仪也可以维持数月之久。

Description

一种基于RC电路的光电血管容积图信号处理方法

技术领域

本发明属于信号处理技术领域，具体涉及一种基于RC电路的光电血管容积图信号处理方法的设计。

背景技术

光电血管容积图(PPG)信号可以描绘出人体丰富的心血管疾病和呼吸系统疾病的物理信息。光通过生物组织后可以被皮肤色素、骨骼、动脉和静脉血液等不同的物质吸收。但是大部分吸收的光在动脉和小动脉的血流中变化剧烈，例如在心动周期的不同阶段。为了获得一个人的心率(HR)和血红蛋白的血氧饱和度(SpO₂)，光源常常为红光或红外光。红光或红外光通过组织传输，由于受到皮肤反射以及肌肉、骨骼、脂肪、静脉和动脉血液等散射，在检测器接收到的光将缓慢衰减。

脉搏血氧仪是用于测量一个人的血液氧饱和度、心率和呼吸速率的一种非侵入性装置。基于PPG信号的医疗设备需要干净且增强的信号用于特征提取、分析和检测，然而从设备输出的PPG信号是易受噪声污染的。许多医疗设备，如运动手环、移动手表等，它们的物理参数需要在人们日常活动中进行监控。通常，运动伪影是劣化PPG信号的主要噪声。运动伪影造成的不准确的PPG信号读数和解释不仅会导致医疗成本增加，还会使病人的治疗变得效率低下甚至造成危险。

在近期的研究中对于运动伪影的减少方法可以大致分为时域和频域处理的方法，时域处理方法的计算过程更加高效。其中有文献基于提取人类脉搏波形的特征并与衍生方案相结合，产生了基于生理规则选择的PPG信号和分离的有用信号。由于PPG信号的运动伪影干扰主要由随机的低频干扰组成，所以运动伪影的减少由信号处理的频域来实现。另有文献简要回顾了电流信号处理中的使PPG运动伪影减少时使用的方法。

对于嵌入式手腕脉搏血氧仪，电源一般由纽扣电池或小型锂电池提供。由信号处理算法所消耗的能量大约占用了一半的总耗能。常规的PPG信号处理方法复杂度高，使得以电池供电的脉搏血氧仪使用时间非常有限，给人们日常心率、脉搏监测带来很多不便。为了实现腕部装置的低功率和长期监测的目标，一个简单的、低复杂度的和有效的PPG信号处理算法是一个挑战。

发明内容

本发明的目的是为了解决现有技术中常规的光电血管容积图信号处理方法复杂度高，使得以电池供电的脉搏血氧仪使用时间非常有限，给人们日常心率、脉搏监测带来很多不便的问题，提出了一种基于RC电路的光电血管容积图信号处理方法。

本发明的技术方案为：一种基于RC电路的光电血管容积图信号处理方法，包括以下步骤：

S1、设计高通滤波器，对光电血管容积图信号进行高通滤波；

S2、设计低通滤波器，对高通滤波后的信号进行低通滤波；

S3、对经过高通滤波和低通滤波后的信号进行平滑滤波处理，使信号平滑化；

S4、对平滑滤波处理后的信号求导，以更准确地计算心率。

进一步地，步骤S1具体为：

导出一阶RC高通滤波器，设输入序列和输出序列分别为(x₁,x₂,…,x_n)和(y₁,y₂,…,y_n)，则有：

y_i＝αy_i-1+α(x_i-x_i-1) (1)

式中，τ为RC的时间常数，一阶RC高通滤波器的截止频率f_c满足关系f_c＝1/(2πτ)；ΔT为输入序列和输出序列的产生时间间隔，i＝1,2,3…n。

进一步地，步骤S2具体为：

导出一阶RC低通滤波器，设输入序列和输出序列分别为(y₁,y₂,…,y_n)和(z₁,z₂,…,z_n)，则有：

z_i＝βy_i+(1-β)z_i-1 (2)

式中，τ为RC的时间常数，一阶RC低通滤波器的截止频率f_e满足关系f_e＝1/(2πτ)；ΔT为输入序列和输出序列的产生时间间隔，i＝1,2,3…n。

进一步地，步骤S3具体为：

设输入序列和输出序列分别为(z₁,z₂,…,z_n)和(v₁,v₂,…,v_n)，则有：

v_i+width/2-1＝sum(z_i,z_i+1,...z_width+i-1)/width (3)

其中width为平滑宽度，i＝1,2,…,n-width。

进一步地，步骤S4具体为：

令x＝(z₁,z₂,…,z_n)，y＝(v₁,v₂,…,v_n)，则通过中央差分法对平滑滤波处理后的信号求导：

$\left\{\begin{array}{c}{y}_j^{\&prime;}=\frac{y_{j+1}-y_{j-1}}{2\mathrm{\&Delta;}x}\\ \begin{array}{cc}{x}_j^{\&prime;}=x_j,& 2<j<n-1\end{array}\end{array}\right.---\left(4\right)$

式中x'_j和y'_j是第j个点的x和y值的导数，Δx为相邻x点的差值；根据公式(4)检测到光电血管容积图信号的峰值，并且通过计算光电血管容积图信号的峰值的间距得到心率。

本发明的有益效果是：本发明提出了一种低复杂度且可靠的PPG信号处理方法来分析和处理脉搏波，并计算心率，不仅能准确地计算出心率，同时也大大节约了能源，即便使用纽扣电池，使用本发明的脉搏血氧仪也可以维持数月之久。

附图说明

图1为本发明提供的一种基于RC电路的光电血管容积图信号处理方法流程图。

图2为本发明实施例的一阶RC高通滤波器示意图。

图3为本发明实施例的一阶RC低通滤波器示意图。

图4为本发明实施例的原始PPG信号的波形示意图。

图5为本发明实施例的原始PPG信号的FFT波形示意图。

图6为本发明实施例的原始PPG信号分别通过高通滤波器，低通滤波器，平滑滤波器的波形示意图。

图7为本发明实施例的原始PPG信号分别通过高通滤波器，低通滤波器，平滑滤波器的FFT波形示意图

图8为本发明实施例的静止和运动状态的PPG采样信号波形示意图。

图9为本发明实施例的真实心率与本发明的计算结果的比较示意图。

具体实施方式

下面结合附图对本发明的实施例作进一步的说明。

根据PPG信号的特点，通常有基线漂移、运动伪影、电源线接口和环境光噪声等污染。通常一个人的心率范围从30bpm到240bpm，其对应频率为0.5Hz至4Hz的PPG信号。因此使用过滤设备的PPG信号算法一般集中在这个频率范围。本发明提供的PPG信号处理方法是基于传统的过滤思想，其中使用了高通和低通滤波器。

本发明基于RC电路理论，利用RC电路硬件特性，提供了一种基于RC电路的光电血管容积图信号处理方法，如图1所示，包括以下步骤：

S1、设计高通滤波器，对光电血管容积图信号进行高通滤波。

高通滤波器的目的是过滤比0.5Hz低的基线漂移频率，我们设计的高通滤波器是基于RC电路理论，结构简单，适合于嵌入式实时进程。在通常情况下，使用第一或第二阶高通滤波器。本发明实施例中，为了简单起见，我们在这里仅导出一阶高通滤波器，如图2所示。

结合基尔霍夫定律和电容器的充电和放电原理，可以建立下面的公式：

$\left\{\begin{array}{c}{u}_o\left(t\right)=i\left(t\right)R\\ Q_c\left(t\right)=C\left(u_i\right(t)-u_o(t\left)\right)\\ i\left(t\right)=\frac{{\mathrm{dQ}}_c\left(t\right)}{dt}\end{array}\right.---\left(5\right)$

其中，u_o(t)是脉搏血氧仪基于图2所示的一阶RC高通滤波器电路，在时间t，电阻R上流经电流i(t)时的电压输出，输入电压为u_i(t)；Q_c(t)是电容C在时间t的电量；依据电容充放电原理，流经电容C的电流i(t)为电容C上电量的导数dQ_c(t)/dt。

根据这三个方程，我们可以得到：

$u_o\left(t\right)=RC(\frac{{\mathrm{du}}_i\left(t\right)}{dt}-\frac{{\mathrm{du}}_o\left(t\right)}{dt})---\left(6\right)$

为了满足PPG信号处理时的输入输出关系，我们应该使公式(6)离散化。假定脉搏血氧计的采样在时间t是u_i(t)，通过RC电路，相应地输出u_o(t)。假定输入的u_i(t)样本在由时间间隔为ΔT的固定间距点产生一个序列(x₁,x₂,…,x_n)，则与其对应于同一个时刻的输出u_o(t)可由序列(y₁,y₂,…,y_n)表示，由此可以导出以下递推方程：

$y_i=RC(\frac{x_i-x_{i-1}}{\mathrm{\&Delta;}T}-\frac{y_i-y_{i-1}}{\mathrm{\&Delta;}T})---\left(7\right)$

重新排列上式，即可得到公式(1)：

y_i＝αy_i-1+α(x_i-x_i-1) (1)

式中，τ为RC的时间常数，一阶RC高通滤波器的截止频率f_c满足关系f_c＝1/(2πτ)；ΔT为输入序列和输出序列的产生时间间隔，对应于PPG信号的采样周期，i＝1,2,3…n。

本发明实施例中，PPG信号的采样频率为50Hz，则ΔT＝0.02s。为了过滤频率低于0.5Hz的基线漂移和噪声，时间常数可通过τ＝1/(2πf_c)＝1/(2*π*0.5)＝0.3183s计算得到。由于即可求得参数α。

S2、设计低通滤波器，对高通滤波后的信号进行低通滤波。

由于大多数的物理信息位于4Hz和0.5Hz之间的频率，并且比4Hz高的信号往往伴随着噪声，因此要筛除4Hz以上的高频噪声，通常使用截止频率为4Hz的低通滤波器。本发明实施例中，简单起见，我们在这里仅导出一阶RC低通滤波器，如图3所示。

结合基尔霍夫定律和电容器的充电和放电原理，可以建立下面的公式：

$\left\{\begin{array}{c}{u}_i\left(t\right)-u_o\left(t\right)=i\left(t\right)R\\ Q_c\left(t\right)={\mathrm{Cu}}_o\left(t\right)\\ i\left(t\right)=\frac{{\mathrm{dQ}}_c\left(t\right)}{dt}\end{array}\right.---\left(8\right)$

其中，u_o(t)是脉搏血氧仪基于图3所示的一阶RC低通滤波器电路，在时间t，电阻R上流经电流i(t)时的电压输出，输入电压为u_i(t)；Q_c(t)是电容C在时间t的电量；依据电容充放电原理，流经电容C的电流i(t)为电容C上电量的导数dQ_c(t)/dt。

通过替换和重新排列，我们得到：

$u_i\left(t\right)-u_o\left(t\right)=RC\frac{{\mathrm{du}}_o\left(t\right)}{dt}---\left(9\right)$

同理，为了离散化公式(9)，假设输入序列和输出序列产生的时间间隔为ΔT，输入序列和输出序列分别表示为(y₁,y₂,…,y_n)和(z₁,z₂,…,z_n)，则有：

$y_i-z_i=RC\frac{z_i-z_{i-1}}{\mathrm{\&Delta;}T}---\left(10\right)$

重新排列上式，即可得到公式(2)：

z_i＝βy_i+(1-β)z_i-1 (2)

式中τ为RC的时间常数，一阶RC低通滤波器的截止频率f_e满足关系f_e＝1/(2πτ)；ΔT为输入序列和输出序列的产生时间间隔，i＝1,2,3…n。

要过滤4Hz以上的噪声，时间常数的值可以通过τ＝1/(2πf_e)＝1/(2*π*4)＝0.0398s计算。本发明实施例中，PPG信号的采样频率为50Hz，ΔT＝0.02s秒，由于即可求得参数β。

S3、对经过高通滤波和低通滤波后的信号进行平滑滤波处理，使信号平滑化。

经过高通和低通滤波之后，得到的PPG信号大部分为频率为0.5Hz至4Hz的心率等物理信号。然而，在计算出的心率之前，PPG的波形应当进行平滑滤波处理，这样更利于信号周期的计算，以导出心率。平滑后的PPG波形，其峰值的高度可以通过使用最大化函数，或过零检测来确定。

目前有许多平滑滤波器信号处理的方法，如Savitzky-Golay滤波，移动平均和中值滤波等。Savitzky-Golay平滑滤波能有效保持原始信号的形状。然而，该算法复杂度比其他平滑滤波更大，且计算时间长。因此，参考目前的信号处理工具，我们设计了一个简单而快速的平滑滤波器。平滑滤波器的一个重要参数是平滑宽度，平滑宽度越大，噪声降低越大，但信号将被平滑操作扭曲的可能性也越大。平滑宽度的最佳选择取决于信号和数字化区间的宽度和形状。由于PPG信号的波形是一个峰型信号，关键的因素是形成平滑宽度和峰的半宽度的点的数量之间的平滑比率。在一般情况下，增加平滑比率能改善信噪比，但会导致在幅度减少并且在峰的带宽增加。

平滑滤波处理的最终目标是要测量的峰值的位置(峰值的x轴的值)，我们可以使用一个更大的平滑比率(大于0.2)，以获得一个突出的峰型同时平滑比率对峰值位置的影响不大。本发明实施例中，对于50Hz采样率的PPG信号，根据75BMP的心率，峰的半宽度约为20。因此，可以选择10为平滑宽度(通常平滑宽度被设定为一个偶数)。

设输入序列和输出序列分别为(z₁,z₂,…,z_n)和(v₁,v₂,…,v_n)，则有：

v_i+width/2-1＝sum(z_i,z_i+1,...z_width+i-1)/width (3)

其中width为平滑宽度，i＝1,2,…,n-width。

在该平滑算法中，样本中最初和最末的几个点不能被滤波。具体来说，一个m宽度的平滑，会有开始和结尾的(m-1)/2个点不能被计算。通常这些点由0替代，平滑这些点的更好的方法是在开始和结束信号时令平滑宽度逐渐变小，当信号边缘包含关键信息时可以使用该方法。

S4、对平滑滤波处理后的信号求导，以更准确地计算心率。

经过平滑滤波，PPG信号通常有很好的波形，因此可以通过信号更好地判断和估计心率。但出于实际应用考虑，应采取更加可靠的判断算法，在这里，我们使用一阶微分方法。

令x＝(z₁,z₂,…,z_n)，y＝(v₁,v₂,…,v_n)，则信号的一阶导数是随x变化的y的变化率dy/dx，表示在信号上每个点的斜率。假设相邻的点之间的x间隔是固定的，一阶导数表示为：

$\left\{\begin{array}{c}{y}_j^{\&prime;}=\frac{y_{j+1}-y_j}{x_{j+1}-x_j}=\frac{y_{j+1}-y_j}{\mathrm{\&Delta;}x}\\ \begin{array}{cc}{x}_j^{\&prime;}=\frac{x_{j+1}-x_j}{2},& 1<j<n-1\end{array}\end{array}\right.---\left(11\right)$

利用此算法计算三个相邻点的平均斜率可表示为：

$\left\{\begin{array}{c}{y}_j^{\&prime;}=\frac{y_{j+1}-y_{j-1}}{2\mathrm{\&Delta;}x}\\ \begin{array}{cc}{x}_j^{\&prime;}=x_j,& 2<j<n-1\end{array}\end{array}\right.---\left(4\right)$

式中x'_j和y'_j是第j个点的x和y值的导数，Δx为相邻x点的差值。

该算法为中央差分法，因为不涉及导数的x轴位置偏移，所以在此处更适用。

根据公式(4)检测到光电血管容积图信号的峰值，并且通过计算光电血管容积图信号的峰值的间距得到心率。

本发明实施例中，如图4所示，原始PPG信号采用红外光电压信号，运放放大倍数为10⁶。图5为原始PPG信号的幅频特性，即信号的傅里叶变换(FFT)，从图5中可以清晰地看出该信号的频率能量尖峰在1.403Hz，即信号的频率集中在1.403Hz附近，其主要频率包含在[0,5Hz]。

如图6所示，原始PPG信号经步骤S1的高通滤波后，滤除了基线漂移及低于0.5Hz的噪声信号，接着经步骤S2的低通滤波后，滤除了不符合人体生理特征的高于4Hz的高频信号，最后经步骤S3的平滑滤波处理后，PPG信号更利于后期计算。

如图7所示，经过高通、低通、平滑滤波后PPG信号的频率特性，即信号的傅里叶变换(FFT)没有改变，频率峰值均集中在1.403Hz附近，而幅值的改变不影响PPG信号的频率特性。经过三次滤波后，PPG信号的主要频率基本符合我们的预设要求，即集中在[0.5,4Hz]之间。

如图8、图9所示，分别在平稳状态和运动状态下对实际心率值和PPG信号处理后得到的心率值进行比较，可见，无论是平稳状态还是运动状态，本发明所计算得到的心率值均具有较高的准确度。

本领域的普通技术人员将会意识到，这里所述的实施例是为了帮助读者理解本发明的原理，应被理解为本发明的保护范围并不局限于这样的特别陈述和实施例。本领域的普通技术人员可以根据本发明公开的这些技术启示做出各种不脱离本发明实质的其它各种具体变形和组合，这些变形和组合仍然在本发明的保护范围内。

Claims (5)

1.一种基于RC电路的光电血管容积图信号处理方法，其特征在于，包括以下步骤：

S1、设计高通滤波器，对光电血管容积图信号进行高通滤波；

S2、设计低通滤波器，对高通滤波后的信号进行低通滤波；

S3、对经过高通滤波和低通滤波后的信号进行平滑滤波处理，使信号平滑化；

S4、对平滑滤波处理后的信号求导，以更准确地计算心率。

2.根据权利要求1所述的基于RC电路的光电血管容积图信号处理方法，其特征在于，所述步骤S1具体为：

导出一阶RC高通滤波器，设输入序列和输出序列分别为(x₁,x₂,…,x_n)和(y₁,y₂,…,y_n)，则有：

y_i＝αy_i-1+α(x_i-x_i-1) (1)式中，τ为RC的时间常数，一阶RC高通滤波器的截止频率f_c满足关系f_c＝1/(2πτ)；ΔT为输入序列和输出序列的产生时间间隔，i＝1,2,3…n。

3.根据权利要求2所述的基于RC电路的光电血管容积图信号处理方法，其特征在于，所述步骤S2具体为：

导出一阶RC低通滤波器，设输入序列和输出序列分别为(y₁,y₂,…,y_n)和(z₁,z₂,…,z_n)，则有：

z_i＝βy_i+(1-β)z_i-1 (2)

式中，τ为RC的时间常数，一阶RC低通滤波器的截止频率f_e满足关系f_e＝1/(2πτ)；ΔT为输入序列和输出序列的产生时间间隔，i＝1,2,3…n。

4.根据权利要求3所述的基于RC电路的光电血管容积图信号处理方法，其特征在于，所述步骤S3具体为：

设输入序列和输出序列分别为(z₁,z₂,…,z_n)和(v₁,v₂,…,v_n)，则有：

v_i+width/2-1＝sum(z_i,z_i+1,...z_width+i-1)/width (3)

其中width为平滑宽度，i＝1,2,…,n-width。

5.根据权利要求4所述的基于RC电路的光电血管容积图信号处理方法，其特征在于，所述步骤S4具体为：

令x＝(z₁,z₂,…,z_n)，y＝(v₁,v₂,…,v_n)，则通过中央差分法对平滑滤波处理后的信号求导：

$\{\begin{array}{c}{y}_j^{\&prime;}=\frac{y_{j+1}-y_{j-1}}{2\mathrm{\&Delta;}x}\\ \begin{array}{cc}{x}_j^{\&prime;}=x_j,& 2<j<n-1\end{array}\end{array}---\left(4\right)$

式中x'_j和y'_j是第j个点的x和y值的导数，Δx为相邻x点的差值；根据公式(4)检测到光电血管容积图信号的峰值，并且通过计算光电血管容积图信号的峰值的间距得到心率。