新型離散週期轉換方法如何處理生理訊號

作者:ADI Dennis Bahr博士 及應用技術人員 Marc Smith


問題:

滑動離散週期轉換(DPT)演算法如何處理即時應用的心率和訊號品質變化?

How the New Discrete Period Transform Method Can Process Physiological Signals

答案:

DPT演算法的目的在於適應即時應用中心率和訊號品質的變化。其採用滑動視窗方法,每次有新樣本時,都會分析一個新的資料視窗。透過如此的演算法,就能追蹤訊號隨時間的變化,並能相應地調整其分析。此外,DPT演算法結合了自相關與總合平均技術,可增強訊號雜訊比並減輕偽影的影響,有助於確保即使在惡劣條件下也能獲得精準可靠的結果。

簡介

生理來源的訊號可能受到雜訊和運動偽影的干擾,這些干擾的通帶可能與訊號本身相互重疊1。生物訊號具有準靜態性,其週期和幅度會隨時間而變化2 。對於此類訊號,無法透過簡單的資料濾波進行處理。為了擷取資訊,一種常用方法是使用與目標訊號時間同步的訊號,作為時間參考,來進行總合平均。依靠ECG 源的外部心臟觸發訊號,總合平均方法可以有效地處理血氧訊號3 ,但在許多情況下,可能無法獲取ECG源。本文在沒有ECG觸發訊號的前提下成功地處理了血氧訊號,並獲得了類似的結果。

最初,我們開發了一種演算法來執行某種形式的自相關和總合平均操作4。然而,我們很快地發現,時域中的總合平均並無必要,因為所有相關的資訊都可以在週期域資料本身中找到。心率和血氧飽和度可以直接根據滑動離散週期變換(DPT)產生的結果來計算。

這項工作始於對離散傅立葉變換(DFT)的回顧,因為DFT能夠產生訊號的頻譜,然後可以利用頻譜確定其週期5,6。該研究的另一個目標,是以非常高的解析度進行資料採樣。為了利用DFT實現高解析度,需要收集大量資料樣本。由於生物訊號具有準靜態性,使用DFT收集大量樣本常常會導致頻譜模糊7,因此我們需要一種解析度高,且所需樣本量少於DFT的演算法,由於目的是將該演算法用於長度不確定的即時資料,因此採用了類似於滑動DFT的滑動變換形式。

方法

演算法要求

最初的目標是找到一種演算法,即使資料本質上是隨機且非平穩的,也能確定資料的潛在基波週期。初始演算法要求如下:

  • 能夠確定任何生物醫學訊號(如ECG和SpO2)的基波週期。
  • 回應時間夠快,能夠即時追蹤心臟心率週期和幅度的變化。
  • 遭遇訊號中斷、雜訊過大或運動偽影時,能夠迅速恢復運行。
  • 運算速度夠快,以免成為確定採樣速率的限制因素。
  • 對於儲存空間的要求較低或適中,能夠在低功耗和可攜式裝置中應用。

演算法開發

從DFT開始,目標是找到週期,因此DFT方程式中的頻率項被替換為週期,並且不是像DFT逐步增加頻率,而是逐步增加週期。DFT以線性方式增加頻率,例如(1f0, 2f0, 3f0, …),其中f0是第一諧波,而DPT則以採樣週期 T0的倍數為單位,線性增大週期。儘管兩種演算法的方程式相似,但DFT無法產生與DPT相同的結果,因為兩種演算法有其本質上的區別。透過分析描述其實現的方程式,我們可以比較DFT和DPT。對於採樣頻率 fS,N點DFT的頻點k對應頻率 fKS / N Hz,公式1是樣本序列XI... XI + N - 1的第k個頻點的頻譜運算式。

Equation 1.

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

DFT的第i個樣本按照公式2進行計算。

Equation 2.

如圖1所示,對於每個諧波,DFT基函數的縱座標值與其之前第N個縱座標值相同。產生此種情況的原因是DFT中的所有諧波之間存在倍數關係,高次諧波是低次諧波的整數倍。

Figure 1. Fourier transform sine basis functions where the red is the first harmonic (1 Hz), blue is the second (2 Hz), and green is the third (3 Hz).
圖1. 傅立葉變換正弦基函數,紅色所示為第一諧波(1 Hz),藍色為第二諧 波(2 Hz),綠色為第三諧波(3 Hz)。

