生理來(lái)源的信號(hào)可能受到噪聲和運(yùn)動(dòng)偽影的干擾，這些干擾的通帶可能與信號(hào)本身相重疊。生物信號(hào)具有準(zhǔn)平穩(wěn)性，其周期和幅度會(huì)隨時(shí)間而變化。對(duì)于此類信號(hào)，無(wú)法通過(guò)簡(jiǎn)單的數(shù)據(jù)濾波進(jìn)行處理。為了提取信息，一種常用方法是使用與目標(biāo)信號(hào)時(shí)間同步的信號(hào)，作為時(shí)間參考，來(lái)進(jìn)行系綜平均。依靠ECG源的外部心臟觸發(fā)信號(hào)，系綜平均方法可以有效地處理血氧信號(hào)，但在許多情況下，可能無(wú)法獲取ECG源。本文在沒(méi)有ECG觸發(fā)信號(hào)的前提下成功地處理了血氧信號(hào)，并獲得了類似的結(jié)果。

最初，我們開(kāi)發(fā)了一種算法來(lái)執(zhí)行某種形式的自相關(guān)和系綜平均操作。然而，我們很快發(fā)現(xiàn)，時(shí)域中的系綜平均并無(wú)必要，因?yàn)樗邢嚓P(guān)的信息都可以在周期域數(shù)據(jù)本身中找到。心率和血氧飽和度可以直接根據(jù)滑動(dòng)離散周期變換(DPT)產(chǎn)生的結(jié)果計(jì)算出來(lái)。

這項(xiàng)工作始于對(duì)離散傅里葉變換(DFT)的回顧，因?yàn)镈FT能夠生成信號(hào)的頻譜，然后可以利用頻譜確定其周期。該研究的另一個(gè)目標(biāo)是以非常高的分辨率進(jìn)行數(shù)據(jù)采樣。為了利用DFT實(shí)現(xiàn)高分辨率，需要收集大量數(shù)據(jù)樣本。由于生物信號(hào)具有準(zhǔn)平穩(wěn)性，使用DFT收集大量樣本常常會(huì)導(dǎo)致頻譜模糊。我們需要一種分辨率高，且所需樣本量少于DFT的算法。我們希望將該算法用于長(zhǎng)度不確定的實(shí)時(shí)數(shù)據(jù)，因此采用了類似于滑動(dòng)DFT的滑動(dòng)變換形式。

方法

算法要求

我們最初的目標(biāo)是找到一種算法，即使數(shù)據(jù)本質(zhì)上是隨機(jī)且非平穩(wěn)的，也能確定數(shù)據(jù)的潛在基波周期。初始算法要求如下：

能夠確定任何生物醫(yī)學(xué)信號(hào)（如ECG和SpO2）的基波周期。

響應(yīng)時(shí)間足夠快，能夠?qū)崟r(shí)跟蹤心臟心率周期和幅度的變化。

遭遇信號(hào)中斷、噪聲過(guò)大或運(yùn)動(dòng)偽影時(shí)，能夠迅速恢復(fù)運(yùn)行。

計(jì)算速度足夠快，以免成為確定采樣速率的限制因素。

對(duì)存儲(chǔ)空間的要求較低或適中，能夠在低功耗和便攜式設(shè)備中應(yīng)用。

算法開(kāi)發(fā)

從DFT開(kāi)始，目標(biāo)是找到周期，因此DFT方程中的頻率項(xiàng)被替換為周期，并且不是像DFT那樣逐步增加頻率，而是逐步增加周期。DFT以線性方式增加頻率，例如(1f0, 2f0, 3f0, …)，其中f0是第一諧波，而DPT則以采樣周期T0的倍數(shù)為單位，線性增大周期。盡管兩種算法的方程相似，但DFT無(wú)法產(chǎn)生與DPT相同的結(jié)果，因?yàn)閮煞N算法有本質(zhì)區(qū)別。通過(guò)分析描述其實(shí)現(xiàn)的方程，我們可以比較DFT和DPT。對(duì)于采樣頻率 fS，N點(diǎn)DFT的頻點(diǎn)k對(duì)應(yīng)頻率fK = k × fS / N Hz，公式1是樣本序列XI... XI + N - 1的第k個(gè)頻點(diǎn)的頻譜表達(dá)式。

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

