滑动离散周期变换算法在脉搏血氧测定中的应用

生理来源的信号可能受到噪声和运动伪影的干扰，这些干扰的通带可能与信号本身相重叠。生物信号具有准平稳性，其周期和幅度会随时间而变化。对于此类信号，无法通过简单的数据滤波进行处理。为了提取信息，一种常用方法是使用与目标信号时间同步的信号，作为时间参考，来进行系综平均。依靠ECG源的外部心脏触发信号，系综平均方法可以有效地处理血氧信号，但在许多情况下，可能无法获取ECG源。本文在没有ECG触发信号的前提下成功地处理了血氧信号，并获得了类似的结果。

最初，我们开发了一种算法来执行某种形式的自相关和系综平均操作。然而，我们很快发现，时域中的系综平均并无必要，因为所有相关的信息都可以在周期域数据本身中找到。心率和血氧饱和度可以直接根据滑动离散周期变换(DPT)产生的结果计算出来。

这项工作始于对离散傅里叶变换(DFT)的回顾，因为DFT能够生成信号的频谱，然后可以利用频谱确定其周期。该研究的另一个目标是以非常高的分辨率进行数据采样。为了利用DFT实现高分辨率，需要收集大量数据样本。由于生物信号具有准平稳性，使用DFT收集大量样本常常会导致频谱模糊。我们需要一种分辨率高，且所需样本量少于DFT的算法。我们希望将该算法用于长度不确定的实时数据，因此采用了类似于滑动DFT的滑动变换形式。

方法

算法要求

我们最初的目标是找到一种算法，即使数据本质上是随机且非平稳的，也能确定数据的潜在基波周期。初始算法要求如下：

能够确定任何生物医学信号（如ECG和SpO2）的基波周期。

响应时间足够快，能够实时跟踪心脏心率周期和幅度的变化。

遭遇信号中断、噪声过大或运动伪影时，能够迅速恢复运行。

计算速度足够快，以免成为确定采样速率的限制因素。

对存储空间的要求较低或适中，能够在低功耗和便携式设备中应用。

算法开发

从DFT开始，目标是找到周期，因此DFT方程中的频率项被替换为周期，并且不是像DFT那样逐步增加频率，而是逐步增加周期。DFT以线性方式增加频率，例如(1f0, 2f0, 3f0, …)，其中f0是第一谐波，而DPT则以采样周期T0的倍数为单位，线性增大周期。尽管两种算法的方程相似，但DFT无法产生与DPT相同的结果，因为两种算法有本质区别。通过分析描述其实现的方程，我们可以比较DFT和DPT。对于采样频率 fS，N点DFT的频点k对应频率fK = k × fS / N Hz，公式1是样本序列XI... XI + N - 1的第k个频点的频谱表达式。

其中， k = 0, 1, 2, ... N – 1

DFT的第i个样本按照公式2进行计算。

如图1所示，对于每个谐波，DFT基函数的纵坐标值与其之前第N个纵坐标值相同。发生这种情况的原因是，DFT中的所有谐波之间存在倍数关系，高次谐波是低次谐波的整数倍。

图1. 傅里叶变换正弦基函数，红色所示为第一谐波(1 Hz)，蓝色为第二谐 波(2 Hz)，绿色为第三谐波(3 Hz)。

DPT中的项N必须针对每个周期进行修改，因为周期之间并非简单的倍数关系，而是相差一个采样周期，如图2所示。

图2. 三个相邻正弦函数和三个相邻余弦函数的周期变换复正弦基函数。 此示例假设这些函数的采样时间间隔为10 ms。

滑动形式的DFT和DPT都需要实现循环或递归缓冲区，用于保存数量固定的最新样本。当输入数据为实数时，使用一个缓冲区；而当输入数据为复数时，则使用两个缓冲区。DPT变换的第i个样本可以套入公式3。

其中，RBS为递归缓冲区大小，TL为最长周期的长度，TN为当前正在处理的基元的周期。这样做可以使每个基础周期的起始和终止纵坐标值相同。周期s从最小周期延伸到所选的最大周期，以覆盖采样数据中的周期。该实现利用了一组基函数，这些基函数代表了图2中复正弦波的增量相位角。

DPT的实现之所以有些困难，是因为基函数由多组复函数组成，这些复函数之间大多不是倍数关系，而且采样周期不同。高效的DPT变换需使用图3所示的基础相位角。这也是本文所采用的实现形式。

图3. 周期变换基础相位角，展示了复相位角的值如何随着每分钟采样周期数的增加而变化。上升曲线表示余弦相位角，下降曲线表示正弦相位角。

使用公式4可以轻松得出相量，其中“s”是以采样周期为步长，从最小选定周期到最大选定周期的周期集。

算法实现