DPT中的項N必須針對每個週期進行修改,因為週期之間並非簡單的倍數關係,而是相差一個採樣週期,如圖2所示。

Figure 2. Period transform complex sinusoidal basis functions for three adjacent sin and three adjacent cos functions. This example assumes that these functions are spaced by a sample period of 10 ms.
圖2. 三個相鄰正弦函數和三個相鄰餘弦函數的週期變換複正弦基函數。 此示例假設這些函數的採樣時間間隔為10 ms。

滑動形式的DFT和DPT都需要實現迴圈或遞迴緩衝區,用於保存數量固定的最新樣本。當輸入資料為實數時使用一個緩衝區;而當輸入資料為複數時,則使用兩個緩衝區。DPT變換的第i個樣本可以套入公式3。

Equation 3.

其中,RBS為遞迴緩衝區大小, TL為最長週期的長度,TN為目前正在處理的基元的週期。如此做法可以使每個基礎週期的起始和終止縱座標值相同。週期s從最小週期延伸到所選的最大週期,以覆蓋採樣資料中的週期。該實現利用了一組基函數,這些基函數代表了圖2中複正弦波的增量相位角。

DPT的實現之所以有些困難,是因為基函數由多組復函陣列成,這些複函數之間大多不是倍數關係,而且採樣週期不同。高效的DPT變換需使用圖3所示的基礎相位角。這也是本文所採用的實現形式。

Figure 3. Period transform basis phase angles showing the values of the complex phase angles for increasing number of sample periods per minute. The ascending curve shows the cosine phase angles and the descending curve shows the sine phase angles.
圖3. 週期變換基礎相位角,展示了複相位角的值如何隨著每分鐘採樣週期數的增加而變化。上升曲線表示餘弦相位角,下降曲線表示正弦相位角。

使用公式4可以輕鬆得出相量,其中“s”是以採樣週期為步驟,從最小選定週期到最大選定週期的週期集。

Equation 4.

演算法實現

滑動DPT變換使用IIR濾波器實現,其訊號流程圖中在一個梳狀濾波器後接了一個諧振器,這與滑動形式DFT的實現類似。N個樣本的梳狀濾波器延遲導致瞬態響應的長度為N-1個樣本。已經有人使用心率調諧的梳狀濾波器並取得了一定的成功8。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。在此顯而易見的是,此種變換實現的解析度非常高。

Figure 4. Amplitude spectrum showing the values from three sets of input sinusoidal data that are incommensurate relative to each other.
圖4. 幅度譜,展示了彼此不成倍數關係的三組輸入正弦資料的值。

圖5為每分鐘73個週期、幅度為4.5的正弦餘弦波的結果。此示例使用了長度為1500個資料點的遞迴緩衝區。請注意其中存在一些較小誤差,幅度誤差為0.366%,週期誤差為0.234%。對於生物醫學應用而言,這些誤差的大小一般是可以接受的。在周邊毛細孔血管血氧飽和度(SpO2)測量中,這些誤差並不重要,因為SpO2是根據紅光和紅外光譜訊號的比率之比來計算的9,10。參見公式6和公式7。

Figure 5. The sliding period transform for a cosine waveform with a period of 73 periods per minute and amplitude of 4.5. The amplitude error is less than 0.37% and the period error is less than 0.24%.
圖5. 餘弦波形的滑動週期變換,每分鐘73個週期,振幅為4.5。幅度誤差小於0.37%,週期誤差小於0.24%。

結果

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

將滑動視窗演算法應用於脈搏血氧測定時,為使演算法正常運行,需要兩個遞迴陣列:一個用於儲存紅外線歷史記錄,另一個用於儲存紅外線歷史記錄。為完成滑動變換,還需根據相應週期的基函數,旋轉遞迴緩衝區(其長度與正在處理的週期點相同)中更新的內容。該緩衝區的長度決定了整體解析度,一旦有夠多的資料進入處理流程以填充這些緩衝區,變換結果就會達到一個穩定的極限,只有幅度或週期會隨著輸入資料的變化而改變。對於所報告的資料處理,遞迴緩衝區保存最後10秒的資料。

