【導讀】生理信號固有的準周期性特征表現為時(shí)變幅頻特性，其有效信息常與運動(dòng)偽影及環(huán)境噪聲在頻域高度重疊。傳統頻域濾波手段因信號非平穩性面臨失效風(fēng)險，業(yè)界普遍采用時(shí)間鎖相式系綜平均法——通過(guò)ECG心電觸發(fā)信號建立時(shí)域參考基準，實(shí)現血氧參數的有效提取。然而臨床場(chǎng)景中ECG信號的可及性限制催生了技術(shù)盲區。本研究突破性開(kāi)發(fā)自主周期定位算法，在不依賴(lài)外部心電參考的條件下，達成與ECG引導方案相當的信號重構精度，為可穿戴醫療設備提供了去ECG依賴(lài)的生理參數解析新范式。
摘要
本文介紹了新型滑動(dòng)離散周期變換(DPT)算法，可設計用于處理生理信號，尤其是脈搏血氧儀采集的光電容積脈搏波(PPG)信號。該算法采用正弦基函數進(jìn)行周期域分析，可解決隨機噪聲和非平穩數據等難題。DPT在MATLAB?中作為滑動(dòng)變換實(shí)現，結合了自相關(guān)與系綜平均。文中將詳細介紹在A(yíng)DI MAX30101器件上開(kāi)發(fā)和實(shí)現的一種算法，并與采用Signal Extraction Technology? (SET)的Masimo血氧儀進(jìn)行比較。
簡(jiǎn)介
生理信號固有的準周期性特征表現為時(shí)變幅頻特性，其有效信息常與運動(dòng)偽影及環(huán)境噪聲在頻域高度重疊。傳統頻域濾波手段因信號非平穩性面臨失效風(fēng)險，業(yè)界普遍采用時(shí)間鎖相式系綜平均法——通過(guò)ECG心電觸發(fā)信號建立時(shí)域參考基準，實(shí)現血氧參數的有效提取。然而臨床場(chǎng)景中ECG信號的可及性限制催生了技術(shù)盲區。本研究突破性開(kāi)發(fā)自主周期定位算法，在不依賴(lài)外部心電參考的條件下，達成與ECG引導方案相當的信號重構精度，為可穿戴醫療設備提供了去ECG依賴(lài)的生理參數解析新范式。
最初，我們開(kāi)發(fā)了一種算法來(lái)執行某種形式的自相關(guān)和系綜平均處理4。然而，我們很快發(fā)現，時(shí)域中的系綜平均并無(wú)必要，因為所有相關(guān)的信息都可以在周期域數據本身中找到。心率和血氧飽和度可以直接根據滑動(dòng)離散周期變換(DPT)產(chǎn)生的結果計算出來(lái)。
這項工作始于對離散傅里葉變換(DFT)的回顧，因為DFT能夠生成信號的頻譜，然后可以利用頻譜確定其周期5,6。該研究的另一個(gè)目標是以非常高的分辨率進(jìn)行數據采樣。為了利用DFT實(shí)現高分辨率，需要收集大量數據樣本。由于生物信號具有準平穩性，使用DFT收集大量樣本常常會(huì )導致頻譜模糊7。我們需要一種分辨率高，且所需樣本量少于DFT的算法。
我們的意圖是將該算法用于長(cháng)度不確定的實(shí)時(shí)數據，因此采用了類(lèi)似于滑動(dòng)DFT的滑動(dòng)變換形式。
方法 算法要求
我們最初的目標是找到一種算法，即使數據本質(zhì)上是隨機且非平穩的，也能確定數據的潛在基波周期。初始算法要求如下：
? 能夠確定任何生物醫學(xué)信號（如ECG和SpO2）的基波周期。
? 響應時(shí)間足夠快，能夠實(shí)時(shí)跟蹤心臟心率周期和幅度的變化。
? 遭遇信號中斷、噪聲過(guò)大或運動(dòng)偽影時(shí)，能夠迅速恢復運行。
? 計算速度足夠快，以免成為確定采樣速率的限制因素。
? 對存儲空間的要求較低或適中，能夠在低功耗和便攜式設備中應用。
算法開(kāi)發(fā)
從DFT開(kāi)始，目標是找到周期，因此DFT方程中的頻率項被替換為周期，并且不是像DFT那樣逐步增加頻率，而是逐步增加周期。DFT以線(xiàn)性方式增加頻率，例如(1f0, 2f0, 3f0, …)，其中f0是第一諧波，而DPT則以采樣周期T0的倍數為單位，線(xiàn)性增大周期。盡管兩種算法的方程相似，但DFT無(wú)法產(chǎn)生與DPT相同的結果，因為兩種算法有本質(zhì)區別。通過(guò)分析描述其實(shí)現的方程，我們可以比較DFT和DPT。對于采樣頻率fS，N點(diǎn)DFT的頻點(diǎn)k對應頻率fK = k × fS / N Hz，公式1是樣本序列XI ... XI + N - 1的第k個(gè)頻點(diǎn)的頻譜表達式。
其中，k = 0, 1, 2, ...N - 1
DFT的第i個(gè)樣本按照公式2進(jìn)行計算。
如圖1所示，對于每個(gè)諧波，DFT基函數的縱坐標值與其之前第N個(gè)縱坐標值相同。發(fā)生這種情況的原因是，DFT中的所有諧波之間存在倍數關(guān)系，高次諧波是低次諧波的整數倍。
圖1.傅里葉變換正弦基函數，紅色所示為第一諧波(1 Hz)，藍色為第二諧波(2 Hz)，綠色為第三諧波(3 Hz)。
DPT中的項N必須針對每個(gè)周期進(jìn)行修改，因為周期之間并非簡(jiǎn)單的倍數關(guān)系，而是相差一個(gè)采樣周期，如圖2所示。
滑動(dòng)形式的DFT和DPT都需要實(shí)現循環(huán)或遞歸緩沖區，用于保存數量固定的最新樣本。當輸入數據為實(shí)數時(shí)，使用一個(gè)緩沖區；而當輸入數據為復數時(shí)，則使用兩個(gè)緩沖區。DPT變換的第i個(gè)樣本可以套入公式3。
其中，RBS為遞歸緩沖區大小，TL為最長(cháng)周期的長(cháng)度，TN為當前正在處理的基元的周期。這樣做可以使每個(gè)基礎周期的起始和終止縱坐標值相同。周期s從最小周期延伸到所選的最大周期，以覆蓋采樣數據中的周期。該實(shí)現利用了一組基函數，這些基函數代表了圖2中復正弦波的增量相位角。
圖2.三個(gè)相鄰正弦函數和三個(gè)相鄰余弦函數的周期變換復正弦基函數。此示例假設這些函數的采樣時(shí)間間隔為10 ms。
DPT的實(shí)現之所以有些困難，是因為基函數由多組復函數組成，這些復函數之間大多不是倍數關(guān)系，而且采樣周期不同。高效的DPT變換需使用圖3所示的基礎相位角。這也是本文所采用的實(shí)現形式。
圖3.周期變換基礎相位角，展示了復相位角的值如何隨著(zhù)每分鐘采樣周期數的增加而變化。上升曲線(xiàn)表示余弦相位角，下降曲線(xiàn)表示正弦相位角。
使用公式4可以輕松得出相量，其中“s”是以采樣周期為步長(cháng)，從最小選定周期到最大選定周期的周期集。
算法實(shí)現
滑動(dòng)DPT變換使用IIR濾波器實(shí)現，其信號流圖中在一個(gè)梳狀濾波器后接了一個(gè)諧振器，這與滑動(dòng)形式DFT的實(shí)現類(lèi)似。N個(gè)樣本的梳狀濾波器延遲導致瞬態(tài)響應的長(cháng)度為N-1個(gè)樣本。已經(jīng)有人使用心率調諧的梳狀濾波器并取得了一定的成功8。DPT復基函數或相位角的分量并非總是諧波相關(guān)，因此這些函數的端點(diǎn)不會(huì )在樣本空間中形成連續函數，這與DFT不同。然而，如果將DPT實(shí)現為滑動(dòng)變換，那么基函數就會(huì )被“包裹”起來(lái)，從而使基函數的分量變成連續的。當數據和基函數滑動(dòng)時(shí)，計算它們的相關(guān)性，基函數連續性得以保持。
在滑動(dòng)窗口算法中，長(cháng)度為N的窗口在長(cháng)度不確定的數據數組上滑動(dòng)。對于DPT而言，由于DPT可以處理實(shí)部和虛部?jì)深?lèi)輸入數據，因此需要維護兩個(gè)遞歸緩沖區。如果輸入只有一個(gè)實(shí)部（通常情況如此），則只需使用一個(gè)遞歸緩沖區。然而，根據輸入和基函數之間的相位關(guān)系，結果仍然可能是復數。結果存儲在兩個(gè)系綜緩沖區中，每個(gè)緩沖區的長(cháng)度為所選的最大周期。
MATLAB概念驗證模型
我們通過(guò)MATLAB腳本實(shí)現了公式4。圖4使用正弦和余弦函數作為輸入，幅度為±1，周期為45 ms、79 ms和175 ms。MATLAB腳本的周期限定在400 ms（200個(gè)周期/分鐘）到2 s（40個(gè)周期/分鐘）之間。本例總共處理了5000個(gè)數據樣本，樣本數量固定不變。由于輸入數據是幅度為1的正弦波形，因此每個(gè)周期的幅度也為1。很容易看出，這種變換實(shí)現的分辨率非常高。
圖4.幅度譜，展示了彼此不成倍數關(guān)系的三組輸入正弦數據的值。
圖5為每分鐘73個(gè)周期、幅度為4.5的正弦余弦波的結果。此示例使用了長(cháng)度為1500個(gè)數據點(diǎn)的遞歸緩沖區。請注意，存在一些較小誤差，幅度誤差為0.366%，周期誤差為0.234%。對于生物醫學(xué)應用而言，這些誤差的大小一般是可以接受的。在外周毛細血管血氧飽和度(SpO2)測量中，這些誤差無(wú)關(guān)緊要，因為SpO2是根據紅光和紅外光譜信號的比率之比來(lái)計算的9,10。參見(jiàn)公式6和公式7。
圖5.余弦波形的滑動(dòng)周期變換，每分鐘73個(gè)周期，振幅為4.5。幅度誤差小于0.37%，周期誤差小于0.24%。
結果 滑動(dòng)窗口DPT在脈搏血氧測定中的應用
將滑動(dòng)窗口算法應用于脈搏血氧測定時(shí)，為使算法正常運行，需要兩個(gè)遞歸數組：一個(gè)用于存儲紅光歷史記錄，另一個(gè)用于存儲紅外歷史記錄。為完成滑動(dòng)變換，還需根據相應周期的基函數，旋轉遞歸緩沖區（其長(cháng)度與正在處理的周期點(diǎn)相同）中更新的內容。該緩沖區的長(cháng)度決定了整體分辨率，一旦有足夠多的數據進(jìn)入處理流程以填充這些緩沖區，變換結果就會(huì )達到一個(gè)穩定的極限，只有幅度或周期會(huì )隨著(zhù)輸入數據的變化而改變。對于所報告的數據處理，遞歸緩沖區保存最后10秒的數據。
原始數據由ADI公司的研究人員收集，用于處理數據的軟件是MATLAB腳本中的滑動(dòng)DPT。圖6為從某位受試者獲取的原始數據；經(jīng)過(guò)1 Hz至4 Hz帶通濾波的數據，以及利用總寬度為200 ms的平坦光滑移動(dòng)平均濾波器處理后的數據。圖7為填充遞歸緩沖區之后頻譜達到穩定幅度的頻譜。隨著(zhù)新數據被采樣，DPT將持續跟蹤原始數據中的所有變化，頻譜也會(huì )相應地更新。
圖6.使用MAX30101 PPG AFE器件從某位受試者獲取的原始光電容積脈搏波數據、經(jīng)濾波的數據和經(jīng)平滑處理后的數據。上方波形表示原始紅外和紅光信號，而下方波形表示經(jīng)過(guò)濾波和平滑處理的數據。
圖7.此圖為采用滑動(dòng)窗口DPT處理的紅光和紅外光譜。兩個(gè)波峰中較大的是紅外光譜；較小的是紅光光譜。
為了估算SpO2，先需使用比率之比公式。交流分量使用圖7所示頻譜的峰值，直流分量使用圖6所示未濾波信號的平均值。
比較從Masimo血氧儀使用SET算法收集的SpO2和心率數據，與使用ADI MAX30101脈搏血氧儀傳感器同時(shí)獲取的數據。隨機選擇某位受試者的數據，并將結果繪制在圖8和圖9中。
圖8.DPT處理的光電容積脈搏波數據比較。
圖9.比較來(lái)自MAX30101血氧儀（采用離散周期變換進(jìn)行處理）和Masimo血氧儀的心率數據。
評估兩種不同儀器測量同一參數所產(chǎn)生的數值，是常見(jiàn)的醫學(xué)做法。其中一種儀器被認為能夠產(chǎn)生正確的結果，用作標準儀器。
Bland和Altman開(kāi)發(fā)了一種用于評估兩種定量測量結果一致性的方法11,12。他們通過(guò)分析平均差異和構建一致性界限來(lái)判斷一致性。Bland-Altman圖分析是評估平均差異之間的偏差和估計一致性區間的一種簡(jiǎn)單方法。如果對兩臺醫療儀器開(kāi)展此項測試，其中一臺被視為標準，則另一臺儀器的結果必須在標準儀器結果的兩個(gè)標準差或95%范圍內，才能認為其在臨床應用上與標準儀器效果相當。
與相關(guān)分析研究?jì)蓚€(gè)變量之間的關(guān)系不同，Bland-Altman方法是一種統計學(xué)方法，關(guān)注的是兩個(gè)變量之間的差異。
我們利用MAX30101脈搏血氧儀傳感器收集了26名健康成年受試者的數據，并將其與Masimo血氧儀（其中融合了新型信號提取技術(shù)Signal Extraction Technology?）的測量結果進(jìn)行比較，從而評估DPT算法的準確性和精確度13。研究對象包括15名男性和11名女性受試者，年齡在20至40歲之間。這項研究的目的是比較同一受試者使用兩種血氧儀的測量結果，而不是男性和女性之間的差異。請注意，兩性之間的SpO2確實(shí)略有不同。一項研究表明，對于年輕健康成年人，男性的平均SpO2為97.1±1.2%，而女性的平均SpO2為98.6±1.0%14。
圖10和圖11位使用Bland-Altman標準的結果，每個(gè)圓圈代表一位受試者的Bland-Altman結果。所有SpO2比較均符合Bland-Altman標準。
圖10.Masimo血氧儀與使用DPT算法的ADI血氧儀的SpO2百分比差異。滿(mǎn)足Bland-Altman標準。
圖11.Masimo血氧儀與使用DPT算法的ADI血氧儀的每分鐘心率差異。除一個(gè)案例外，其他所有案例均滿(mǎn)足Bland-Altman標準。箭頭標示了超出兩個(gè)標準差范圍的分析結果。
在圖11中，箭頭指向的心率比較值超出了兩個(gè)標準差范圍。該受試者的心率與時(shí)間關(guān)系圖如圖12所示，其中Masimo血氧儀的標準差為1.7892，而使用DPT算法的MAX30101血氧儀的標準差為0.8935。在這種情況下，我們很難確定哪種儀器更準確，但可以從標準差中找到一些線(xiàn)索。
圖12.Masimo血氧儀和ADI血氧儀的心率與時(shí)間的關(guān)系圖。在25秒周期內，Masimo血氧儀的標準差為1.7892，而MAX30101的標準差為0.8935。階梯波形是來(lái)自Masimo血氧儀的信號；平滑信號來(lái)自運行DPT算法的ADI血氧儀。
采用SDPT算法的血氧儀系統原型
最后，我們采用Arm?微處理器（運行裸機操作系統），設計了一個(gè)血氧儀原型。使用樹(shù)莓派Zero作為計算機平臺，MAX30102集成電路用作傳感器。操作系統和滑動(dòng)窗口DPT采用標準C語(yǔ)言實(shí)現。圖13即為該原型。整個(gè)血氧儀由USB 3.0連接供電。兩個(gè)數模轉換器根據監控軟件的判斷，通過(guò)帶狀電纜將數據發(fā)送到Tektronix DPO-4034示波器，在其中繪制圖像。然后，圖像通過(guò)網(wǎng)絡(luò )連接發(fā)送到臺式計算機。圖15為大約9秒的時(shí)間內從單個(gè)受試者獲得的結果，之后用10秒的時(shí)間來(lái)填充遞歸緩沖區。
圖13.基于樹(shù)莓派Zero的脈搏血氧儀原型。MAX30102 SpO2傳感器位于圖片左上角所示的指夾中。
通過(guò)一階低通IIR濾波器從原始信號中提取紅光和紅外直流信號；通過(guò)一階高通IIR濾波器提取交流信號。參見(jiàn)圖14。這些濾波器的時(shí)間常數設置為大約1秒。數據以100 SPS的速率采樣，并以MAX30102的中斷作為定時(shí)信號。對于紅光和紅外信號，該器件的輸出均為12位定點(diǎn)數字格式。
圖14.使用無(wú)限脈沖響應(IIR)濾波器從原始光譜數據中提取交流信號和直流信號。
圖15.樹(shù)莓派Zero血氧儀原型產(chǎn)生的PPG波形，上方為紅光脈沖，下方為紅外脈沖。心率約為58 bpm。圖中所示為倒置的波形，以便更準確地表示手指中的實(shí)際動(dòng)脈壓力。
紅光和紅外交流信號通過(guò)濾波器提取出來(lái)之后，就交由DPT處理，而無(wú)需任何進(jìn)一步的信號預處理。光譜信號的第一諧波產(chǎn)生的峰值如圖16所示。心率由橫坐標上數據峰值的位置決定，而SpO2通過(guò)比率之比公式使用紅光和紅外數據峰值的幅度直接計算。
圖16.樹(shù)莓派Zero使用滑動(dòng)窗口離散周期變換生成的頻譜，SpO2值為97%，心率為58 bpm。光標b（中心垂直藍線(xiàn)）顯示測量的心跳周期為1.03秒。左上角的矩形信號指向橫坐標上400 ms周期的位置；右上角的矩形信號指向橫坐標上2000 ms周期的位置。
討論
血氧儀產(chǎn)生的原始光信號包含較大的穩定直流分量和較小的振蕩交流分量，后者約為直流信號的1%。這些振蕩分量反映的是毛細血管中的脈動(dòng)活動(dòng)。任何運動(dòng)或其他偽影都可能輕易覆蓋這些信號，使讀數不準確。多年來(lái)人們花費了大量時(shí)間來(lái)研究將這些信號與偽影分離的方法。事實(shí)證明，這些方法通常非常復雜且難以實(shí)施16,17。
正是出于這些原因，我們才開(kāi)展了這項研究。DPT算法采用的變換只需少量樣本，但卻能實(shí)現準確的測量，許多挑戰因此迎刃而解。在周期域內進(jìn)行測量，并按采樣周期將每個(gè)周期點(diǎn)分隔開(kāi)來(lái)，便能提供所需的分辨率。然后，我們可以利用來(lái)自DPT的周期和幅度信息直接計算心率和血氧飽和度，而無(wú)需返回時(shí)域。結論
采用增量DPT算法的周期域分析，是處理周期性生物醫學(xué)信號以獲得頻譜成分的有效方法。該方法支持頻域分析，而且在實(shí)現上也有優(yōu)勢。研究表明，運行DPT算法的ADI MAX30101集成電路傳感器足夠精確，在醫療實(shí)踐中能夠取代Masimo血氧儀。
參考文獻
1 Amal Jubran?！癙ulse Oximetry”?！禖ritical Care》，第3卷，第2期，1999年2月。
2 Han-Wook Lee、Ju-Won Lee、Won-Geun Jun和Gun-Ki Lee?！癟he Periodic Moving Average Filter for Removing Motion Artifacts from PPG Signals”?！秶H控制自動(dòng)化與系統雜志》，第5卷，第6期，2007年12月。
3 Brendan Conlon、James A. Devine和James A. Dittmar?！癊CG Synchronized Pulse Oximeter”。美國專(zhuān)利4,960,126，1990年10月。
4 James Reuss和Dennis Bahr?！癙eriod Domain Analysis in Fetal Pulse Oximetry”。Proceedings of the Second Joint 24th Annual Conference and the Annual Fall Meeting of the Biomedical Engineering Society，2002年。
5 Eric Jacobsen和Richard Lyons?！癟he Sliding DFT”?！禝EEE信號處理雜志》，第20卷，第2期，2003年3月。
6 Eric Jacobsen和Richard Lyons?！癆n Update to the Sliding DFT”?！禝EEE信號處理雜志》，第21卷，第1期，2004年1月。
7 Lawrence R. Rabiner和Bernard Gold。Theory and Application of Digital Signal Processing。Prentice-Hall，1975年1月。
8 Ludvik Alkhoury、Ji-Won Choi、Chizhong Wang、Arjun Rajasekar、Sayandeep Acharya、Sean Mahoney、Barry S.Shender、Leonid Hrebien和Mose Kam?！癏eart-Rate Tuned Comb Filters For Processing Photoplethysmogram (PPG) Signals in Pulse Oximetry”?！杜R床監測與計算雜志》，第35卷，第4期，2021年8月。
9 “Recommended Configurations and Operating Profiles for MAX30101/MAX30102 EV Kits ”。 Maxim Integrated，2018年3月。
10 Sang-Soo Oak和Praveen Aroul?！癏ow to Design Peripheral Oxygen Saturation (SpO2) and Optical Heart Rate Monitoring (OHRM) Systems Using the AFE4403”。德州儀器，2015年3月。
11 Douglas Altman和J. Martin Bland?！癕easurement in Medicine:The Analysis of Method Comparison Studies”?！痘始医y計學(xué)會(huì )志，D輯:統計學(xué)家》，第32卷，第3期，1983年9月。
12 J.Martin Bland和Douglas G. Altman?！癝tatistical Methods for Assessing Agreement Between Two Methods of Clinical Measurement”?！读~刀》，1986年2月。
13 Julian M. Goldman、Michael T. Petterson、Robert J. Kopotic和Steven J. Barker?！癕asimo Signal Extraction Pulse Oximetry”?！杜R床監測與計算雜志》，第16卷，第7期，2000年。
14 Sagi Levental、Elie Picard、Francis Mimouni、Leon Joseph、Tal Y Samuel、Reuben Bromiker、Dror Mandel、Nissim Arish和Shmuel Goldberg?！癝ex-Linked Difference in Blood Oxygen Saturation”?！杜R床呼吸雜志》，第12卷，第5期，2018年5月。
15 Sam Koblenski?！癊veryday DSP for Programmers: DC and Impulsive NoiseRemoval”。2015年11月。
16 Surekha Palreddy?！癝ignal Processing Algorithms”。Design of Pulse Oximeters，第一版，1997年10月
17 Terry L. Rusch、Ravi Sankar和John E. Scharf?！癝ignal Processing Methods for Pulse Oximetry ”。 《生物學(xué)和醫學(xué)中的計算機雜志》，第26卷，第2期，1996年3月。
（來(lái)源：ADI公司，作者：Dennis E. Bahr博士，Bahr Management, Inc.總裁兼生物醫學(xué)工程師，Marc Smith，首席工程師）
免責聲明：本文為轉載文章，轉載此文目的在于傳遞更多信息，版權歸原作者所有。本文所用視頻、圖片、文字如涉及作品版權問(wèn)題，請聯(lián)系小編進(jìn)行處理。
推薦閱讀：
無(wú)纜智能終端的能源進(jìn)化論：破解微型設備供能困局
L Nanopower革新智能家居能源架構：nA級功耗技術(shù)破解無(wú)線(xiàn)終端續航困境
學(xué)子專(zhuān)區論壇 - ADALM2000實(shí)驗：Hartley振蕩器
