2023年2月9日

脈搏感應器PPG的應用(6)：心率變異性的頻域分析與pythnon程式影片說明

心率變異性做頻域的分析，其中一個目的就是在評估自律神經的功能，而頻域分析最重要的工具就是使用FFT。

而使用FFT做頻域分析有幾個限制：

(1)數據點需要相同取樣頻率。

假設每秒取樣一次，代表要知道該秒時，心跳的BPM是多少，但實際上不可能知道每秒時的BPM，因為你那秒不一定有心跳啊。所以這部分要進行interpolate插值，然後再重新取樣(resampling)

那麼應該重新取樣成多少數據點呢？牽涉到FFT的下一個限制

(2)在FFT的演算法中，需要用二的冪次個數據點來運算才會最有效率，如果超過要刪除，不夠的話就補零。(二的冪次如2、4、8、...)

我看了幾篇論文都是用 7.11 Hz來重新取樣(google 7.11Hz HRV)，我推測原因是這樣：每個量測時間段約5分鐘(300秒)，取288 秒，用7.11Hz重新取樣，則可以得到 7.11 x 288 = 2047.68個數據點，至少有 2048個數據點可以做FFT。使用 7.11 Hz的取樣頻率，則Nyquist frequency為 7.11/2 =  3.55 Hz，能表示的頻率範圍為 0 ~ 3.55 Hz。
數據長度為2048，而FFT之後信號在頻域的表現是對稱的，所以只取一半。總頻帶是 0 ~ 3.55Hz，中間分布 2048/2=1024個點，所以每個bins的頻率是 3.55/1024 = 0.00346679687 Hz

在這個指引中，HRV的頻域分析上，會分成四個頻段來分析，而在這篇【An Overview of Heart Rate Variability Metrics and Norms】說明了不同頻段的研究限制和意義。

以下是我的某次測量的分析結果圖，分成了三個區塊。由左至右為VLF、LF、HF




ULF 頻段 (ultra low frequency，≤0.003 Hz) 
  • 週期為 5 分鐘至 24 小時，要用 24 小時記錄來測量，無法用5分鐘的短期測量來觀察。影響這個頻段的可能包含晝夜節律、核心體溫變化、新陳代謝和腎素-血管緊張素系統。
VLF 頻段  (very low frequency，0.0033–0.04 Hz) 
  • 周期在 25 到 300 秒之間。需要至少 5 分鐘的記錄時間，但最好在 24 小時內進行監控。此頻段可能受到身體活動 、體溫調節、腎素-血管緊張素和內皮對心臟的影響，另外副交感神經的活動可能有助於 VLF 功率。
LF  頻段(low frequency，0.04–0.15 Hz) 
  • 代表自律神經(交感+副交感)的活性
  • 周期為 7 至 25 秒的節奏組成，會受每分鐘 ~3 至 9 次的呼吸頻率影響 。在 5 分鐘內有 12-45 個完整的振盪週期 。此頻段主要反映休息狀態的壓力感受器活動。此頻段的功率可能由交感和副交感產生。
  • 血壓調節的方式是藉由感壓受器，主要由副交感神經或單獨通過壓力反射而活動。
  • 交感神經似乎不會產生遠高於 0.1 Hz的頻率，而副交感神經會影響低至 0.05  Hz（20 秒節律）的心律。
  • 在休息狀態下，LF 頻段反映出壓力反射活動不是受交感神經支配。
HF 頻段 (high frequency，0.15–0.40 Hz) 
  • 代表副交感的活性
  • 受 9 至 24 bpm 的呼吸影響 
低頻與高頻功率之比（LF/HF ratio）
  • 可以估計交感神經系統和副交感神經系統活動之間的比率。
  • 交感和副交感的活動都對 LF 功率有貢獻，而副交感則主要對 HF 功率有貢獻。
  • 有些研究會將此比當作交感與副交感的活動比例，但因為LF的功率同時包含交感和副交感，所以這點有爭議。

這篇也同時提到幾點
  • 副交感神經發揮作用比交感神經（>5 秒）更快（<1 秒）
  • 短期測量的HRV受到以下幾個因素影響：1.交感神經和副交感神經的交互影響，2.RSA(呼吸性竇性心律不整)控制心跳的調節機制、壓力感受器反射（血壓的負反饋控制）和血管張力的節律變化 。
  • 交感神經可以抑制副交感神經活動，但也可能增加副交感神經的反應。
  • 交感和副交感的關係並不像翹翹板那樣一個高一個就低。副交感的活動增加可能與交感的活動減少、增加或無變化有關
