rPPG訊號(脈搏波)獲取測試程式_透過臉部影片額頭部位亮度變化_來擷取生理訊號以測量心率

 
脈搏波是隨心臟跳動而變化的血流量。回憶過去若曾去過醫院做健檢,你應該有看過類似感測器夾在手指上,測量出像山峰一般的波形,那就是脈搏波。脈搏波是利用血液中的血紅素會大量吸收綠光的特性。
市面上很多運動手錶,如今也已涵蓋心律、血氧測量等諸多身體指數量測,多半採用到所謂脈搏感測技術,常見的rPPG藉由反射行脈波測量手段。


用綠光的原因在於能被淺層靜脈中的血紅素部分吸收,剩餘沒被吸收的綠光則會反射回來,有利於被手錶感測器量測。
至於紅光部分則穿透能力較強,能反射回來被感測器分析的比例較少,而不利於運動時全程心率量測。

運動手錶主要就是利用血紅素密度變化,去間接推斷人體心率變化。每當人體心臟收縮跟舒張,血管中的血紅素密度就會有密度起伏。若經過手錶等感測裝置的血紅素密度高,被偵測到的綠光就少。反之,若經過的血紅素密度低,則被偵測到的綠光就多。



而從視訊或是影片中測量脈搏波的原理也是一樣的。一般電腦中視訊、影片有分為RGB三種顏色通道不同的光線強度,以 0 到 255 個階度的亮度值進行記錄。
我們可以拍攝或擷取皮膚表面ROI所得視訊影像中的綠色成分光線強度,根據心臟跳動引起的血紅素變化量而產生變化。
同理可得知,拍攝皮膚表面的視訊影像中,綠色成分的變動即為脈搏波。

為了讓血紅素增加時訊號隨之變大,脈搏波訊號會使用從最大亮度值(255)來減去綠色亮度值後所得的數值。透過計算一分鐘內出現多少次脈搏波,也就是測量一分鐘內心臟跳動的次數,就能從視訊影像中測量出心率。

拿一個10秒的公開測試影片作測試
Stage1.先定位額頭部位ROI

import cv2
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
movie = "face.mp4"
video_path = ("C:/img/{}".format(movie))
cap = cv2.VideoCapture(video_path)
cap.set(cv2.CAP_PROP_POS_FRAMES, 0)
roi = (147, 69, 340, 178)
ret, frame_bgr = cap.read() #先只讀取第一frame
if ret:
    #OpenCV 讀取的影像色彩順序是 BGR,但 matplotlib.pyplot.imshow() 預期的是 RGB
    # OpenCV 的 BGR 轉為 Matplotlib 使用的 RGB
    frame_rgb = cv2.cvtColor(frame_bgr, cv2.COLOR_BGR2RGB)
    # 取得 ROI 座標
    x1, y1, x2, y2 = roi
    # 裁切 ROI
    roi_rgb = frame_rgb[y1:y2, x1:x2]
    # 複製影像,避免直接修改原始 frame_rgb
    frame_with_roi = frame_rgb.copy()
    # 在完整影像上畫出 ROI 方框
    cv2.rectangle(frame_with_roi,(x1, y1),(x2, y2),(255, 0, 0),2)

    # 建立並排顯示畫面
    fig, axes = plt.subplots(1, 2, figsize=(12, 6))
    # 顯示完整影像
    axes[0].imshow(frame_with_roi)
    axes[0].set_title("Full Frame with ROI")
    axes[0].axis("off")

    # 顯示 ROI
    axes[1].imshow(roi_rgb)
    axes[1].set_title("Selected ROI")
    axes[1].axis("off")

    plt.tight_layout()
    plt.show()



cap.release()
cv2.destroyAllWindows()

Stage2.形成 300 個時間序列數值=>尋找相鄰脈搏低谷,低谷間隔換算 BPM。
影片長度為 10 秒,幀率為 30 FPS
實際心率可能包括:
  • 靜止或運動員:低於 60 BPM
  • 一般靜止:60~100 BPM
  • 運動、緊張:高於 100 BPM

FIR濾波器係數實測調整到161會有比較好的波形訊號
階段2程式

import cv2
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

#Phase 1:確認 ROI 擷取位置
movie = "face.mp4"
video_path = ("C:/img/{}".format(movie))
cap = cv2.VideoCapture(video_path)
cap.set(cv2.CAP_PROP_POS_FRAMES, 0)#將影片移到第 0 幀
roi = (147, 69, 340, 178)
ret, frame_bgr = cap.read()
if ret:
    frame_rgb = cv2.cvtColor(frame_bgr, cv2.COLOR_BGR2RGB)
    select_roi = frame_rgb[roi[1]: roi[3], roi[0]: roi[2]]
    plt.imshow(select_roi)
    plt.show()
