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的输出进行比较时，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方法起作用。

最佳答案

这里的问题是您没有正确处理信号采样。在欢迎 session 中，您将考虑采样频率为4Hz的定期采样信号。如果您查看时间向量t

In [1]: dt = t[1:]-t[:-1]

In [2]: dt.mean(), np.median(dt)
Out[2]: 0.76693059125964014, 0.75600000000000023

In [3]: dt.min(), dt.max()
Out[3]: (0.61599999999998545, 1.0880000000000081)

因此，您的信号不会定期采样。因此，您需要考虑这一点，否则您将无法正确估计PSD，这将给您带来错误的估计。

第一个更正应该是在welsch中正确使用参数fs。此参数指示给定信号的采样频率。设置为ot 4意味着您的时间向量应为常规[0, .25, .5, .75, .1, ....]。更好的估计是dt或len(t)/(t.max()-t.min())的中位数，所以在4/3左右。
这可以为某些常数提供更好的PSD估计和正确的阶数，但与Kubios值相比仍然有所不同。

为了获得正确的PSD估算值，您应该使用non uniform DFT。
可以找到一个实现这种转换的包here。该软件包的文档非常晦涩难懂，但是您需要使用adjoint方法来获得傅立叶变换而不会出现缩放问题:

N = 128    # Number of frequency you will get
M = len(t) # Number of irregular samples you have
plan = NFFT(N, M)

# Give the sample times and precompute variable for
# the NFFT algorithm.
plan.x = t
plan.precompute()

# put your signal in `plan.f` and use the `plan.adjoint`
# to compute the Fourier transform of your signal
plan.f = ibi
Fxx = plan.adjoint()
plt.plot(abs(Fxx))

这里的估计似乎与Kubios的估计不一致。估计可能是因为您对整个信号进行了psd估计。您可以尝试将welch technique与此nfft结合使用，方法是对窗口信号的估计值求平均，因为它不依赖于FFT，而是依赖于PSD的任何估计。

关于Python频谱分析，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/40086784/