滑动DPT变换使用IIR滤波器实现，其信号流图中在一个梳状滤波器后接了一个谐振器，这与滑动形式DFT的实现类似。N个样本的梳状滤波器延迟导致瞬态响应的长度为N-1个样本。已经有人使用心率调谐的梳状滤波器并取得了一定的成功。DPT复基函数或相位角的分量并非总是谐波相关，因此这些函数的端点不会在样本空间中形成连续函数，这与DFT不同。然而，如果将DPT实现为滑动变换，那么基函数就会被“包裹”起来，从而使基函数的分量变成连续的。当数据和基函数滑动时，计算它们的相关性，基函数连续性得以保持。

在滑动窗口算法中，长度为N的窗口在长度不确定的数据数组上滑动。对于DPT而言，由于DPT可以处理实部和虚部两类输入数据，因此需要维护两个递归缓冲区。如果输入只有一个实部（通常情况如此），则只需使用一个递归缓冲区。然而，根据输入和基函数之间的相位关系，结果仍然可能是复数。结果存储在两个系综缓冲区中，每个缓冲区的长度为所选的最大周期。

MATLAB概念验证模型

我们通过MATLAB 脚本实现了公式4。图4使用正弦和余弦函数作为输入，幅度为±1，周期为45 ms、79 ms和175 ms。MATLAB脚本的周期限定在400 ms（200个周期/分钟）到2 s（40个周期/分钟）之间。本例总共处理了5000个数据样本，样本数量固定不变。由于输入数据是幅度为1的正弦波形，因此每个周期的幅度也为1。很容易看出，这种变换实现的分辨率非常高。

图4. 幅度谱，展示了彼此不成倍数关系的三组输入正弦数据的值。

图5为每分钟73个周期、幅度为4.5的正弦余弦波的结果。此示例使用了长度为1500个数据点的递归缓冲区。请注意，存在一些较小误差，幅度误差为0.366%，周期误差为0.234%。对于生物医学应用而言，这些误差的大小一般是可以接受的。在外周毛细血管血氧饱和度(SpO2)测量中，这些误差无关紧要，因为SpO2是根据红光和红外光谱信号的比率之比来计算的。参见公式6和公式7。

图5. 余弦波形的滑动周期变换，每分钟73个周期，振幅为4.5。幅度误差小于0.37%，周期误差小于0.24%。

结果

滑动窗口DPT在脉搏血氧测定中的应用

将滑动窗口算法应用于脉搏血氧测定时，为使算法正常运行，需要两个递归数组：一个用于存储红光历史记录，另一个用于存储红外历史记录。为完成滑动变换，还需根据相应周期的基函数，旋转递归缓冲区（其长度与正在处理的周期点相同）中更新的内容。该缓冲区的长度决定了整体分辨率，一旦有足够多的数据进入处理流程以填充这些缓冲区，变换结果就会达到一个稳定的极限，只有幅度或周期会随着输入数据的变化而改变。对于所报告的数据处理，递归缓冲区保存最后10秒的数据。

原始数据由ADI公司的研究人员收集，用于处理数据的软件是MATLAB脚本中的滑动DPT。图6为从某位受试者获取的原始数据；经过1 Hz至4 Hz带通滤波的数据，以及利用总宽度为200 ms的平坦光滑移动平均滤波器处理后的数据。图7为填充递归缓冲区之后频谱达到稳定幅度的频谱。随着新数据被采样，DPT将持续跟踪原始数据中的所有变化，频谱也会相应地更新。

图6. 使用MAX30101 PPG AFE器件从某位受试者获取的原始光电容积脉搏波数据、经滤波的数据和经平滑处理后的数据。上方波形表示原始红外和红光信号，而下方波形表示经过滤波和平滑处理的数据。

图7. 此图为采用滑动窗口DPT处理的红光和红外光谱。两个波峰中较大的是红外光谱；较小的是红光光谱。

为了估算SpO2，先需使用比率之比公式。交流分量使用图7所示频谱的峰值，直流分量使用图6所示未滤波信号的平均值。

比较从Masimo血氧仪使用SET算法收集的SpO2和心率数据，与使用ADI MAX30101脉搏血氧仪传感器同时获取的数据。随机选择某位受试者的数据，并将结果绘制在图8和图9中。

图8. DPT处理的光电容积脉搏波数据比较。

图9. 比较来自MAX30101血氧仪（采用离散周期变换进行处理）和Masimo血 氧仪的心率数据。

评估两种不同仪器测量同一参数所产生的数值，是常见的医学做法。其中一种仪器被认为能够产生正确的结果，用作标准仪器。

Bland和Altman开发了一种用于评估两种定量测量结果一致性的方法。他们通过分析平均差异和构建一致性界限来判断一致性。Bland-Altman图分析是评估平均差异之间的偏差和估计一致性区间的一种简单方法。如果对两台医疗仪器开展此项测试，其中一台被视为标准，则另一台仪器的结果必须在标准仪器结果的两个标准差或95%范围内，才能认为其在临床应用上与标准仪器效果相当。

与相关分析研究两个变量之间的关系不同，Bland-Altman方法是一种统计学方法，关注的是两个变量之间的差异。

