我正在尝试估计心电信号心率变异性的PSD.为了测试我的代码,我从
fantasia ECG database提取了R-R间隔.我提取了可以访问的信号
here.要计算PSD,我使用的是如下所示的welch方法:
import matplotlib.pyplot as plt
import numpy as np
from scipy.signal import welch
ibi_signal = np.loadtxt('fantasia-f1y01-RR.txt')
t = np.array(ibi_signal[:, 0]) # time index in seconds
ibi = np.array(ibi_signal[:, 1]) # the IBI in seconds
# Convert the IBI in milliseconds
ibi = ibi * 1000
# Calculate the welch estimate
Fxx, Pxx = welch(ibi, fs=4.0, window='hanning', nperseg=256, noverlap=128)
接下来,计算曲线下面积以估计不同HRV频带的功率谱,如下所示
ulf = 0.003
vlf = 0.04
lf = 0.15
hf = 0.4
Fs = 250
# find the indexes corresponding to the VLF, LF, and HF bands
ulf_freq_band = (Fxx <= ulf)
vlf_freq_band = (Fxx >= ulf) & (Fxx <= vlf)
lf_freq_band = (Fxx >= vlf) & (Fxx <= lf)
hf_freq_band = (Fxx >= lf) & (Fxx <= hf)
tp_freq_band = (Fxx >= 0) & (Fxx <= hf)
# Calculate the area under the given frequency band
dy = 1.0 / Fs
ULF = np.trapz(y=abs(Pxx[ulf_freq_band]), x=None, dx=dy)
VLF = np.trapz(y=abs(Pxx[vlf_freq_band]), x=None, dx=dy)
LF = np.trapz(y=abs(Pxx[lf_freq_band]), x=None, dx=dy)
HF = np.trapz(y=abs(Pxx[hf_freq_band]), x=None, dx=dy)
TP = np.trapz(y=abs(Pxx[tp_freq_band]), x=None, dx=dy)
LF_HF = float(LF) / HF
HF_LF = float(HF) / LF
HF_NU = float(HF) / (TP - VLF)
LF_NU = float(LF) / (TP - VLF)
然后我绘制PSD并获得以下图
起初我坚韧的输出看起来不错.然而,当我将我的输出与Kubios的输出进行比较时,Kubios是比分析HRV的软件,我注意到有差异.下图显示了Kubios计算出的PSD的预期值
名称上,两个图在视觉上不同,它们的值是不同的.为了确认这一点,我的数据中的打印清楚地表明我的计算是错误的
ULF 0.0
VLF 13.7412277853
LF 45.3602063444
HF 147.371442221
TP 239.521363002
LF_HF 0.307795090152
HF_LF 3.2489147228
HF_NU 0.652721029154
LF_NU 0.200904328012
因此,我想知道:
有人可以提出我应该阅读的文件来提高我对光谱分析的理解吗?
>我的做法有什么问题?
>如何选择最适合的焊接功能参数?
>虽然这两个地块的形状不一样,但数据完全不一样.如何改善这个?
有没有更好的方法来解决这个问题?我正在考虑使用Lomb-Scargle的估计,但我正在等待至少让Welch方法工作.