Python频谱分析

Question

我正在尝试估计 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 方法起作用。

Answer 1

这里的问题是您没有正确处理信号的采样。在您的 welsch 调用中，您考虑采样频率为 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，这会给您带来错误的估计。

第一个修正应该是正确使用威尔士语中的参数 fs。该参数表示给定信号的采样频率。把它放在 4 意味着你的时间向量应该是一个常规的 [0, .25, .5, .75, .1, ....]。更好的估计是 dt 或 len(t)/(t.max()-t.min()) 的中值，所以大约 4/3。 这为一些常数提供了更好的 PSD 估计和正确的顺序，但与 Kubios 值相比仍然不同。

要获得 PSD 的正确估计值，您应该使用 non uniform DFT . 可以找到实现这种转换的包 here .该包的文档非常神秘，但您需要使用伴随方法来获得没有缩放问题的傅里叶变换:

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 的任何估计。