I’ve gotten python routines to calculate time domain values for Heart Rate Variability (HRV). Now I’m working on doing spectrum analysis with Frequency domain.
Frequency Domain step 1. Interpolating beats into fixed time segments.
FFT (Fast Fourier Transform) and Welch are the two most common methods for getting a frequency power spectrum. Both of those require data with evenly spaced time points. Because heart beats are not evenly spaced they need to be interpolated.
Some time back I had captured links to a number of articles that have example code for frequency power spectrum hrv calculations. Looking at the different ways that they use for handling the interpolation. One of the articles shows how to compare a graph of the interpolated values with your original values so I have been doing that.
Going slowly and just doing one step at a time but am at this point just comparing interpolation methods before moving on to actually creating power spectrums.
I download my heart beat intervals into an array called “beattimes”. “beattimes” is an array of beats, in milliseconds, from my database for one session. I will leave it to you to figure out how to fill the beattimes array whether from database or importing beat data from a file.
Let’s start by graphing the data.
import numpy as np
import matplotlib.pyplot as plt
#beattimes is magically filled by a routine that gets the intervals out of the database
t = np.cumsum(beattimes) / 1000.0
t -= t[0] #alternately t = t[1:]#that created an array t of timestamps for each of your heartbeat intervals … milliseconds divided by 1000 to be seconds
plt.title(“Original Inter Beat Intervals”)
plt.plot(t, beattimes, label=”Original”, color=’blue’)
plt.legend()
plt.show()
If you haven’t done much with graphing thought you might like a very simple graph routine. You do have to fill the “beattimes” arrray with your own code but that will graph any of your sessions of heart rate data.
So now the interpolation. I’ve looked at three sources of code so far and they all interpolate with somewhat different methods.
First article I looked at.
from scipy.interpolate import splrep, splev
tx = np.arange(t[0], t[-1], 1.0 / 4.0)
tck = splrep(t, beattimes, s=0)
rrix = splev(tx, tck, der=0)
tx is the new timestamps for each of the values. In this case they are each 1/4 of a second. (Which greatly increases the number of data points, not quite sure why 1/4 and not 1/2 or even 1 second.) “rrix” is the new array of interpolated values that match the evenly spaced time stamps.
Second article I looked at did it a little differently.
from scipy.interpolate import interp1d
tx = np.linspace(t[0],t[-1],t[-1])
f = interp1d(t, beattimes, kind=’cubic’)
rrix = f(tx)
He created the evenly spaced timestamp array, tx, a little differently and did the interpolation from a different library in scipy.
Third source was the python hrv module. They did the tx timestamp array the same as the first method but used numpy for the interpolation instead of scipy.
tx = np.arange(t[0], t[-1], 1.0 / 4.0)
rrix = np.interp(tx, t, beattimes)
So now I have three methods of interpolating my heart beat interval data into evenly spaced values. How do I compare them? Graphing is one method. I can see how the interpolated values look compared to the original.
plt.title(“Original and Interpolated Signal”)
plt.plot(t, beattimes, label=”Original”, color=’blue’)
plt.plot(tx, rrix, label=”Interpolated”, color=’red’)
plt.legend()
plt.show()