Python Spectrum Analysis

Question

I am trying to estimate the PSD of the heart rate variability of an ECG signal. To test my code,I have extracted the R-R interval from from the fantasia ECG database. I have extracted the signal can be accessed here. To calculate the PSD, I am using the welch method as shown below:

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)

Next,the area under the curve is calculated to estimate the power spectrum of the different HRV bands as shown below

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)

I then plot the PSD and get the following plot

At first I tough the output looks okay. However, when I compare my output with that of Kubios, which is a software than analyze HRV, I noticed that there are differences. The following chart shows the expected value for the PSD as calculated by Kubios Namely, the two plots are visually different and their values are way different. To confirm this, a print out of my data clearly shows that my calculation are wrong

    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

I am thus, wondering:

• Can someone suggest a document I should read to improve my understanding on spectra analysis?
• What's wrong with my approach?
• How do I choose the most suitable parameters for the welch function?
• While the two plots somehow have the same shape, the data is completely different. How can I improve this?
• Is there a better approach to solve this? I am thinking about using the the Lomb-Scargle estimate but I am waiting to get at least the Welch method to work.

Answer 1

The problem here is that you do not handle correctly the sampeling of your signal. In your welsch call, you consider a regularly sampled signal with sample frequency 4Hz. If you look at the time vector 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)

Your signal is thus not regularly sampled. You thus need to take that into acount, else you do not estimate correctly the PSD and this gives you bad estimates.

A first correction should be to use correctly the parameter fs in welsch. This parameter indicates the sampling frequency of the given signal. Putting it ot 4 means that your time vector should be a regular [0, .25, .5, .75, .1, ....]. A better estimate would be either the median of dt or len(t)/(t.max()-t.min()), so around 4/3. This gives better PSD estimate and correct order for some of the constant but it is still different compared to Kubios values.

To get correct estimate of the PSD, you should use a non uniform DFT. A package that implement such transform can be found here. The documentation is quite cryptic for this package but you need to use the adjoint method to get the Fourier Transform without scaling issue:

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))

Here the estimates do not seems to be in line with the one from Kubios. It is possible the estimation is probably of because you do a psd estimate on the whole signal. You can try to use the welch technique combined with this nfft by averaging estimates on windowed signals as it do not rely on FFT but on any estimation of the PSD.