PSD estimation of a RRi serie with Python

The new release o Scipy comes with the  modified periodogram method called Welch’s Periodogram. It is very similar to the pwelch function provided by Matlab and you’ll find it in scipy.signal.

Because some indices of heart rate variability are calculated in the frequency domain we first need to estimate the power spectral density (PSD) of the respective RRi signal. The welch function from scipy.signal allows you to choose the size of the segments, choose the overlap size and the window function (hanning, hamming,kaiser etc). It also provides two types of detrend methods: constant and linear. Although the welch function from scipy module is quite similar to the pwelch from Matlab, it does not calculates the Confidence Interval, which is a very important parameter, once RRi signals are stochastic.
So, in this post I’ll show how to estimate the PSD of a given RRi signal, calculate the 95% Confidence Interval of the respective PSD and then calculate its area under the curve: Very Low Frequency (VLF), Low Frequency (LF) and High Frequency (HF). The code will be written in Python with some parallel comparison to the Matlab pwelch function.

To implement this script we need just some functions from numpy and scipy:

from  numpy import cumsum, interp, arange
from scipy.signal import welch
from scipy.interpolate import splrep, splev

Then we need to load a RRi signal. In this example we’ll open a text file with the RRi values in a single column.

#Open RRi file
rri = [int(rri.strip())
       for rri in open('filename.txt') if rri.strip()]

From the RRi values we need to create a time array. If your RRi values are in milliseconds you have to divide it by 1000.0, because later we’ll interpolate the RRi signal in 4 Hz.

#Create time array
# Divide by 1000.0 if RRi are in miliseconds.
t = cumsum(rri) / 1000.0 
#Force the time array to start from zero.
t -= t[0]

Because RRi are unevenly spaced we must interpolate it for accurate PSD estimation. In this case we’ll use the cubic spline interpolation.

#Create evenly spaced time array (4Hz).

#After the creation of a evenly spaced time
#array we are able to interpolate the RRi signal.
tx = arange(t[0], t[-1], 1.0 / 4.0)
tck = splrep(t, rri, s=0)
rrix = splev(tx, tck, der=0)

And Finally the PSD estimation with the non parametric Welch’s Method. As you can see here, the function provided by scipy for welch periodogam estimation look very similar to Matlab.

The Welch periodogram is an adaptation of the Bartlett method, which uses windows and data segmentation with overlap. The idea of the windows is to reduce the spectral leakage resulting in a decrease in spectral bias. With the overlaps more estimatives are made, so when all estimatives are averaged the spectral variance decreases.

Let’s consider a signal with N samples . The N samples signal is divided in P segments of D samples each with a shift S between adjacents segments. This way, the number of estimates are P, which is the integer part of (N-D)/S + 1.

The Welch periodogram can be calculated using the following equations.

where is the discrete Fourier Transform of the pth segment:

and U is the discrete window energy:

using a hanning window with 256 points its discrete energy is ~= 0.373.

Finally, averaging all P estimates we have the Welch’s periodogram:

#PSD esitmation with hanning window, 256 points each segment with 50% overlap.
Fxx, Pxx = welch(x=rrix, fs=4.0, window="hanning",
                 nperseg=256, noverlap=128,
                 detrend="linear")

When we use pwelch in Matlab it’s possible to give two other parameters to this function: a string ‘ConfidenceLevel’ and the probability of the confidence interval. Giving these extra arguments a thirth output is returned:

[Pxx, Fxx, Pxxc] = pwelch(signal, window(segment),...
                           overlap, nfft, fs)

The Pxxc contains the 95% probability confidence interval of the PSD.

The welch function provided in scipy do not support theses argurments to calculate the confidence interval, but we can do it with few more lines of code using the chi2 function from scipy.stats

We can calculate the confidence interval of the actual PSD using a chi-square relationship. Given the estimate, the actual PSD lies between the limits with a desired probability. These relationship can described as follows:

With these relation we can say that the confidence interval for the true Pxx is then:

where v is the degree of freedom and its equal to P (number of estimates) and p is the desired significance (p=0.05 for the 95% confidence interval). So to calculate the confidence interval of the PSD we’ll use the Percent Point Function of the chi-square with p = 0.05.

from scipy.stats import chi2

probability = 0.95

#P is the number of estimates in welch function
#and also the degree of freedom.
alfa = 1 - probability
v = 2 * P
c = chi2.ppf([1 - alfa / 2, alfa / 2], v)
c = v / c
Pxxc_lower = Pxx * c[0]
Pxxc_upper = Pxx * c[1]

With our PSD we can now calculate the area under the curve from 0-0.04 Hz to estimate VLF, 0.04-0.15 to LF and 0.15-0.4 for HF. We’ll create a function called psd_auc that receives Fxx, Pxx and the frequency interval of each paramter. This function calculates the AUC with the sum of the respective PSD values multiplied by de frequency resolution:

from numpy import logical_and

def psdauc(Fxx, Pxx, vlf=0.04, lf = 0.15, hf = 0.4):
    df = Fxx[1] - Fxx[0]
    vlf = sum(Pxx(Fxx <=  vlf)) * df
    lf = sum(Pxx[logical_and(Fxx >= vlf, Fxx < lf)]) * df
    hf = sum(Pxx[logical_and(Fxx >= lf, Fxx < hf)]) * df
    return vlf, lf, hf

You can find the code of this post here