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

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

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

方法

算法要求

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

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

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

遭遇信號中斷、噪聲過大或運動偽影時，能夠迅速恢復運行。

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

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

算法開發(fā)

從DFT開始，目標是找到周期，因此DFT方程中的頻率項被替換為周期，并且不是像DFT那樣逐步增加頻率，而是逐步增加周期。DFT以線性方式增加頻率，例如(1f0, 2f0, 3f0, …)，其中f0是第一諧波，而DPT則以采樣周期T0的倍數(shù)為單位，線性增大周期。盡管兩種算法的方程相似，但DFT無法產(chǎn)生與DPT相同的結(jié)果，因為兩種算法有本質(zhì)區(qū)別。通過分析描述其實現(xiàn)的方程，我們可以比較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基函數(shù)的縱坐標值與其之前第N個縱坐標值相同。發(fā)生這種情況的原因是，DFT中的所有諧波之間存在倍數(shù)關(guān)系，高次諧波是低次諧波的整數(shù)倍。

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

DPT中的項N必須針對每個周期進行修改，因為周期之間并非簡單的倍數(shù)關(guān)系，而是相差一個采樣周期，如圖2所示。

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

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

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

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

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

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

算法實現(xiàn)

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

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

MATLAB概念驗證模型

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

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

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

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

結(jié)果

滑動窗口DPT在脈搏血氧測定中的應用

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

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

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

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

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

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

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

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

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

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

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

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

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

圖10. Masimo血氧儀與使用DPT算法的ADI血氧儀的SpO2百分比差異。滿足Bland-Altman標準。

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

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

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

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

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

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

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

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

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

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

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

討論

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

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

結(jié)論

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