程式分析的方式
最後終於講到分析方式

在Kaggle上，有一篇是利用python對 ECG進行分析的文章

我參考這篇文章和前述的文獻，改寫了程式碼。
整個分析程序分以下幾個部分
  1. micro:bit 的程式只負責讀取類比pin，然後印出。要注意指定傳輸速度為 115200


  2. 先用這個python的程式去讀取micro:bit的資料，並加上以微秒為單位的時間戳記，並存檔。此程式也同時展示脈搏波形PPG與BPM變化


  3. 使用另外一個python程式(用jupter notebook)，先讀取方才存檔的文字檔後，前段程式負責劃出段測量時間的脈搏圖，使用者須自行看圖，藉此用來決定用哪一段時間來進行分析。因為通常開頭會有雜訊，而結尾會有一些資料無法讀取。
  4. 偵測峰值之後，會在脈搏的PPG上標示峰值點，可藉此判斷是否都標示出正確峰值。如果這部分有異常，將影響後面的時域分析和頻域分析。


  5. 確認峰值沒錯之後，直接執行後半的程式，即可產出分析的結果報告。圖上同時會有時域分析、頻域分析、非線性分析、幾何分析等結果。
分析結果圖片

HRV變化圖
灰色水平線為平均HR，兩條紅線為標準差。多數都位於一個標準差內

龐加萊圖Poincaré plot
也稱為Lorenz (Poincaré) plots (LPs)，是將此次心跳間隔和下次心跳間隔畫在圖上。假設每次心跳間隔都相同，則數據點都會落在同一點上，但因為有心跳變異性，所以健康人應該會是分散的。

我在搜尋這圖片資料時，經常看到一張標示出處是台北榮總(Taipei Veterans General Hospital, Taiwan)的圖片，是觀察不同年齡和身體狀況的病人的心跳間隔。新生兒位於左下角，而右上角是加護病房的病人。




心跳間隔的直方圖
這是以50 ms為間隔，藉由這張圖可以分析壓力，例如 baevsky stress index這個壓力指數。處於壓力之下的人，交感神經活性強，心率變異會降低，所以直方圖會變窄而高，大多數心跳間隔位在一個小小區間內。而壓力低的人，有比較高的心率變異，所以直方圖會矮而寬。
方法是：

  • 分子為發生頻率最高區間次數佔全部心率次數百分比，例如下圖高度約為169，全部心跳數為296，分子為 57%
  • Mo為發生頻率最高區間的心率值，單位為秒。此圖為750-800之間，取0.775秒
  • var 心率分佈區間，最高心率與最低心率區間相差，單位為秒。此圖為(900-650)/1000=0.25
  • SI即為壓力指數147



我在國內一家作心率感測的舒壓儀廠商的說明書看到的壓力指數與壓力狀態的對應關係如下
  • 小於 50 紀錄異常導致訊號讀取不正確或長期接受紓壓訓練
  • 50-150 為一般正常狀態
  • 150-500 心理壓力偏大,或是過度運動導致的疲勞狀態
  • 500-900 表示嚴重心理壓力，長時間處於高壓狀態，或是心臟跳動有異常
  • >900 表示高度心理壓力或是存在其他生理疾病，需尋求醫師協助
其參考的文獻如下
  • Roman М. Baevsky, A.G.C., Heart rate variability analysis: physiological foundations and main methods. Cardiometry, 2017. 10: p. 66-76.

  • Bayevsky, R.M., et al., HRV Analysis under the usage of different electrocardiography systems ((Methodical recommendations) The Committee of Clinic Diagnostic Apparatus and The Committee of New Medical Techniques of Ministry of Health of Russia, protocol, 2002. 4.



以下是整個時域分析和頻域分析的程式說明和結果影片


【用micro:bit與脈搏感應器偵測脈搏。再以python做時域和頻域分析HRV，用以分析自律神經活性與身體狀況】

程式位置


最後可以產生的報告如下圖


-
Labels: ,