我正在尝试估计心电信号心率变异性的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方法工作.