DFT的第i個(gè)樣本按照公式2進(jìn)行計(jì)算。

如圖1所示，對(duì)于每個(gè)諧波，DFT基函數(shù)的縱坐標(biāo)值與其之前第N個(gè)縱坐標(biāo)值相同。發(fā)生這種情況的原因是，DFT中的所有諧波之間存在倍數(shù)關(guān)系，高次諧波是低次諧波的整數(shù)倍。

圖1. 傅里葉變換正弦基函數(shù)，紅色所示為第一諧波(1 Hz)，藍(lán)色為第二諧 波(2 Hz)，綠色為第三諧波(3 Hz)。

DPT中的項(xiàng)N必須針對(duì)每個(gè)周期進(jìn)行修改，因?yàn)橹芷谥g并非簡(jiǎn)單的倍數(shù)關(guān)系，而是相差一個(gè)采樣周期，如圖2所示。

圖2. 三個(gè)相鄰正弦函數(shù)和三個(gè)相鄰余弦函數(shù)的周期變換復(fù)正弦基函數(shù)。 此示例假設(shè)這些函數(shù)的采樣時(shí)間間隔為10 ms。

滑動(dòng)形式的DFT和DPT都需要實(shí)現(xiàn)循環(huán)或遞歸緩沖區(qū)，用于保存數(shù)量固定的最新樣本。當(dāng)輸入數(shù)據(jù)為實(shí)數(shù)時(shí)，使用一個(gè)緩沖區(qū)；而當(dāng)輸入數(shù)據(jù)為復(fù)數(shù)時(shí)，則使用兩個(gè)緩沖區(qū)。DPT變換的第i個(gè)樣本可以套入公式3。

其中，RBS為遞歸緩沖區(qū)大小，TL為最長(zhǎng)周期的長(zhǎng)度，TN為當(dāng)前正在處理的基元的周期。這樣做可以使每個(gè)基礎(chǔ)周期的起始和終止縱坐標(biāo)值相同。周期s從最小周期延伸到所選的最大周期，以覆蓋采樣數(shù)據(jù)中的周期。該實(shí)現(xiàn)利用了一組基函數(shù)，這些基函數(shù)代表了圖2中復(fù)正弦波的增量相位角。

DPT的實(shí)現(xiàn)之所以有些困難，是因?yàn)榛瘮?shù)由多組復(fù)函數(shù)組成，這些復(fù)函數(shù)之間大多不是倍數(shù)關(guān)系，而且采樣周期不同。高效的DPT變換需使用圖3所示的基礎(chǔ)相位角。這也是本文所采用的實(shí)現(xiàn)形式。

圖3. 周期變換基礎(chǔ)相位角，展示了復(fù)相位角的值如何隨著每分鐘采樣周期數(shù)的增加而變化。上升曲線表示余弦相位角，下降曲線表示正弦相位角。

使用公式4可以輕松得出相量，其中“s”是以采樣周期為步長(zhǎng)，從最小選定周期到最大選定周期的周期集。

算法實(shí)現(xiàn)

