生理來源的信號可能受到噪聲和運動偽影的干擾，這些干擾的通帶可能與信號本身相重疊。生物信號具有準平穩性，其周期和幅度會隨時間而變化。對于此類信號，無法通過簡單的數據濾波進行處理。為了提取信息，一種常用方法是使用與目標信號時間同步的信號，作為時間參考，來進行系綜平均。依靠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血氧儀。