原始資料由ADI研究人員收集,用於處理資料的軟體是MATLAB腳本中的滑動DPT。圖6為從某位受試者獲取的原始資料;經過1 Hz至4 Hz帶通濾波的資料,以及利用總寬度為200 ms的平坦光滑移動平均濾波器處理後的資料。圖7為填充遞迴緩衝區之後頻譜達到穩定幅度的頻譜。隨著新資料被採樣,DPT將持續追蹤原始資料中的所有變化,頻譜也會相應地更新。

Figure 6. Raw, filtered, and smoothed photoplethysmographic pulse data from a subject using a MAX30101 PPG AFE device. The top waveforms show the raw infrared and red signals while the bottom waveforms show the filtered and smoothed data.
圖6. 使用 MAX30101 PPG AFE元件從某位受試者獲取的原始光電容積脈搏波數 據、經濾波的資料和經平滑處理後的資料。上方波形表示原始紅外線和紅色波段訊號,而下方波形表示經過濾波和平滑處理的資料。
Figure 7. This graph shows the red and infrared spectrums using the sliding window DPT. The larger of the two peaks is the infrared spectrum and the smaller is the red spectrum.
圖7. 此圖為採用滑動視窗DPT處理的紅外線和紅色波段。兩個波峰中較大的是紅外線光譜;較小的是紅色光譜。

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

Equation 5.

Equation 6.

Equation 7.

比較從Masimo血氧計使用SET演算法收集的SpO2和心率資料,與使用ADI MAX30101脈搏血氧計感測器同時擷取的資料。圖8和圖9繪製的結果為隨機選擇某位受試者的資料。

Figure 8. Comparison of DPT processed photoplethysmographic data.
圖8. DPT處理的光電容積脈搏波資料比較。
Figure 9. Comparison of processed heart rate data from the MAX30101 oximeter processed with the discrete period transform and a Masimo oximeter.
圖9. 比較來自MAX30101血氧計(採用離散週期變換進行處理)和Masimo血氧計的心率數據。

評估由兩種不同儀器來測量同一參數所產生的數值,是常見的醫學做法。其中一種儀器被認為能夠產生正確的結果,主要用於標準儀器。

Bland和Altman開發了一種用於評估兩種定量測量結果一致性的方法11,12。其透過分析平均差異和建構一致性界限來判斷一致性。Bland-Altman圖分析是評估平均差異之間的偏差和估計一致性區間的一種簡單方法。如果對兩台醫療儀器開展此項測試,其中一台被視為標準,則另一台儀器的結果,必須在標準儀器結果的兩個標準差或95%範圍內才能視為在臨床應用上與標準儀器效果相符。

與相關分析研究兩個變數之間的關係不同,Bland-Altman方法是一種統計學方法,所關注的是兩個變數之間的差異。

我們利用MAX30101脈搏血氧計感測器收集了26名健康成年受試者的資料,並將其與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標準。

Figure 10. The percentage differences of SpO2 between a Masimo oximeter and an ADI oximeter using the DPT algorithm. The Bland-Altman criteria was met.
圖10. Masimo血氧計與使用DPT演算法的ADI血氧計的SpO2 百分比差異。滿足Bland-Altman標準。
Figure 11. Heart rate differences in beats per minute between a Masimo oximeter and an ADI oximeter using the DPT algorithm. The Bland-Altman algorithm was met in all cases but one. An arrow shows where the analysis falls outside of two standard deviations.
圖11. Masimo血氧計與使用DPT演算法的ADI血氧計的每分鐘心率差異。除一個案例外,其他所有案例均滿足Bland-Altman標準。箭頭標示了超出兩個標準差範圍的分析結果。

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

Figure 12. Plot of heart rate vs. time for a Masimo oximeter and an ADI oximeter. For the 25 s period, the Masimo oximeter had a standard deviation of 1.7892 and the MAX30101 had a standard deviation of 0.8935. The stepped waveform is the signal from the Masimo oximeter and the smooth signal is from the ADI oximeter running the DPT algorithm.
圖12. Masimo血氧計和ADI血氧計的心率與時間的關係圖。在25秒週期內, Masimo血氧計的標準差為1.7892,而MAX30101的標準差為0.8935。階梯波形是來自Masimo血氧計的訊號;平滑訊號來自運行DPT演算法的ADI血氧計。