滑動(dòng)DPT變換使用IIR濾波器實(shí)現(xiàn)，其信號(hào)流圖中在一個(gè)梳狀濾波器后接了一個(gè)諧振器，這與滑動(dòng)形式DFT的實(shí)現(xiàn)類似。N個(gè)樣本的梳狀濾波器延遲導(dǎo)致瞬態(tài)響應(yīng)的長(zhǎng)度為N-1個(gè)樣本。已經(jīng)有人使用心率調(diào)諧的梳狀濾波器并取得了一定的成功。DPT復(fù)基函數(shù)或相位角的分量并非總是諧波相關(guān)，因此這些函數(shù)的端點(diǎn)不會(huì)在樣本空間中形成連續(xù)函數(shù)，這與DFT不同。然而，如果將DPT實(shí)現(xiàn)為滑動(dòng)變換，那么基函數(shù)就會(huì)被“包裹”起來(lái)，從而使基函數(shù)的分量變成連續(xù)的。當(dāng)數(shù)據(jù)和基函數(shù)滑動(dòng)時(shí)，計(jì)算它們的相關(guān)性，基函數(shù)連續(xù)性得以保持。

在滑動(dòng)窗口算法中，長(zhǎng)度為N的窗口在長(zhǎng)度不確定的數(shù)據(jù)數(shù)組上滑動(dòng)。對(duì)于DPT而言，由于DPT可以處理實(shí)部和虛部?jī)深愝斎霐?shù)據(jù)，因此需要維護(hù)兩個(gè)遞歸緩沖區(qū)。如果輸入只有一個(gè)實(shí)部（通常情況如此），則只需使用一個(gè)遞歸緩沖區(qū)。然而，根據(jù)輸入和基函數(shù)之間的相位關(guān)系，結(jié)果仍然可能是復(fù)數(shù)。結(jié)果存儲(chǔ)在兩個(gè)系綜緩沖區(qū)中，每個(gè)緩沖區(qū)的長(zhǎng)度為所選的最大周期。

MATLAB概念驗(yàn)證模型

我們通過(guò)MATLAB 腳本實(shí)現(xiàn)了公式4。圖4使用正弦和余弦函數(shù)作為輸入，幅度為±1，周期為45 ms、79 ms和175 ms。MATLAB腳本的周期限定在400 ms（200個(gè)周期/分鐘）到2 s（40個(gè)周期/分鐘）之間。本例總共處理了5000個(gè)數(shù)據(jù)樣本，樣本數(shù)量固定不變。由于輸入數(shù)據(jù)是幅度為1的正弦波形，因此每個(gè)周期的幅度也為1。很容易看出，這種變換實(shí)現(xiàn)的分辨率非常高。

圖4. 幅度譜，展示了彼此不成倍數(shù)關(guān)系的三組輸入正弦數(shù)據(jù)的值。

圖5為每分鐘73個(gè)周期、幅度為4.5的正弦余弦波的結(jié)果。此示例使用了長(zhǎng)度為1500個(gè)數(shù)據(jù)點(diǎn)的遞歸緩沖區(qū)。請(qǐng)注意，存在一些較小誤差，幅度誤差為0.366%，周期誤差為0.234%。對(duì)于生物醫(yī)學(xué)應(yīng)用而言，這些誤差的大小一般是可以接受的。在外周毛細(xì)血管血氧飽和度(SpO2)測(cè)量中，這些誤差無(wú)關(guān)緊要，因?yàn)镾pO2是根據(jù)紅光和紅外光譜信號(hào)的比率之比來(lái)計(jì)算的。參見(jiàn)公式6和公式7。

圖5. 余弦波形的滑動(dòng)周期變換，每分鐘73個(gè)周期，振幅為4.5。幅度誤差小于0.37%，周期誤差小于0.24%。

結(jié)果

滑動(dòng)窗口DPT在脈搏血氧測(cè)定中的應(yīng)用