cap.release()
cv2.destroyAllWindows()


#Phase 2:逐幀擷取綠色訊號
cap = cv2.VideoCapture(video_path)
g_list = []
while (cap.isOpened()):
    ret, frame = cap.read()
    if ret == False:
        break
    # 計測額頭ROI區域
    roi = (147, 69, 340, 178)
    roi = frame[roi[1]: roi[3], roi[0]: roi[2]]
    roi = cv2.medianBlur(roi, 5)#每個像素都參考周圍 5×5 區域的像素,把中間值當作新結果。,降低雜訊。
    #通道分割
    b, g, r = cv2.split(roi)
    #g_list.append(g.mean())#綠色平均亮度值
    g_list.append(255-g.mean())#將整個額頭 ROI 中所有綠色像素取平均,然後做亮度反轉。(原本的波峰 ↔ 反轉後的波谷)
cap.release()


g_array = np.array(g_list)#將list轉換為 ndarray
#FIR濾波器變數
sample_rate = 30#取樣率:每秒取得 30 個綠色亮度樣本。(寫比較動態一點是cap.get(cv2.CAP_PROP_FPS))
dt = 1 / 30 #如果兩個脈搏低谷相差 20 幀,就可以去用 20* (1/30)換算經過的秒數0.6667秒,心率去用60/0.6667去算90BPM
nyq = sample_rate / 2 #Nyquist 頻率:最高只能正確辨認到 15 Hz 的變化。至少需要兩個樣本才能描述一個週期。
#正常的心率速度通常是60-100 次/分
#使用帶通濾波器(Band-pass filter),僅允許頻率介於 0.7Hz 至 3Hz 之間的訊號通過。
hf = 3 #最高保留頻率3*60=180BPM
lf = 0.7 #最低保留頻率0.7*60=42BPM
#截止頻率(正規化至 0 到 1 之間)
wh = hf / nyq
wl = lf / nyq
numtaps = 161 #濾波器係數的數量(FIR 濾波器階數)
#FIR濾波器響應延遲
delay = (numtaps - 1) / 2 * dt #(161 - 1) / 2 / 30 = 2.67 秒

fir = signal.firwin(numtaps, cutoff=[wl, wh], window="hann", pass_zero=False)
gf = signal.lfilter(fir, 1, g_array)


#為了計算心跳數,檢測脈搏的低谷,然後計算脈搏與脈搏之間的間隔,即心跳間期。
minimal_idx_gf = signal.argrelmin(gf, order=10)
bottom_number = minimal_idx_gf[0]
#bottom_number = np.delete(bottom_number, slice(0, 5))
bottom_count = len(bottom_number)

#心率變異性
i = 0
bottom_interval = []
while i < bottom_count - 1:
    mi = (bottom_number[i + 1] - bottom_number[i])
    bottom_interval.append(mi)
    i += 1

heart_rate_interval = np.array(bottom_interval) * dt
#心跳數就是每分鐘心臟跳動的次數,可以通過將 60 除以心跳間期來計算。
hri_ave = np.average(heart_rate_interval)
heart_rate = 60 / hri_ave
print(heart_rate)
print(heart_rate_interval)


x = np.linspace(1, len(gf), len(gf))
dt = 1 / 30#每兩個資料點之間的時間,每隔約 33.3 毫秒量測一次額頭綠色亮度。
t = (x * dt) - dt
tf = t - delay

plt.plot(tf, gf, label="fir", c="tomato")
plt.plot(tf[minimal_idx_gf], gf[minimal_idx_gf], 'bo', label='peak_minimal')
plt.xlabel("time(s)")
plt.ylabel("green_average_brightness_value")
plt.xlim(0, 10)#只顯示第 0~10 秒
plt.ylim(-1.0, 1.0)
plt.show()


運行效果








Ref:
高速脈搏感測IC可量測壓力和血管年齡
https://www.eettaiwan.com/20180305np21/
脉搏传感器
脈搏與脈搏感測器
【教學】為什麼心率用綠光、血氧用紅光?智慧手錶感測原理一次看懂


留言

這個網誌中的熱門文章

SAP物料主數據(Material Master Data)

何謂淨重(Net Weight)、皮重(Tare Weight)與毛重(Gross Weight)

外貿Payment Term 付款條件(方式)常見的英文縮寫與定義