採用SDPT演算法的血氧計系統原型

最後,我們採用Arm®微處理器(運行裸機作業系統),設計了一個血氧計原型。使用樹莓派Zero作為電腦平台,並以 MAX30102 IC作為感測器。作業系統和滑動視窗DPT採用標準C語言實現。圖13即為該原型。整個血氧計由USB 3.0連接供電。兩個數位類比轉換器根據監控軟體的判斷,透過帶狀電纜將資料發送到Tektronix DPO-4034示波器,在其中繪製圖像。然後,圖像透過網路連接發送到桌上型電腦。圖15為大約是在9秒的時間內從單一受試者獲得的結果,之後則是用10秒的時間來填充遞迴緩衝區。

Figure 13. The prototype Raspberry Pi Zero-based pulse oximeter. The MAX30102 SpO2 sensor is in the finger clip shown in the upper left corner of the picture.
圖13. 基於樹莓派Zero的脈搏血氧計原型。MAX30102 SpO2感測器位於圖片左上角所示的指夾中。

透過一階低通IIR濾波器從原始訊號中擷取紅光和紅外線直流訊號,另外則是透過一階高通IIR濾波器提取交流訊號。參見圖14。這些濾波器的時間常數設定為大約1秒。資料以100 SPS的速率採樣,並以MAX30102的中斷作為定時訊號。對於紅光和紅外線訊號而言,該元件的輸出均為12位元定點數位格式。

Figure 14. The infinite impulse response (IIR) filters used to extract the AC signals and the DC signals from the raw spectrographic data.
圖14. 使用無限脈衝回應(IIR)濾波器從原始光譜資料中擷取交流訊號和直流訊號。
Figure 15. PPG waveforms from the Raspberry Pi Zero oximeter prototype showing the red pulses on the top and infrared pulses on the bottom. The heart rate is approximately 58 bpm. Inverted waveforms are shown to more accurately represent the actual arterial pressure in the finger.
圖15. 樹莓派Zero血氧計原型產生的PPG波形,上方為紅光脈衝,下方為紅外線脈衝。心率約為58 bpm。圖中所示為倒置的波形,以便更精準地表示手指中的實際動脈壓力。

紅光和紅外線交流訊號透過濾波器擷取出來之後,就可交由DPT處理,而無需任何進一步的訊號預處理。光譜訊號的第一諧波產生的峰值如圖16所示。心率由橫座標上資料峰值的位置決定,而 SpO2則是透過比率之比公式,使用紅光和紅外線資料峰值的幅度來直接計算。

Figure 16. The spectrums generated by Raspberry Pi Zero using the sliding window discrete period transform with an SpO2 value of 97% and a heart rate of 58 bpm. The cursor b (center vertical blue line) shows where the heart period is measured as 1.03 s. The rectangular signal at the top left points to where the 400 ms period is located on the abscissa and the rectangular signal at the top right points to where 2000 ms period is located on the abscissa.
圖16. 樹莓派Zero使用滑動視窗離散週期變換生成的頻譜,SpO2值為97%,心率為58 bpm。游標b(中心垂直藍線)顯示測量的心跳週期為1.03秒。左上角的矩形訊號指向橫座標上400 ms週期的位置;右上角的矩形訊號指向橫座標上2000 ms週期的位置。

討論

血氧計產生的原始光訊號包含較大的穩定直流成份和較小的振盪交流成份,後者約為直流訊號的1%。這些振盪分量反映的是毛細孔血管中的脈動活動。任何運動或其他偽影都可能輕易覆蓋這些訊號,影響讀數的精準度。多年來,人們花費了大量時間來研究將這些訊號與偽影分離的方法。事實證明,這些方法通常非常複雜且難以執行16,17

正是出於這些原因,我們才開展了這項研究。DPT演算法採用的變換只需少量樣本,但卻能實現精準的測量,許多挑戰因此迎刃而解。在週期域內進行測量,並按採樣週期將每個週期點分隔開來,便能提供所需的解析度。然後,我們可以利用來自DPT的週期和幅度資訊直接計算心率和血氧飽和度,而無需返回時域。