將滑動(dòng)窗口算法應(yīng)用于脈搏血氧測(cè)定時(shí)，為使算法正常運(yùn)行，需要兩個(gè)遞歸數(shù)組：一個(gè)用于存儲(chǔ)紅光歷史記錄，另一個(gè)用于存儲(chǔ)紅外歷史記錄。為完成滑動(dòng)變換，還需根據(jù)相應(yīng)周期的基函數(shù)，旋轉(zhuǎn)遞歸緩沖區(qū)（其長(zhǎng)度與正在處理的周期點(diǎn)相同）中更新的內(nèi)容。該緩沖區(qū)的長(zhǎng)度決定了整體分辨率，一旦有足夠多的數(shù)據(jù)進(jìn)入處理流程以填充這些緩沖區(qū)，變換結(jié)果就會(huì)達(dá)到一個(gè)穩(wěn)定的極限，只有幅度或周期會(huì)隨著輸入數(shù)據(jù)的變化而改變。對(duì)于所報(bào)告的數(shù)據(jù)處理，遞歸緩沖區(qū)保存最后10秒的數(shù)據(jù)。

原始數(shù)據(jù)由ADI公司的研究人員收集，用于處理數(shù)據(jù)的軟件是MATLAB腳本中的滑動(dòng)DPT。圖6為從某位受試者獲取的原始數(shù)據(jù)；經(jīng)過(guò)1 Hz至4 Hz帶通濾波的數(shù)據(jù)，以及利用總寬度為200 ms的平坦光滑移動(dòng)平均濾波器處理后的數(shù)據(jù)。圖7為填充遞歸緩沖區(qū)之后頻譜達(dá)到穩(wěn)定幅度的頻譜。隨著新數(shù)據(jù)被采樣，DPT將持續(xù)跟蹤原始數(shù)據(jù)中的所有變化，頻譜也會(huì)相應(yīng)地更新。

圖6. 使用MAX30101 PPG AFE器件從某位受試者獲取的原始光電容積脈搏波數(shù)據(jù)、經(jīng)濾波的數(shù)據(jù)和經(jīng)平滑處理后的數(shù)據(jù)。上方波形表示原始紅外和紅光信號(hào)，而下方波形表示經(jīng)過(guò)濾波和平滑處理的數(shù)據(jù)。

圖7. 此圖為采用滑動(dòng)窗口DPT處理的紅光和紅外光譜。兩個(gè)波峰中較大的是紅外光譜；較小的是紅光光譜。

為了估算SpO2，先需使用比率之比公式。交流分量使用圖7所示頻譜的峰值，直流分量使用圖6所示未濾波信號(hào)的平均值。

比較從Masimo血氧儀使用SET算法收集的SpO2和心率數(shù)據(jù)，與使用ADI MAX30101脈搏血氧儀傳感器同時(shí)獲取的數(shù)據(jù)。隨機(jī)選擇某位受試者的數(shù)據(jù)，并將結(jié)果繪制在圖8和圖9中。

圖8. DPT處理的光電容積脈搏波數(shù)據(jù)比較。

圖9. 比較來(lái)自MAX30101血氧儀（采用離散周期變換進(jìn)行處理）和Masimo血 氧儀的心率數(shù)據(jù)。

評(píng)估兩種不同儀器測(cè)量同一參數(shù)所產(chǎn)生的數(shù)值，是常見(jiàn)的醫(yī)學(xué)做法。其中一種儀器被認(rèn)為能夠產(chǎn)生正確的結(jié)果，用作標(biāo)準(zhǔn)儀器。

Bland和Altman開(kāi)發(fā)了一種用于評(píng)估兩種定量測(cè)量結(jié)果一致性的方法。他們通過(guò)分析平均差異和構(gòu)建一致性界限來(lái)判斷一致性。Bland-Altman圖分析是評(píng)估平均差異之間的偏差和估計(jì)一致性區(qū)間的一種簡(jiǎn)單方法。如果對(duì)兩臺(tái)醫(yī)療儀器開(kāi)展此項(xiàng)測(cè)試，其中一臺(tái)被視為標(biāo)準(zhǔn)，則另一臺(tái)儀器的結(jié)果必須在標(biāo)準(zhǔn)儀器結(jié)果的兩個(gè)標(biāo)準(zhǔn)差或95%范圍內(nèi)，才能認(rèn)為其在臨床應(yīng)用上與標(biāo)準(zhǔn)儀器效果相當(dāng)。

