【導讀】生理信號固有的準周期性特征表現(xiàn)為時變幅頻特性，其有效信息常與運動偽影及環(huán)境噪聲在頻域高度重疊。傳統(tǒng)頻域濾波手段因信號非平穩(wěn)性面臨失效風險，業(yè)界普遍采用時間鎖相式系綜平均法——通過ECG心電觸發(fā)信號建立時域參考基準，實現(xiàn)血氧參數(shù)的有效提取。然而臨床場景中ECG信號的可及性限制催生了技術盲區(qū)。本研究突破性開發(fā)自主周期定位算法，在不依賴外部心電參考的條件下，達成與ECG引導方案相當?shù)男盘栔貥嬀?，為可穿戴醫(yī)療設備提供了去ECG依賴的生理參數(shù)解析新范式。
摘要
本文介紹了新型滑動離散周期變換(DPT)算法，可設計用于處理生理信號，尤其是脈搏血氧儀采集的光電容積脈搏波(PPG)信號。該算法采用正弦基函數(shù)進行周期域分析，可解決隨機噪聲和非平穩(wěn)數(shù)據(jù)等難題。DPT在MATLAB?中作為滑動變換實現(xiàn)，結合了自相關與系綜平均。文中將詳細介紹在ADI MAX30101器件上開發(fā)和實現(xiàn)的一種算法，并與采用Signal Extraction Technology? (SET)的Masimo血氧儀進行比較。
簡介
生理信號固有的準周期性特征表現(xiàn)為時變幅頻特性，其有效信息常與運動偽影及環(huán)境噪聲在頻域高度重疊。傳統(tǒng)頻域濾波手段因信號非平穩(wěn)性面臨失效風險，業(yè)界普遍采用時間鎖相式系綜平均法——通過ECG心電觸發(fā)信號建立時域參考基準，實現(xiàn)血氧參數(shù)的有效提取。然而臨床場景中ECG信號的可及性限制催生了技術盲區(qū)。本研究突破性開發(fā)自主周期定位算法，在不依賴外部心電參考的條件下，達成與ECG引導方案相當?shù)男盘栔貥嬀?，為可穿戴醫(yī)療設備提供了去ECG依賴的生理參數(shù)解析新范式。
最初，我們開發(fā)了一種算法來執(zhí)行某種形式的自相關和系綜平均處理4。然而，我們很快發(fā)現(xiàn)，時域中的系綜平均并無必要，因為所有相關的信息都可以在周期域數(shù)據(jù)本身中找到。心率和血氧飽和度可以直接根據(jù)滑動離散周期變換(DPT)產生的結果計算出來。
這項工作始于對離散傅里葉變換(DFT)的回顧，因為DFT能夠生成信號的頻譜，然后可以利用頻譜確定其周期5,6。該研究的另一個目標是以非常高的分辨率進行數(shù)據(jù)采樣。為了利用DFT實現(xiàn)高分辨率，需要收集大量數(shù)據(jù)樣本。由于生物信號具有準平穩(wěn)性，使用DFT收集大量樣本常常會導致頻譜模糊7。我們需要一種分辨率高，且所需樣本量少于DFT的算法。
我們的意圖是將該算法用于長度不確定的實時數(shù)據(jù)，因此采用了類似于滑動DFT的滑動變換形式。
方法 算法要求
我們最初的目標是找到一種算法，即使數(shù)據(jù)本質上是隨機且非平穩(wěn)的，也能確定數(shù)據(jù)的潛在基波周期。初始算法要求如下：
? 能夠確定任何生物醫(yī)學信號（如ECG和SpO2）的基波周期。
? 響應時間足夠快，能夠實時跟蹤心臟心率周期和幅度的變化。
? 遭遇信號中斷、噪聲過大或運動偽影時，能夠迅速恢復運行。
? 計算速度足夠快，以免成為確定采樣速率的限制因素。
? 對存儲空間的要求較低或適中，能夠在低功耗和便攜式設備中應用。
算法開發(fā)
從DFT開始，目標是找到周期，因此DFT方程中的頻率項被替換為周期，并且不是像DFT那樣逐步增加頻率，而是逐步增加周期。DFT以線性方式增加頻率，例如(1f0, 2f0, 3f0, …)，其中f0是第一諧波，而DPT則以采樣周期T0的倍數(shù)為單位，線性增大周期。盡管兩種算法的方程相似，但DFT無法產生與DPT相同的結果，因為兩種算法有本質區(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ù)關系，高次諧波是低次諧波的整數(shù)倍。
圖1.傅里葉變換正弦基函數(shù)，紅色所示為第一諧波(1 Hz)，藍色為第二諧波(2 Hz)，綠色為第三諧波(3 Hz)。
DPT中的項N必須針對每個周期進行修改，因為周期之間并非簡單的倍數(shù)關系，而是相差一個采樣周期，如圖2所示。
滑動形式的DFT和DPT都需要實現(xiàn)循環(huán)或遞歸緩沖區(qū)，用于保存數(shù)量固定的最新樣本。當輸入數(shù)據(jù)為實數(shù)時，使用一個緩沖區(qū)；而當輸入數(shù)據(jù)為復數(shù)時，則使用兩個緩沖區(qū)。DPT變換的第i個樣本可以套入公式3。
其中，RBS為遞歸緩沖區(qū)大小，TL為最長周期的長度，TN為當前正在處理的基元的周期。這樣做可以使每個基礎周期的起始和終止縱坐標值相同。周期s從最小周期延伸到所選的最大周期，以覆蓋采樣數(shù)據(jù)中的周期。該實現(xiàn)利用了一組基函數(shù)，這些基函數(shù)代表了圖2中復正弦波的增量相位角。
圖2.三個相鄰正弦函數(shù)和三個相鄰余弦函數(shù)的周期變換復正弦基函數(shù)。此示例假設這些函數(shù)的采樣時間間隔為10 ms。
DPT的實現(xiàn)之所以有些困難，是因為基函數(shù)由多組復函數(shù)組成，這些復函數(shù)之間大多不是倍數(shù)關系，而且采樣周期不同。高效的DPT變換需使用圖3所示的基礎相位角。這也是本文所采用的實現(xiàn)形式。
圖3.周期變換基礎相位角，展示了復相位角的值如何隨著每分鐘采樣周期數(shù)的增加而變化。上升曲線表示余弦相位角，下降曲線表示正弦相位角。
使用公式4可以輕松得出相量，其中“s”是以采樣周期為步長，從最小選定周期到最大選定周期的周期集。
算法實現(xiàn)
滑動DPT變換使用IIR濾波器實現(xiàn)，其信號流圖中在一個梳狀濾波器后接了一個諧振器，這與滑動形式DFT的實現(xiàn)類似。N個樣本的梳狀濾波器延遲導致瞬態(tài)響應的長度為N-1個樣本。已經有人使用心率調諧的梳狀濾波器并取得了一定的成功8。DPT復基函數(shù)或相位角的分量并非總是諧波相關，因此這些函數(shù)的端點不會在樣本空間中形成連續(xù)函數(shù)，這與DFT不同。然而，如果將DPT實現(xiàn)為滑動變換，那么基函數(shù)就會被“包裹”起來，從而使基函數(shù)的分量變成連續(xù)的。當數(shù)據(jù)和基函數(shù)滑動時，計算它們的相關性，基函數(shù)連續(xù)性得以保持。
在滑動窗口算法中，長度為N的窗口在長度不確定的數(shù)據(jù)數(shù)組上滑動。對于DPT而言，由于DPT可以處理實部和虛部兩類輸入數(shù)據(jù)，因此需要維護兩個遞歸緩沖區(qū)。如果輸入只有一個實部（通常情況如此），則只需使用一個遞歸緩沖區(qū)。然而，根據(jù)輸入和基函數(shù)之間的相位關系，結果仍然可能是復數(shù)。結果存儲在兩個系綜緩沖區(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ù)關系的三組輸入正弦數(shù)據(jù)的值。
圖5為每分鐘73個周期、幅度為4.5的正弦余弦波的結果。此示例使用了長度為1500個數(shù)據(jù)點的遞歸緩沖區(qū)。請注意，存在一些較小誤差，幅度誤差為0.366%，周期誤差為0.234%。對于生物醫(yī)學應用而言，這些誤差的大小一般是可以接受的。在外周毛細血管血氧飽和度(SpO2)測量中，這些誤差無關緊要，因為SpO2是根據(jù)紅光和紅外光譜信號的比率之比來計算的9,10。參見公式6和公式7。
圖5.余弦波形的滑動周期變換，每分鐘73個周期，振幅為4.5。幅度誤差小于0.37%，周期誤差小于0.24%。
結果 滑動窗口DPT在脈搏血氧測定中的應用
將滑動窗口算法應用于脈搏血氧測定時，為使算法正常運行，需要兩個遞歸數(shù)組：一個用于存儲紅光歷史記錄，另一個用于存儲紅外歷史記錄。為完成滑動變換，還需根據(jù)相應周期的基函數(shù)，旋轉遞歸緩沖區(qū)（其長度與正在處理的周期點相同）中更新的內容。該緩沖區(qū)的長度決定了整體分辨率，一旦有足夠多的數(shù)據(jù)進入處理流程以填充這些緩沖區(qū)，變換結果就會達到一個穩(wěn)定的極限，只有幅度或周期會隨著輸入數(shù)據(jù)的變化而改變。對于所報告的數(shù)據(jù)處理，遞歸緩沖區(qū)保存最后10秒的數(shù)據(jù)。
原始數(shù)據(jù)由ADI公司的研究人員收集，用于處理數(shù)據(jù)的軟件是MATLAB腳本中的滑動DPT。圖6為從某位受試者獲取的原始數(shù)據(jù)；經過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ù)、經濾波的數(shù)據(jù)和經平滑處理后的數(shù)據(jù)。上方波形表示原始紅外和紅光信號，而下方波形表示經過濾波和平滑處理的數(shù)據(jù)。
圖7.此圖為采用滑動窗口DPT處理的紅光和紅外光譜。兩個波峰中較大的是紅外光譜；較小的是紅光光譜。
為了估算SpO2，先需使用比率之比公式。交流分量使用圖7所示頻譜的峰值，直流分量使用圖6所示未濾波信號的平均值。
比較從Masimo血氧儀使用SET算法收集的SpO2和心率數(shù)據(jù)，與使用ADI MAX30101脈搏血氧儀傳感器同時獲取的數(shù)據(jù)。隨機選擇某位受試者的數(shù)據(jù)，并將結果繪制在圖8和圖9中。
圖8.DPT處理的光電容積脈搏波數(shù)據(jù)比較。
圖9.比較來自MAX30101血氧儀（采用離散周期變換進行處理）和Masimo血氧儀的心率數(shù)據(jù)。
評估兩種不同儀器測量同一參數(shù)所產生的數(shù)值，是常見的醫(yī)學做法。其中一種儀器被認為能夠產生正確的結果，用作標準儀器。
Bland和Altman開發(fā)了一種用于評估兩種定量測量結果一致性的方法11,12。他們通過分析平均差異和構建一致性界限來判斷一致性。Bland-Altman圖分析是評估平均差異之間的偏差和估計一致性區(qū)間的一種簡單方法。如果對兩臺醫(yī)療儀器開展此項測試，其中一臺被視為標準，則另一臺儀器的結果必須在標準儀器結果的兩個標準差或95%范圍內，才能認為其在臨床應用上與標準儀器效果相當。
與相關分析研究兩個變量之間的關系不同，Bland-Altman方法是一種統(tǒng)計學方法，關注的是兩個變量之間的差異。
我們利用MAX30101脈搏血氧儀傳感器收集了26名健康成年受試者的數(shù)據(jù)，并將其與Masimo血氧儀（其中融合了新型信號提取技術Signal Extraction Technology?）的測量結果進行比較，從而評估DPT算法的準確性和精確度13。研究對象包括15名男性和11名女性受試者，年齡在20至40歲之間。這項研究的目的是比較同一受試者使用兩種血氧儀的測量結果，而不是男性和女性之間的差異。請注意，兩性之間的SpO2確實略有不同。一項研究表明，對于年輕健康成年人，男性的平均SpO2為97.1±1.2%，而女性的平均SpO2為98.6±1.0%14。
圖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算法的血氧儀系統(tǒng)原型
最后，我們采用Arm?微處理器（運行裸機操作系統(tǒng)），設計了一個血氧儀原型。使用樹莓派Zero作為計算機平臺，MAX30102集成電路用作傳感器。操作系統(tǒng)和滑動窗口DPT采用標準C語言實現(xiàn)。圖13即為該原型。整個血氧儀由USB 3.0連接供電。兩個數(shù)模轉換器根據(jù)監(jiān)控軟件的判斷，通過帶狀電纜將數(shù)據(jù)發(fā)送到Tektronix DPO-4034示波器，在其中繪制圖像。然后，圖像通過網(wǎng)絡連接發(fā)送到臺式計算機。圖15為大約9秒的時間內從單個受試者獲得的結果，之后用10秒的時間來填充遞歸緩沖區(qū)。
圖13.基于樹莓派Zero的脈搏血氧儀原型。MAX30102 SpO2傳感器位于圖片左上角所示的指夾中。
通過一階低通IIR濾波器從原始信號中提取紅光和紅外直流信號；通過一階高通IIR濾波器提取交流信號。參見圖14。這些濾波器的時間常數(shù)設置為大約1秒。數(shù)據(jù)以100 SPS的速率采樣，并以MAX30102的中斷作為定時信號。對于紅光和紅外信號，該器件的輸出均為12位定點數(shù)字格式。
圖14.使用無限脈沖響應(IIR)濾波器從原始光譜數(shù)據(jù)中提取交流信號和直流信號。
圖15.樹莓派Zero血氧儀原型產生的PPG波形，上方為紅光脈沖，下方為紅外脈沖。心率約為58 bpm。圖中所示為倒置的波形，以便更準確地表示手指中的實際動脈壓力。
紅光和紅外交流信號通過濾波器提取出來之后，就交由DPT處理，而無需任何進一步的信號預處理。光譜信號的第一諧波產生的峰值如圖16所示。心率由橫坐標上數(shù)據(jù)峰值的位置決定，而SpO2通過比率之比公式使用紅光和紅外數(shù)據(jù)峰值的幅度直接計算。
圖16.樹莓派Zero使用滑動窗口離散周期變換生成的頻譜，SpO2值為97%，心率為58 bpm。光標b（中心垂直藍線）顯示測量的心跳周期為1.03秒。左上角的矩形信號指向橫坐標上400 ms周期的位置；右上角的矩形信號指向橫坐標上2000 ms周期的位置。
討論
血氧儀產生的原始光信號包含較大的穩(wěn)定直流分量和較小的振蕩交流分量，后者約為直流信號的1%。這些振蕩分量反映的是毛細血管中的脈動活動。任何運動或其他偽影都可能輕易覆蓋這些信號，使讀數(shù)不準確。多年來人們花費了大量時間來研究將這些信號與偽影分離的方法。事實證明，這些方法通常非常復雜且難以實施16,17。
正是出于這些原因，我們才開展了這項研究。DPT算法采用的變換只需少量樣本，但卻能實現(xiàn)準確的測量，許多挑戰(zhàn)因此迎刃而解。在周期域內進行測量，并按采樣周期將每個周期點分隔開來，便能提供所需的分辨率。然后，我們可以利用來自DPT的周期和幅度信息直接計算心率和血氧飽和度，而無需返回時域。結論
采用增量DPT算法的周期域分析，是處理周期性生物醫(yī)學信號以獲得頻譜成分的有效方法。該方法支持頻域分析，而且在實現(xiàn)上也有優(yōu)勢。研究表明，運行DPT算法的ADI MAX30101集成電路傳感器足夠精確，在醫(yī)療實踐中能夠取代Masimo血氧儀。
參考文獻
1 Amal Jubran?！癙ulse Oximetry”。《Critical 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控制自動化與系統(tǒng)雜志》，第5卷，第6期，2007年12月。
3 Brendan Conlon、James A. Devine和James A. Dittmar?！癊CG Synchronized Pulse Oximeter”。美國專利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”。《臨床監(jiān)測與計算雜志》，第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(tǒng)計學會志，D輯:統(tǒng)計學家》，第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床監(jiān)測與計算雜志》，第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 ”。 《生物學和醫(yī)學中的計算機雜志》，第26卷，第2期，1996年3月。
（來源：ADI公司，作者：Dennis E. Bahr博士，Bahr Management, Inc.總裁兼生物醫(yī)學工程師，Marc Smith，首席工程師）
免責聲明：本文為轉載文章，轉載此文目的在于傳遞更多信息，版權歸原作者所有。本文所用視頻、圖片、文字如涉及作品版權問題，請聯(lián)系小編進行處理。
推薦閱讀：
破解運放穩(wěn)定性謎題：工程師必備的穩(wěn)定性設計手冊
L Nanopower革新智能家居能源架構：nA級功耗技術破解無線終端續(xù)航困境
學子專區(qū)論壇 - ADALM2000實驗：Hartley振蕩器
