Python-hrvanalysis庫 使用總結

最近接手一個心電專案，用到了hrvanalysis庫分析心電資料特徵。網上的關於hrvanalysis中文資源不多，今天總結一下。

參考 https://github.com/Aura-healthcare/hrvanalysis

1 hrvanalysis介紹和安裝

hrvanalysis是一個基於SciPy, AstroPy, Nolds and NumPy的python模組，用於構建RR-interval心率變異分析的工具。

python3下，使用下面?命令安裝這個庫

$ pip install hrv-analysis

hrvanalysis對其他庫的依賴如下：

Python (>= 3.5)
astropy >= 3.0.4
future >= 0.16.0
nolds >= 0.4.1
numpy >= 1.15.1
scipy >= 1.1.0

2 離群值和異常值刪除

# 匯入庫和對應的函式
from hrvanalysis import remove_outliers
from hrvanalysis import remove_ectopic_beats
from hrvanalysis import interpolate_nan_values

# 此為RR-interval長度列表，每個元素為每個RR-interval的間隔時間
rr_intervals_list = [1000, 1050, 1020, 1080, ..., 1100, 1110, 1060]
# 刪除離群值，被刪除的元素置為nan，自low_rri和high__rri分別為最小和為最大的RR-interval
rr_intervals_without_outliers = remove_outliers(rr_intervals=rr_intervals_list,  
                                                low_rri=300, high_rri=2000)
# 將離群的rr-interval值刪除，並用線性插值的方法新增新值
interpolated_rr_intervals = interpolate_nan_values(rr_intervals=rr_intervals_without_outliers,
                                                   interpolation_method="linear")
# 從RR-interval心電訊號中，採用malik方法刪除異常值，即將值置為nan
nn_intervals_list = remove_ectopic_beats(rr_intervals=interpolated_rr_intervals, method="malik")
# 使用線性插值的方法替換nan
interpolated_nn_intervals = interpolate_nan_values(rr_intervals=nn_intervals_list)

通過這三個函式可以清洗RR-interval列表的資料，降低離群值和異常值的出現，利於後續分析。

3 RR-interval特徵提取

寫在前面 RR-interval和nn-interval本質上是一樣的。

get_time_domain_features()

#返回包含HRV分析的時域特徵字典，通常用於24h的長期記錄，也有些研究在2-5min的短期記錄
from hrvanalysis import get_time_domain_features

nn_intervals_list = [1000, 1050, 1020, 1080, ..., 1100, 1110, 1060]
time_domain_features = get_time_domain_features(nn_intervals_list)
# 返回time_domain_features是一個字典，包含如下鍵-值：
mean_nni     RR間期平均值
sdnn         RR間期標準差
sdsd         相鄰RR之間差異的標準差
rmssd        相鄰RR間期的差方和的均值的平方根，反映高頻(快速或副交感神經)對hrV的影響
median_nni   RR間隔差的中位數的絕對值
nni_50       連續RR間期間隔差大於50ms的數量
pnni_50      用nni_50 / RR間期數量
nni_20       連續RR間期間隔差大於20ms的數量
pnni_20      用nni_20 / RR間期數量
range_nni    RR間期最大值和最小值的差異
cvcd         連續差異的變化係數 = rmssd / mean_nni
cvnni        變異係數 = sdnn / mean_nni
mean_hr      平均心率
max_hr       最大心率
min_hr       最小心率
std_hr       心率標準差

get_geometrical_features()

#返回返回幾何時域特徵字典，記錄的資料必須在20min-24h內
from hrvanalysis import get_geometrical_features

nn_intervals_list = [1000, 1050, 1020, 1080, ..., 1100, 1110, 1060]
geometrical_features = get_geometrical_features(nn_intervals_list)
#返回geometrical_features字典，包含如下鍵值：
triangular_index    HRV三角指數測量值 = 密度分佈的積分 = RR間期個數 / 密度分佈最大值
tinn        		測量以三角形為底的分佈的基線寬度

get_frequency_domain_features()

# 返回包含用於HRV分析的頻域特徵的字典，必須在2-5min的視窗內使用
from hrvanalysis import get_frequency_domain_features

nn_intervals_list = [1000, 1050, 1020, 1080, ..., 1100, 1110, 1060]
frequency_domain_features = get_frequency_domain_features(nn_intervals_list)
#返回字典，包含如下鍵值：
total_power   總功率密度譜
vlf           極低頻時（0.003-0.04hz）HRV的方差，反映心臟產生的內在節律，受交感神經活動調節
lf            低頻時（0.04-0.15hz）HRV的方差，反映交感神經和副交感神經活動
hf            高頻時（0.15-40hz）HRV的方差，反映了由副交感神經(迷走神經)活動引起的一拍一拍變異性的快速變化。
lf_hf_ratio   等於 lf / hf
lfnu          歸一化低頻功率
hfnu          歸一化高頻功率

get_csi_cvi_features()

# 返回來自非線性域的3個用於HRV分析的特徵字典，必須在30、50、100秒的RR-interval中使用
from hrvanalysis import get_csi_cvi_features

nn_intervals_list = [1000, 1050, 1020, 1080, ..., 1100, 1110, 1060]
csi_cvi_features = get_csi_cvi_features(nn_intervals_list)
#返回字典包含如下鍵值：
csi        交感神經指數
cvi        迷走神經指數
Modified_csi        改良的CSI 是癲癇發作檢測研究中的替代方法

get_poincare_plot_features()

# 返回包含非線性域的3個特徵字典，用於HRV分析，必須在5分鐘的短期視窗使用
from hrvanalysis import get_poincare_plot_features

nn_intervals_list = [1000, 1050, 1020, 1080, ..., 1100, 1110, 1060]
poincare_plot_features = get_poincare_plot_features(nn_intervals_list)
#返回如下鍵值：
Sd1      龐加萊圖在垂直於等式的線的直線上投影的標準差
Sd2      龐加萊投影的標準差
ratio_sd2_sd1        sd2 / sd1

get_sampen()

# 計算給定資料的樣本熵，必須在一分鐘內的短期視窗使用
from hrvanalysis import get_sampen

nn_intervals_list = [1000, 1050, 1020, 1080, ..., 1100, 1110, 1060]
sampen = get_sampen(nn_intervals_list)
# 返回如下鍵值：
sampen    資料的樣本熵

4 繪圖函式

plot_psd()

from hrvanalysis import plot_psd

nn_intervals_list = [1000, 1050, 1020, 1080, ..., 1100, 1110, 1060]
# 採用兩種方法繪圖
plot_psd(nn_intervals_list, method="welch")
plot_psd(nn_intervals_list, method="lomb")

這裡附上官方原圖，VLF、LF、HF三個指標在get_frequency_domain_features()中均有輸出。

get_frequency_domain_features()

from hrvanalysis import plot_poincare

nn_intervals_list = [1000, 1050, 1020, 1080, ..., 1100, 1110, 1060]
plot_poincare(nn_intervals_list)
plot_poincare(nn_intervals_list, plot_sd_features=True)

還是官方例圖，SD1和SD2指標為get_poincare_plot_features()的輸出。

總結一下

我的心搏片段切分的比較小，與分佈特徵提取函式建議的最佳輸入片段長度要求不符合，導致使用時有些特徵提取失敗，產生空值。因此在使用hrvanalysis提取心搏特徵時，需要根據自己資料的特點，選用合適的特徵提取函式。