與相關(guān)分析研究?jī)蓚€(gè)變量之間的關(guān)系不同，Bland-Altman方法是一種統(tǒng)計(jì)學(xué)方法，關(guān)注的是兩個(gè)變量之間的差異。

我們利用MAX30101脈搏血氧儀傳感器收集了26名健康成年受試者的數(shù)據(jù)，并將其與Masimo血氧儀（其中融合了新型信號(hào)提取技術(shù)Signal Extraction Technology）的測(cè)量結(jié)果進(jìn)行比較，從而評(píng)估DPT算法的準(zhǔn)確性和精確度。研究對(duì)象包括15名男性和11名女性受試者，年齡在20至40歲之間。這項(xiàng)研究的目的是比較同一受試者使用兩種血氧儀的測(cè)量結(jié)果，而不是男性和女性之間的差異。請(qǐng)注意，兩性之間的SpO2確實(shí)略有不同。一項(xiàng)研究表明，對(duì)于年輕健康成年人，男性的平均SpO2為97.1±1.2%，而女性的平均SpO2為98.6±1.0%。

圖10和圖11位使用Bland-Altman標(biāo)準(zhǔn)的結(jié)果，每個(gè)圓圈代表一位受試者的Bland-Altman結(jié)果。所有SpO2比較均符合Bland-Altman標(biāo)準(zhǔn)。

圖10. Masimo血氧儀與使用DPT算法的ADI血氧儀的SpO2百分比差異。滿足Bland-Altman標(biāo)準(zhǔn)。

圖11. Masimo血氧儀與使用DPT算法的ADI血氧儀的每分鐘心率差異。除一個(gè)案例外，其他所有案例均滿足Bland-Altman標(biāo)準(zhǔn)。箭頭標(biāo)示了超出兩個(gè)標(biāo)準(zhǔn)差范圍的分析結(jié)果。

在圖11中，箭頭指向的心率比較值超出了兩個(gè)標(biāo)準(zhǔn)差范圍。該受試者的心率與時(shí)間關(guān)系圖如圖12所示，其中Masimo血氧儀的標(biāo)準(zhǔn)差為1.7892，而使用DPT算法的MAX30101血氧儀的標(biāo)準(zhǔn)差為0.8935。在這種情況下，我們很難確定哪種儀器更準(zhǔn)確，但可以從標(biāo)準(zhǔn)差中找到一些線索。

圖12. Masimo血氧儀和ADI血氧儀的心率與時(shí)間的關(guān)系圖。在25秒周期內(nèi)， Masimo血氧儀的標(biāo)準(zhǔn)差為1.7892，而MAX30101的標(biāo)準(zhǔn)差為0.8935。階梯波形是 來(lái)自Masimo血氧儀的信號(hào)；平滑信號(hào)來(lái)自運(yùn)行DPT算法的ADI血氧儀。

采用SDPT算法的血氧儀系統(tǒng)原型

最后，我們采用Arm微處理器（運(yùn)行裸機(jī)操作系統(tǒng)），設(shè)計(jì)了一個(gè)血氧儀原型。使用樹莓派Zero作為計(jì)算機(jī)平臺(tái)，MAX30102集成電路用作傳感器。操作系統(tǒng)和滑動(dòng)窗口DPT采用標(biāo)準(zhǔn)C語(yǔ)言實(shí)現(xiàn)。圖13即為該原型。整個(gè)血氧儀由USB 3.0連接供電。兩個(gè)數(shù)模轉(zhuǎn)換器根據(jù)監(jiān)控軟件的判斷，通過(guò)帶狀電纜將數(shù)據(jù)發(fā)送到Tektronix DPO-4034示波器，在其中繪制圖像。然后，圖像通過(guò)網(wǎng)絡(luò)連接發(fā)送到臺(tái)式計(jì)算機(jī)。圖15為大約9秒的時(shí)間內(nèi)從單個(gè)受試者獲得的結(jié)果，之后用10秒的時(shí)間來(lái)填充遞歸緩沖區(qū)。