結論

採用增量DPT演算法的週期域分析是處理週期性生物醫學訊號以獲得頻譜成分的有效方法。該方法支援頻域分析,而且在實踐上也具有優勢。研究顯示,運行DPT演算法的ADI MAX30101 IC感測器擁有足夠精準度,而能在醫療實踐中取代Masimo血氧計。

致謝

特別感謝美國威斯康辛大學生物醫學工程系的Amit Nimunkar博士,他對本文提出了很多建設性建議,並協助完成了校對工作。此外,感謝Everett Smith博士的寶貴協助和鼓勵,以及ADI為本專案提供的受試者原始資料。

參考文獻

  1. 1 Amal Jubran。“Pulse Oximetry.”。《Critical Care》,第3卷,第2期,1999年2月。
  2. 2Han-Wook Lee、Ju-Won Lee、Won-Geun Jun和Gun-Ki Lee。“The Periodic Moving Average Filter for Removing Motion Artifacts from PPG Signals.”。《國際控制自動化與系統雜誌》,第5卷,第6期,2007年12月。
  3. 3 Brendan Conlon、James A. Devine和James A. Dittmar。“ECG Synchronized Pulse Oximeter.”。美國專利4,960,126,1990年10月。
  4. 4James Reuss和Dennis Bahr。“Period 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. 5 Eric Jacobsen和Richard Lyons。“The Sliding DFT.”。《IEEE訊號處理雜誌》,第20卷,第2期,2003年3月。
  6. 6 Eric Jacobsen和Richard Lyons。“An Update to the Sliding DFT.”。《IEEE訊號處理雜誌》,第21卷,第1期,2004年1月。
  7. 7Lawrence R. Rabiner和Bernard Gold。Theory and Application of Digital Signal Processing。Prentice-Hall,1975年1月。
  8. 8 Ludvik Alkhoury、Ji-Won Choi、Chizhong Wang、Arjun Rajasekar、Sayandeep Acharya、Sean Mahoney、Barry S.Shender、Leonid Hrebien和Mose Kam。 “Heart-Rate Tuned Comb Filters For Processing Photoplethysmogram (PPG) Signals in Pulse Oximetry.”。《臨床監測與計算雜誌》,第35卷,第4期,2021年8月。
  9. 9 Recommended Configurations and Operating Profiles for MAX30101/MAX30102 EV Kits.”。ADI,2018年3月。
  10. 10Sang-Soo Oak和Praveen Aroul。“How to Design Peripheral Oxygen Saturation (SpO2) and Optical Heart Rate Monitoring (OHRM) Systems Using the AFE4403.”。德州儀器,2015年3月。
  11. 11 Douglas Altman和J. Martin Bland。“Measurement in Medicine: The Analysis of Method Comparison Studies.”。《皇家統計學會志,D輯:統計學家》,第32卷,第3期,1983年9月。
  12. 12 J. Martin Bland和Douglas G. Altman。“Statistical Methods for Assessing Agreement Between Two Methods of Clinical Measurement.”。《柳葉刀》,1986年2月。
  13. 13 Julian M. Goldman、Michael T. Petterson、Robert J. Kopotic和Steven J. Barker。“Masimo Signal Extraction Pulse Oximetry.”。《臨床監測與計算雜誌》,第16卷,第7期,2000年。
  14. 14Sagi Levental、Elie Picard、Francis Mimouni、Leon Joseph、Tal Y Samuel、Reuben Bromiker、Dror Mandel、Nissim Arish和Shmuel Goldberg。“Sex-Linked Difference in Blood Oxygen Saturation.”。《臨床呼吸雜誌》,第12卷,第5期,2018年5月。
  15. 15 Sam Koblenski。 “Everyday DSP for Programmers: DC and Impulsive Noise Removal.”。2015年11月。
  16. 16 Surekha Palreddy。“Signal Processing Algorithms.”。Design of Pulse Oximeters,第一版,1997年10月。
  17. 17 Terry L. Rusch、Ravi Sankar和John E. Scharf。“Signal Processing Methods for Pulse Oximetry.”。《生物學和醫學中的電腦雜誌》,第26卷,第2期,1996年3月。