我们利用MAX30101脉搏血氧仪传感器收集了26名健康成年受试者的数据，并将其与Masimo血氧仪（其中融合了新型信号提取技术Signal Extraction Technology）的测量结果进行比较，从而评估DPT算法的准确性和精确度。研究对象包括15名男性和11名女性受试者，年龄在20至40岁之间。这项研究的目的是比较同一受试者使用两种血氧仪的测量结果，而不是男性和女性之间的差异。请注意，两性之间的SpO2确实略有不同。一项研究表明，对于年轻健康成年人，男性的平均SpO2为97.1±1.2%，而女性的平均SpO2为98.6±1.0%。

图10和图11位使用Bland-Altman标准的结果，每个圆圈代表一位受试者的Bland-Altman结果。所有SpO2比较均符合Bland-Altman标准。

图10. Masimo血氧仪与使用DPT算法的ADI血氧仪的SpO2百分比差异。满足Bland-Altman标准。

图11. Masimo血氧仪与使用DPT算法的ADI血氧仪的每分钟心率差异。除一个案例外，其他所有案例均满足Bland-Altman标准。箭头标示了超出两个标准差范围的分析结果。

在图11中，箭头指向的心率比较值超出了两个标准差范围。该受试者的心率与时间关系图如图12所示，其中Masimo血氧仪的标准差为1.7892，而使用DPT算法的MAX30101血氧仪的标准差为0.8935。在这种情况下，我们很难确定哪种仪器更准确，但可以从标准差中找到一些线索。

图12. Masimo血氧仪和ADI血氧仪的心率与时间的关系图。在25秒周期内， Masimo血氧仪的标准差为1.7892，而MAX30101的标准差为0.8935。阶梯波形是 来自Masimo血氧仪的信号；平滑信号来自运行DPT算法的ADI血氧仪。

采用SDPT算法的血氧仪系统原型

最后，我们采用Arm微处理器（运行裸机操作系统），设计了一个血氧仪原型。使用树莓派Zero作为计算机平台，MAX30102集成电路用作传感器。操作系统和滑动窗口DPT采用标准C语言实现。图13即为该原型。整个血氧仪由USB 3.0连接供电。两个数模转换器根据监控软件的判断，通过带状电缆将数据发送到Tektronix DPO-4034示波器，在其中绘制图像。然后，图像通过网络连接发送到台式计算机。图15为大约9秒的时间内从单个受试者获得的结果，之后用10秒的时间来填充递归缓冲区。

图13. 基于树莓派Zero的脉搏血氧仪原型。MAX30102 SpO2传感器位于图片左上角所示的指夹中。

通过一阶低通IIR滤波器从原始信号中提取红光和红外直流信号；通过一阶高通IIR滤波器提取交流信号。参见图14。这些滤波器的时间常数设置为大约1秒。数据以100 SPS的速率采样，并以MAX30102的中断作为定时信号。对于红光和红外信号，该器件的输出均为12位定点数字格式。

图14. 使用无限脉冲响应(IIR)滤波器从原始光谱数据中提取交流信号和直流信号。

图15. 树莓派Zero血氧仪原型产生的PPG波形，上方为红光脉冲，下方为红外脉冲。心率约为58 bpm。图中所示为倒置的波形，以便更准确地表示手指中的实际动脉压力。

红光和红外交流信号通过滤波器提取出来之后，就交由DPT处理，而无需任何进一步的信号预处理。光谱信号的第一谐波产生的峰值如图16所示。心率由横坐标上数据峰值的位置决定，而SpO2通过比率之比公式使用红光和红外数据峰值的幅度直接计算。

图16.树莓派Zero使用滑动窗口离散周期变换生成的频谱，SpO2值为97%，心率为58 bpm。光标b（中心垂直蓝线）显示测量的心跳周期为1.03秒。左上角的矩形信号指向横坐标上400 ms周期的位置；右上角的矩形信号指向横坐标上2000 ms周期的位置。

讨论

血氧仪产生的原始光信号包含较大的稳定直流分量和较小的振荡交流分量，后者约为直流信号的1%。这些振荡分量反映的是毛细血管中的脉动活动。任何运动或其他伪影都可能轻易覆盖这些信号，使读数不准确。多年来，人们花费了大量时间来研究将这些信号与伪影分离的方法。事实证明，这些方法通常非常复杂且难以实施。

正是出于这些原因，我们才开展了这项研究。DPT算法采用的变换只需少量样本，但却能实现准确的测量，许多挑战因此迎刃而解。在周期域内进行测量，并按采样周期将每个周期点分隔开来，便能提供所需的分辨率。然后，我们可以利用来自DPT的周期和幅度信息直接计算心率和血氧饱和度，而无需返回时域。

结论

采用增量DPT算法的周期域分析，是处理周期性生物医学信号以获得频谱成分的有效方法。该方法支持频域分析，而且在实现上也有优势。研究表明，运行DPT算法的ADI MAX30101集成电路传感器足够精确，在医疗实践中能够取代Masimo血氧仪。