圖13. 基于樹莓派Zero的脈搏血氧儀原型。MAX30102 SpO2傳感器位于圖片左上角所示的指夾中。

通過(guò)一階低通IIR濾波器從原始信號(hào)中提取紅光和紅外直流信號(hào)；通過(guò)一階高通IIR濾波器提取交流信號(hào)。參見(jiàn)圖14。這些濾波器的時(shí)間常數(shù)設(shè)置為大約1秒。數(shù)據(jù)以100 SPS的速率采樣，并以MAX30102的中斷作為定時(shí)信號(hào)。對(duì)于紅光和紅外信號(hào)，該器件的輸出均為12位定點(diǎn)數(shù)字格式。

圖14. 使用無(wú)限脈沖響應(yīng)(IIR)濾波器從原始光譜數(shù)據(jù)中提取交流信號(hào)和直流信號(hào)。

圖15. 樹莓派Zero血氧儀原型產(chǎn)生的PPG波形，上方為紅光脈沖，下方為紅外脈沖。心率約為58 bpm。圖中所示為倒置的波形，以便更準(zhǔn)確地表示手指中的實(shí)際動(dòng)脈壓力。

紅光和紅外交流信號(hào)通過(guò)濾波器提取出來(lái)之后，就交由DPT處理，而無(wú)需任何進(jìn)一步的信號(hào)預(yù)處理。光譜信號(hào)的第一諧波產(chǎn)生的峰值如圖16所示。心率由橫坐標(biāo)上數(shù)據(jù)峰值的位置決定，而SpO2通過(guò)比率之比公式使用紅光和紅外數(shù)據(jù)峰值的幅度直接計(jì)算。

圖16.樹莓派Zero使用滑動(dòng)窗口離散周期變換生成的頻譜，SpO2值為97%，心率為58 bpm。光標(biāo)b（中心垂直藍(lán)線）顯示測(cè)量的心跳周期為1.03秒。左上角的矩形信號(hào)指向橫坐標(biāo)上400 ms周期的位置；右上角的矩形信號(hào)指向橫坐標(biāo)上2000 ms周期的位置。

討論

血氧儀產(chǎn)生的原始光信號(hào)包含較大的穩(wěn)定直流分量和較小的振蕩交流分量，后者約為直流信號(hào)的1%。這些振蕩分量反映的是毛細(xì)血管中的脈動(dòng)活動(dòng)。任何運(yùn)動(dòng)或其他偽影都可能輕易覆蓋這些信號(hào)，使讀數(shù)不準(zhǔn)確。多年來(lái)，人們花費(fèi)了大量時(shí)間來(lái)研究將這些信號(hào)與偽影分離的方法。事實(shí)證明，這些方法通常非常復(fù)雜且難以實(shí)施。

正是出于這些原因，我們才開(kāi)展了這項(xiàng)研究。DPT算法采用的變換只需少量樣本，但卻能實(shí)現(xiàn)準(zhǔn)確的測(cè)量，許多挑戰(zhàn)因此迎刃而解。在周期域內(nèi)進(jìn)行測(cè)量，并按采樣周期將每個(gè)周期點(diǎn)分隔開(kāi)來(lái)，便能提供所需的分辨率。然后，我們可以利用來(lái)自DPT的周期和幅度信息直接計(jì)算心率和血氧飽和度，而無(wú)需返回時(shí)域。

結(jié)論

采用增量DPT算法的周期域分析，是處理周期性生物醫(yī)學(xué)信號(hào)以獲得頻譜成分的有效方法。該方法支持頻域分析，而且在實(shí)現(xiàn)上也有優(yōu)勢(shì)。研究表明，運(yùn)行DPT算法的ADI MAX30101集成電路傳感器足夠精確，在醫(yī)療實(shí)踐中能夠取代Masimo血氧儀。