python信号分析_Python频谱分析

我正在尝试估计ECG信号的心率变异性的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(比分析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

因此，我想知道：

有人可以建议我应该阅读的文件，以提高我对光谱分析的理解吗？

我的方法出了什么问题？

如何为welch功能选择最合适的参数？

虽然这两个图以某种方式具有相同的形状，但数据完全不同。我该如何改进呢？

有没有更好的方法来解决这个问题？我正在考虑使用Lomb-Scargle估计，但我等着至少让Welch方法起作用。