This is a Python module/tutorial for performing heart rate variability (HRV) analysis on electrocardiogram (ECG) time series.
This code is licensed under the MIT License - see the LICENSE.md file for details.
Implements the popular QRS complex detection algorithm introduced in Pan, et al (1985). The algorithm uses filtering, adaptive thresholding, and criteria based on human cardiac physiology to detect QRS complexes in the face of noise and quickly changing and diverse ECG morphologies. This function implements the Pan-Tompkins algorithm as it was originally published, along with two modifications which include additional filtering and eliminating potential QRS detections that occur within the refractory period. Since this algorithm is often used to find R-peak locations (and not just general QRS detection) for applications such as Heart Rate Variability (HRV) analysis, this function also performs a neighborhood search around the final QRS detection locations to find exact R-peak locations.
R_peak_locs = panTompkins(ECG, fs, plot = 1)
Time domain features are perhaps the simplist method computationally for performning HRV analysis. Following R-peak detection (see PanTompkins), ectopic beats are removed to produce a normal-to-normal (NN) interval series. Various statistical measures can then be extracted from the time series that provide insights into heart rate linear dynamics.
timeDomainFeats = timDomain(RR_interval_series)| Label | Description |
|---|---|
| ANN | Average NN interval |
| SDNN | Standard deviation of NN intervals |
| SDSD | Standard deviation of successive NN intervals |
| NN50 | Number of successive NN intervals differing by more than 50 ms |
| pNN50 | Proportion of successive NN intervals differing by more than 50 ms |
| rMMSD | Root mean square of successive NN intervals |
| MedianNN | Median of NN intervals |
Spectral analysis is a standard in HRV analysis. Features are extracted from the power spectral density (PSD) of the NN interval time series. The PSD within various frequency bands is estimated using Welch's method. The most common frequency bands under investigation are the very low frequency (VLF) band, the low frequency (LF) band, and the high frequency (HF) band. The spectral power distributions across these bands has been shown to provide insight into modulations in autnomic nervous system (ANS) activity.
freqDomainFeats = frequencyDomain(RR_interval_series)| Label | Description |
|---|---|
| VLF Power | Log of normalized spectral power between 0.003 Hz and 0.04 Hz |
| LF Power | Log of normalized spectral power between 0.04 Hz and 0.15 Hz |
| HF Power | Log of normalized spectral power between 0.15 Hz and 0.4 Hz |
| LF/HF Ratio | Ratio between LF and HF spectral power |
The frequencyDomain function also has the capability to compute the frequency domain features listed in the table above using the spectral boundary adaptation method described in X. Long, et al (2014). The idea behind using adapted spectral bands is to account for the time varying behavior of these bands (e.g. whilst sleeping) and to potentially reduce the within – and between – subject frequency domain feature variabilities. This method works by adapting the bands according to the peak frequencies found within the traditional spectral bands and setting their lower and upper bounds to values based off specified bandwidths.
Poincare plots are an important visualization technique for quantifying the non-linear characteristics of the RR interval time series. The geometrical descriptors that can be extracted from the poincare plot with this package have been shown to provide insights into short and long term HRV trends. This includes parameters derived from the ellipse fitting method described in Brennan, et al (2001) and the heart rate asymmetry (HRA) method that quantifies accelerations/decelerations in heart rate introduced in Guzik, Przemyslaw, et al (2006).
Multi-fractal detrended fluctuation analysis (MF-DFA) introduced in Kantelhardt, et al. MF-DFA is based on the standard detrended fluctuation analysis (DFA) introduced by Peng, et al (1995).
from DFA import scalingExponent
scaleExp = scalingExponent(timeSeries, lowerScaleLimit = 10, upperScaleLimit = 40, scaleDensity = 30, m = 1, q = 2, plot = 1)
The algorithm involves dividing the integrated RR interval time series into non-overlapping segments of equal length, s.
A polynomial is then fitted to each non-overlapping segment. Linear, quadratic, cubic, and/or other higher polynomials may be used in the fitting procedure. These DFA orders differ in their ability to eliminate various types of trends from the time series. Thus, an estimation of the type of time series trend can be captured via comparison of the different DFA orders, as seen below.
