Анализ спектра Python - программирование

Я пытаюсь оценить PSD вариабельности сердечного ритма сигнала ЭКГ. Чтобы протестировать мой код, я извлек интервал R-R из базы данных ЭКГЭФФЭКТ. Я извлек сигнал, доступ к которому можно получить здесь. Чтобы рассчитать 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, который является программным обеспечением, чем анализирует ВСР, я заметил, что есть различия. Следующая таблица показывает ожидаемое значение для PSD, рассчитанное Kubios А именно, эти два графика визуально отличаются друг от друга, и их значения различны. Чтобы подтвердить это, распечатка моих данных ясно показывает, что мои вычисления ошибочны.

    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.

Ответ 1

Проблема заключается в том, что вы неправильно обрабатываете sampeling вашего сигнала. В вашем вызове welsch вы рассматриваете регулярно дискретизированный сигнал с частотой дискретизации 4 Гц. Если вы посмотрите на вектор времени 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 в welsch. Этот параметр указывает частоту дискретизации данного сигнала. Положив его на 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.