Python频谱分析

Question

我试图估计ECG信号的心率变异性的PSD.为了测试我的代码,我从幻想曲心电图数据库中提取了RR间隔.我提取的信号可以在这里访问.要计算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,这会给您带来不好的估计.

第一个修正应该是正确使用fswelsch中的参数.该参数表示给定信号的采样频率.把它放到4意味着你的时间向量应该是常规的[0, .25, .5, .75, .1, ....].一个更好的估计可能是中位数,dt或者len(t)/(t.max()-t.min())左右4/3.这给出了更好的PSD估计和一些常数的正确顺序,但与Kubios值相比仍然不同.

要获得PSD的正确估计,您应该使用非均匀DFT.可以在此处找到实现此类转换的包.该文档对于这个包非常神秘,但您需要使用伴随方法来获得傅立叶变换而不会出现缩放问题:

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技术结合此nfft,通过对窗口信号进行平均估计,因为它不依赖于FFT,而是依赖于PSD的任何估计.