HRV analysis methods

In this section, the analysis methods available in Kubios HRV software are introduced. The presented methods are mainly based on the guidelines given in Task Force 1996. The presentation of the methods is divided into four categories, i.e. time-domain, frequency-domain, nonlinear, and time-varying methods.

Time-domain HRV analysis methods

Energy expenditure and training load

Heart rate based energy expenditure models, provide a reliable estimates of daily energy expenditure. In Kubios HRV software the energy expenditure is divided into: 1) basal metabolic rate (BMR) computed using the Mifflin-St Jeor equationMifflin et al. 1990 and 2) activity related energy expenditure (EE) computed using the Keytel’s formulaKeytel et al. 2005. The energy expenditure computations are based on heart rate, body weight, height and age. In Kubios HRV, we compute the EE using beat-to-beat HR values, and thus, the instantaneous EE (kcal/min) can be realiably derived. The instantaneous values of EE can be used to assess how energy expenditure is distributed throughout the day.

Training impulse (TRIMP) is a rather commonly used index for training volume. The TRIMP can be computed according to exponential Banister’s model Morton et al. 1990, which is defined for male and female subjects separately as

    \begin{equation*} \rm TRIMP = T \times \Delta HR \times 0.64 e^{1.92\times\Delta HR},\quad {\rm (Male)} \end{equation*}

    \begin{equation*} \rm TRIMP = T \times \Delta HR \times 0.86 e^{1.67\times\Delta HR},\quad {\rm (Female)} \end{equation*}

where T is duration of exercise and \rm\Delta HR=\frac{HR_{ex}-HR_{rest}}{HR_{max}-HR_{rest}} is a heart rate reserve ratio. TRIMP accumulation rate increases exponentially as a function of exercise intensity, modelling lactate accumulation during exercise.

In Kubios HRV, we compute the TRIMP using beat-to-beat HR values, and thus, instantaneous value of TRIMP (TRIMP/min) can be reliably derived. TRIMP/min can be used as an index of training intensity. The cumulative sum of TRIMP/min, on the other hand, gives an index of training load accumulation. The link between TRIMP/min values and training intensity is shown in Fig. 1. The figure also shows for a 60-min training session the link between TRIMP value and training load.

Training impulse (TRIMP) for heart rate (HR) data

Figure 1: Training intensity levels corresponding to different Instantaneous TRIMP/min values (top panel) and training load levels corresponding to different cumulative TRIMP values over a 60-min training session (bottom panel).

Mean RR, SDNN, RMSSD and pNN50

The time-domain methods are derived from the beat-to-beat RR interval values in time domain. Let the RR interval time series include N successive beat intervals, i.e. {\rm RR}=({\rm RR}_1, {\rm RR}_2, \ldots, {\rm RR}_N). The mean RR interval (\overline{\rm RR}) and the mean heart rate (\overline{\rm HR}) are then defined as

    \begin{equation*} \overline{\rm RR} = \frac{1}{N}\sum_{n=1}^N {\rm RR}_n,\quad \overline{\rm HR} &=& \frac{60}{\overline{\rm RR}} \end{equation*}

where {\rm RR}_n denotes the value of n‘th RR interval. Several heart rate variability parameters that measure the variability within the RR time intervals in time-domain exist. The standard deviation of RR intervals (SDNN) is defined as

    \begin{equation*} {\rm SDNN} = \sqrt{\frac{1}{N-1}\sum_{n=1}^N ({\rm RR}_n-\overline{\rm RR})^2} \end{equation*}

SDNN reflects the overall (both short-term and long-term) variation within the RR interval time series, whereas the standard deviation of successive RR interval differences (SDSD) given by

    \begin{equation*} {\rm SDSD} = \sqrt{E\{\Delta{\rm RR}_n^2\} - E\{\Delta{\rm RR}_n\}^2} \end{equation*}

is a measure of short-term (beat-by-beat) variability. For stationary RR series E\{\Delta{\rm RR}_n\} = E\{{\rm RR}_{n+1}\}-E\{{\rm RR}_n\} = 0 and SDSD equals the root mean square of successive differences (RMSSD) given by

    \begin{equation*} {\rm RMSSD} = \sqrt{\frac{1}{N-1}\sum_{n=1}^{N-1} ({\rm RR}_{n+1} - {\rm RR}_n)^2}. \end{equation*}

Another measure calculated from successive RR interval differences is the NN50 which is the number of successive intervals differing more than 50 ms or the corresponding relative amount

    \begin{equation*} {\rm pNN50} = \frac{\rm NN50}{N-1} \times 100\% . \end{equation*}

Geometric measures

In addition to the above statistical measures, there are some geometric measures that are calculated from the RR interval histogram. The HRV triangular index is obtained as the integral of the histogram (i.e. total number of RR intervals) divided by the height of the histogram which depends on the selected bin width. In order to obtain comparable results, a bin width of 1/128 seconds is recommended Task Force 1996. Another geometric measure is the TINN which is the baseline width of the RR histogram evaluated through triangular interpolation, see Fig. 2A.

The Baevsky’s stress index (SI) is computed according to the formula Baevsky 2009

    \begin{equation*} {\rm SI} = \frac{{\rm AMo}\times 100\%}{2{\rm Mo\times MxDMn}} \end{equation*}

where AMo is the so-called mode amplitude presented in percent, Mo is the mode (the most frequent RR interval) and MxDMn is the variation scope reflecting degree of RR interval variability (see Fig. 2B). The mode Mo is simply taken as the median of the RR intervals. The AMo is obtained as the height of the normalised RR interval histogram (bin width 50 msec) and MxDMn as the difference between longest and shortest RR interval values. In order to make SI less sensitive to slow changes in mean heart rate (which would increase the MxDMn and lower AMo), the very low frequency trend is removed from the RR interval time series by using the smoothness priors method Tarvainen et al 2002. In addition, the square root of SI is taken to transform the tailed distribution of SI values towards normal distribution.

Geometric measures of HRV

Figure 2: Computation of geometric measures of HRV: A) triangular index (HRVi) and triangular interpolation of RR intervals (TINN), and B) Baevsky’s stress index.

HR deceleration and acceleration capacity

Deceleration capacity (DC) of heart rate is a measure of cardiac parasympathetic modulation as it captures the lengthening of RR interval within 2-4 successive beats as explained below. Acceleration capacity (AC) of heart rate captures the opposite, i.e. shortening of RR interval within few successive beats. These parameters were originally proposed in Bauer et al. 2006 and refined in Nasario-Junior et al. 2014.

Computation of DC (AC) is performed as follows. First, one RR interval value from every deceleration (acceleration) phase of RR interval time series is identified as an anchor interval. The anchor intervals are those that show steepest change compared to previous RR interval value, i.e. strongest deceleration (acceleration) within every deceleration-acceleration phase (see Fig. 3).

Deceleration capacity (DC) of heart rate

Figure 3: Illustration of heart rate deceleration capacity (DC) and acceleration capacity (AC) computations. The anchor RR intervals (●) and the surrounding RR intervals (○).

Next, all deceleration (acceleration) phases found from the RR interval data are aligned at the anchor points and an ensemble average is computed, i.e. phase rectified signal averaging (PRSA) is applied. Fig. 2 shows the ensemble average (black lines) and the standard deviation intervals (light-blue and light-grey areas) for both deceleration and acceleration phases. DC and AC are calculated from the ensemble averaged deceleration and acceleration phases as follows

    \begin{eqnarray*} {\rm DC} &=& \left[ {\rm RR}(0) + {\rm RR}(1) - {\rm RR}(-1) - {\rm RR}(-2)\right] / 2 \\ {\rm AC} &=& \left[ {\rm RR}(0) + {\rm RR}(1) - {\rm RR}(-1) - {\rm RR}(-2)\right] / 2 \end{eqnarray*}

That is, DC and AC are computed as a four point difference from the deceleration and acceleration PRSA signals, respectively. To better capture the instantaneous deceleration and acceleration, a modified DC and AC parameters are also computed according toNasario-Junior et al. 2014 as follows

    \begin{eqnarray*} {\rm DC}_{mod} &=& \left[ {\rm RR}(0) - {\rm RR}(-1) \right] \\ {\rm AC}_{mod} &=& \left[ {\rm RR}(0) - {\rm RR}(-1) \right]. \end{eqnarray*}

Table 1: Summary of time-domain HRV parameters calculated by Kubios HRV software.
EE [kcal/min] Energy expenditure (EE) divided into Basal Metabolic Rate (BMR) and activity related energy expenditure. BMR is estimated with Mifflin-St Jeor formula and energy expenditure (EE) using Keytel’s modelMifflin et al. 1990Keytel et al. 2005
TRIMP [1/min] Training impulse (TRIMP) is estimated using the exponential Banister’s modelMorton et al. 1990
[ms] The mean of RR intervals
STD RR (SDNN) [ms] Standard deviation of RR intervals
[beats/min] The mean heart rate
Min & Max HR [beats/min] Minimum and maximum HR computed using N beat moving average (default value: N=5)
RMSSD [ms] Square root of the mean squared differences between successive RR intervals
NNxx [beats] Number of successive RR interval pairs that differ more than xx ms (default value: xx=50)
pNNxx [%] NNxx divided by the total number of RR intervals
HRV triangular index The integral of the RR interval histogram divided by the height of the histogram
TINN [ms] Baseline width of the RR interval histogramTask Force 1996
Stress index Square root of Baevsky’s stress indexBaevsky 2009
DC, AC [ms] Deceleration capacity (DC) and acceleration capacity (AC) of heart rate computed as a four-point differenceBauer et al. 2006, Nasario-Junior et al. 2014
DCmod, ACmod [ms] Modified DC and AC computed as a two-point differenceBauer et al. 2006, Nasario-Junior et al. 2014
SDANN [ms] Standard deviation of the averages of RR intervals in 5-min segmentsTask Force 1996
SDNNI [ms] Mean of the standard deviations of RR intervals in 5-min segmentsTask Force 1996

Frequency-domain HRV analysis methods

In the frequency-domain methods, a power spectrum density (PSD) estimate is calculated for the RR interval series. The regular spectrum estimators implicitly assume equidistant sampling and, thus, the RR interval series is converted to equidistantly sampled series by interpolation methods prior to spectrum estimation. In Kubios HRV software, a cubic spline interpolation method is used. In HRV analysis, the spectrum is generally estimated using either a Fast Fourier transformation (FFT) based methods or parametric autoregressive (AR) modeling based methods. For details on these methods see, e.g., Marple 1987. The advantage of FFT based methods is the simplicity of implementation, while the AR spectrum yields improved resolution especially for short samples. Another property of AR spectrum that has made it popular in HRV analysis is that it can be factorized into separate spectral components.

In Kubios HRV software, the HRV spectrum is calculated with FFT based Welch’s periodogram method and with the AR method (see Fig. 4). Spectrum factorization in AR method is optional. In the Welch’s periodogram method the HRV sample is divided into overlapping segments. The spectrum is then obtained by averaging the spectra of these segments, which decreases the variance of the FFT spectrum.

Kubios HRV software includes also the Lomb-Scargle periodogram van Dongen et al. 1999, which differs from the Welch’s periodogram in the sense that it does not assume equidistant sampling and is thus computed directly from the non-interpolated RR interval time series. The variance of the Lomb-Scargle periodogram is decreased by smoothing the periodogram using MA-filering (the window width of the MA-filter can be adjusted in software preferences).

Heart rate variability (HRV) spectrum

Figure 4: HRV spectrum estimates using FFT based Welch’s periodogram method (left) and autoregressive (AR) modelling based spectrum estimation method.

The generalized frequency bands in case of short-term HRV recordings are the very low frequency (VLF, 0–0.04 Hz), low frequency (LF, 0.04–0.15 Hz), and high frequency (HF, 0.15–0.4 Hz). The frequency-domain measures extracted from a spectrum estimate for each frequency band include absolute and relative powers of VLF, LF and HF bands; LF and HF band powers in normalized units; the LF/HF power ratio; and peak frequencies for each band. In the case of FFT spectrum, absolute power values for each frequency band are obtained by simply integrating the spectrum over the band limits. In the case of AR spectrum, on the other hand, if factorization is enabled distinct spectral components emerge for each frequency band with a proper selection of the model order and the absolute power values are obtained directly as the powers of these components. If factorization is disabled the AR spectrum powers are calculated as for the FFT spectrum. The band powers in relative and normalized units are obtained from the absolute values as described in Table 2.

Table 2: Summary of frequency-domain HRV parameters calculated by Kubios HRV software.
Peak frequency [Hz] VLF, LF, and HF band peak frequencies
Absolute power [ms2] Absolute powers of VLF, LF, and HF bands
Absolute power [log] Natural logarithm transformed values of absolute powers of VLF, LF, and HF bands
Relative power [%] Relative powers of VLF, LF, and HF bands:
VLF [%] = VLF [ms2] / total power [ms2] x 100%
LF [%] = LF [ms2] / total power [ms2] x 100%
HF [%] = HF [ms2] / total power [ms2] x 100%
Normalized power [n.u.] Powers of LF and HF bands in normalised units:
LF [n.u.] = LF [ms2] / (total power [ms2] – VLF [ms2]) x 100%
HF [n.u.] = HF [ms2] / (total power [ms2] – VLF [ms2]) x 100%
LF/HF Ratio between LF and HF band powers
EDR [Hz] ECG derived respiration (available only if ECG data is used for HRV analysis)

Nonlinear HRV analysis methods

Considering the complex control systems of the heart it is reasonable to assume that nonlinear mechanisms are involved in heart rate regulation. The nonlinear properties of HRV have been analysed using measures such as Poincaré plot Brennan et al. 2001Carrasco et al. 2001, approximate and sample entropy Richman & Moorman 2000Fusheng et al. 2001, detrended fluctuation analysis Peng et al. 1995Penzel et al. 2003, correlation dimension Guzzetti et al 1996Henry et al. 2001, and recurrence plots Webber et al. 1994Trulla et al. 1996Zbilut et al. 2002. During the last years, the number of studies utilizing such methods have increased substantially. The downside of these methods is still, however, the difficulty of physiological interpretation of the results.

Poincaré plot

One commonly used nonlinear method that is simple to interpret is the so-called Poincaré plot. It is a graphical representation of the correlation between successive RR intervals, i.e. plot of {\rm RR}_{n+1} as a function of {\rm RR}_n. The shape of the plot is essential and a common approach to parameterize the shape is to fit an ellipse to the plot as shown in the figure below.

Poincare plot of heart rate variability (HRV) data

Figure 5: Poincaré plot analysis with the ellipse fitting procedure. SD1 and SD2 are the standard deviations perpendicular to and along the the line-of-identity {\rm RR}_n = {\rm RR}_{n+1}, respectively.

The ellipse is oriented according to the line-of-identity ({\rm RR}_n = {\rm RR}_{n+1}). The standard deviation of the points perpendicular to the line-of-identity denoted by SD1 describes short-term variability which is mainly caused by RSA. It can be shown that SD1 is related to the time-domain measure SDSD according toBrennan et al. 2001

    \begin{equation*} \rm SD1^2 = \frac{1}{2} SDSD^2 \end{equation*}

where SDSD is defined as {\rm SDSD}=\sqrt{E[\Delta {\rm RR}_n^2] - \overline{\Delta {\rm RR}_n}^2}, which is equal to RMSSD for stationary time series. The standard deviation along the line-of-identity denoted by SD2, on the other hand, describes long-term variability and has been shown to be related to time-domain measures SDNN and SDSD by Brennan et al. 2001

    \begin{equation*} \rm SD2^2 = 2\, SDNN^2 - \frac{1}{2} SDSD^2. \end{equation*}

The standard Poincaré plot can be considered to be of the first order. The second order plot would be a three dimensional plot of values ({\rm RR}_n,\RR_{n+1},{\rm RR}_{n+2}). In addition, the lag can be bigger than 1, e.g., the plot ({\rm RR}_n,{\rm RR}_{n+2}).

Approximate entropy

Approximate entropy (ApEn) measures the complexity or irregularity of the signal Richman & Moorman 2000Fusheng et al. 2001. Large values of ApEn indicate high irregularity and smaller values of ApEn more regular signal. The ApEn is computed as follows.

First, a set of length m vectors u_j is formed

    \begin{equation*} u_j = ({\rm RR}_j,{\rm RR}_{j+1},\ldots,{\rm RR}_{j+m-1}),\quad j= 1,2,\ldots N-m+1 \end{equation*}

where m is called the embedding dimension and N is the number of measured RR intervals. The distance between these vectors is defined as the maximum absolute difference between the corresponding elements, i.e.

    \begin{equation*} d(u_j,u_k) = \max\left\{\vert {\rm RR}_{j+l} - {\rm RR}_{k+l}\vert \,\bigl\vert\, l = 0,\ldots,m-1\right\}. \end{equation*}

Next, for each u_j the relative number of vectors u_k for which d(u_j,u_k)\le r is calculated. This index is denoted with C_j^m(r) and can be written in the form

    \begin{equation*} C_j^m(r) = \frac{{\rm nbr\, of\,}\left\{ u_k\,\bigl\vert\, d(u_j,u_k)\le r \right\}}{N-m+1}\quad\forall\, k. \end{equation*}

Due to the normalization, the value of C_j^m(r) is always smaller or equal to 1. Note that the value is, however, at least 1/(N-m+1) since u_j is also included in the count. Then, take the natural logarithm of each C_j^m(r) and average over j to yield

    \begin{equation*} \Phi^m(r) = \frac{1}{N-m+1}\sum_{j=1}^{N-m+1} \ln{C_j^m(r)}. \end{equation*}

Finally, the approximate entropy is obtained as

    \begin{equation*} {\rm ApEn}(m,r,N) = \Phi^m(r)-\Phi^{m+1}(r). \end{equation*}

Thus, the value of the estimate ApEn depends on three parameters, the length m of the vectors u_j, the tolerance value r, and the data length N. In Kubios HRV software the default value of m is set to be m=2. The length N of the data also affects ApEn. When N is increased the ApEn approaches its asymptotic value. The tolerance r has a strong effect on ApEn and it should be selected as a fraction of the standard deviation of the RR interval data (SDNN). This selection enables the comparison of RR data from different subjects. A common selection for r is r = 0.2 {\rm SDNN}, which is also the default value in Kubios HRV software.

Sample entropy

Sample entropy (SampEn) is similar to ApEn, but there are two important differences in its calculation Richman & Moorman 2000, Lake et al. 2002. For ApEn, in the calculation of the number of vectors u_k for which d(u_j,u_k)\le r also the vector u_j itself is included. This ensures that C_j^m(r) is always larger than 0 and the logarithm can be applied, but at the same time it causes bias to ApEn. In sample entropy the self-comparison of u_j is eliminated by calculating C_j^m(r) as

    \begin{equation*} C_j^m(r) = \frac{{\rm nbr\, of\,}\left\{ u_k\,\bigl\vert\, d(u_j,u_k)\le r \right\}}{N-m}\quad\forall\, k\neq j. \end{equation*}

Now the value of C_j^m(r) will be between 0 and 1. Next, the values of C_j^m(r) are averaged to yield

    \begin{equation*} C^m(r) = \frac{1}{N-m+1}\sum_{j=1}^{N-m+1} C_j^m(r) \end{equation*}

and the sample entropy is obtained as

    \begin{equation*} {\rm SampEn}(m,r,N) = \ln{(C^m(r)/C^{m+1}(r))}. \end{equation*}

The default values set for the embedding dimension m and for the tolerance parameter r in Kubios HRV software are the same as those for the ApEn computation. Both ApEn and SampEn are estimates for the negative natural logarithm of the conditional probability that a data of length N, having repeated itself within a tolerance r for m points, will also repeat itself for m+1 points. SampEn was designed to reduce the bias of ApEn and has a closer agreement with the theory for data with known probabilistic content Lake et al. 2002.

Multiscale entropy (MSE)

Multiscale entropy (MSE) is an extension of SampEn in the sense that it incorporates two procedures Costa et al. 2005

  1. A coarse-graining process is applied to the RR interval time series. Multiple coarse-grained time series are constructed for the time series by averaging the data points within non-overlapping windows of increasing length \tau, where \tau represents the scale factor and is selected to range between \tau=1,2,\ldots,20. The length of each coarse-grained time series is N/\tau, where N is the number of RR intervals in the data. For scale \tau=1, the coarse-grained time series is simply the original beat-to-beat RR interval time series.
  2. SampEn is calculated for each coarse-grained time series. SampEn as a function of the scale factor produces the MSE. MSE for scale factor \tau=1 returns standard SampEn (computed from the original data points).
Detrended fluctuation analysis

Detrended fluctuation analysis (DFA) measures the correlation within the signal. The correlation is extracted for different time scales as follows Peng et al. 1995. First, the RR interval time series is integrated

    \begin{equation*} y(k) = \sum_{j=1}^k ({\rm RR}_j - \overline{\rm RR}),\quad k = 1,\ldots,N \end{equation*}

where \overline{\rm RR} is the average RR interval. Next, the integrated series is divided into segments of equal length n. Within each segment, a least squares line is fitted into the data. Let y_n(k) denote these regression lines. Next the integrated series y(k) is detrended by subtracting the local trend within each segment and the root-mean-square fluctuation of this integrated and detrended time series is calculated by

    \begin{equation*} F(n) = \sqrt{\frac{1}{N}\sum_{k=1}^N (y(k)-y_n(k))^2}. \end{equation*}

This computation is repeated over different segment lengths to yield the index F(n) as a function of segment length n. Typically F(n) increases with segment length. A linear relationship on a double log graph indicates presence of fractal scaling and the fluctuations can be characterized by scaling exponent \alpha (the slope of the regression line relating \log{F(n)} to \log{n}. Different values of \alpha indicate the following

  • \alpha = 1.5: Brown noise (integral of white noise)
  • 1 < \alpha < 1.5: Different kinds of noise
  • \alpha = 1: 1/f noise
  • 0.5 < \alpha < 1:Large values are likely to be followed by large value and vice versa
  • \alpha = 0.5: White noise
  • 0 < \alpha < 0.5: Large value is likely to be followed by small value and vice versa

Typically, in DFA the correlations are divided into short-term and long-term fluctuations. In Kubios HRV software, the short-term fluctuations are characterized by the slope \alpha_1 obtained from the (\log{n},\log{F(n)}) graph within range 4\le n \le 12 (default values). Correspondingly, the slope \alpha_2 obtained by default from the range 13\le n \le 64 characterizes long-term fluctuations, see Fig. 6.

Detrended fluctuation analysis (DFA) of heart rate variability (HRV) data

Figure 6: Detrended fluctuation analysis. A double log plot of the index F(n) as a function of segment length n. \alpha_1 and \alpha_2 are the short term and long term fluctuation slopes, respectively.

Correlation dimension

Another method for measuring the complexity or strangeness of the time series is the correlation dimension which was proposed in Grassberger et al. 1983. The correlation dimension is expected to give information on the minimum number of dynamic variables needed to model the underlying system and it can be obtained as follows.

Similarly as in the calculation of approximate and sample entropies, form length m vectors u_j

    \begin{equation*} u_j = ({\rm RR}_j,{\rm RR}_{j+1},\ldots,{\rm RR}_{j+m-1}),\quad j= 1,2,\ldots, N-m+1 \end{equation*}

and calculate the number of vectors u_k for which d(u_j,u_k)\le r, that is

    \begin{equation*} C_j^m(r) = \frac{{\rm nbr\, of\,}\left\{ u_k\,\bigl\vert\, d(u_j,u_k)\le r \right\}}{N-m+1}\quad\forall\, k \end{equation*}

where the distance function d(u_j,u_k) is now defined as

    \begin{equation*} d(u_j,u_k) = \sqrt{\sum_{l=1}^m \left( u_j(l)-u_k(l)\right)^2}. \end{equation*}

Next, an average of the term C_j^m(r) is taken

    \begin{equation*} C^m(r) = \frac{1}{N-m+1}\sum_{j=1}^{N-m+1} C_j^m(r) \end{equation*}

which is the so-called correlation integral. The correlation dimension D_2 is defined as the limit value

    \begin{equation*} D_2(m) = \lim_{r\rightarrow 0}\lim_{N\rightarrow \infty} \frac{\log{C^m(r)}}{\log{r}}. \end{equation*}

In practice this limit value is approximated by the slope of the regression curve (\log{r},\log{C^m(r)}) Henry et al. 2001. The slope is calculated from the linear part of the log-log plot (see Fig. 7). The slope of the regression curves tend to saturate on the finite value of D_2 when m is increased. In the software, a default value of m=10 was selected for the embedding.

Correlation dimension (D2) for heart rate variability (HRV) data

Figure 7: Approximation of the correlation dimension D_2 from the (\log r, \log C^m(r)) plot.

Recurrence plot analysis

Yet another approach, included in the software, for analyzing the complexity of the time series is the so-called recurrence plot (RP) analysis. In this approach, vectors

    \begin{equation*} u_j = ({\rm RR}_j,{\rm RR}_{j+\tau},\ldots,{\rm RR}_{j+(m-1)\tau}),\quad j = 1,2,\ldots,N-(m-1)\tau \end{equation*}

where m is the embedding dimension and \tau the embedding lag. The vectors u_j then represent the {\rm RR} interval time series as a trajectory in m dimensional space. A recurrence plot is a symmetrical [N-(m-1)\tau]\times [N-(m-1)\tau] matrix of zeros and ones. The element in the j‘th row and k‘th column of the RP matrix, i.e. RP(j,k), is 1 if the point u_j on the trajectory is close to point u_k. That is

    \begin{equation*} {\rm RP}(j,k) = \left\{\begin{array}{ll} 1, & d(u_j-u_k)\le r \\ 0, & {\rm otherwise} \end{array}\right. \end{equation*}

where d(u_j,u_k) is the Euclidean distance (see above) and r is a fixed threshold. The structure of the RP matrix usually shows short line segments of ones parallel to the main diagonal. The lengths of these diagonal lines describe the duration of which the two points are close to each other. An example RP for heart rate variability time series is presented in figure below.

Recurrence plot analysis (RPA) of heart rate variability (HRV) data

Figure 8: Recurrence plot matrix for HRV time series (black = 1 and white = 0).

In Kubios HRV software, the embedding dimension and lag were selected to be m=10 (default value) and \tau=1 (fixed), respectively. The threshold distance r was selected to be \sqrt{m}\, {\rm SD} (default value), where SD is the standard deviation of the RR time series. These selection are similar to those made inDabire et al 1998. Methods for quantifying recurrence plots were proposed in Webber et al. 1994. The methods included in Kubios HRV software are introduced below.

The first quantitative measure of RP is the recurrence rate (REC) which is simply the ratio of ones and zeros in the RP matrix. The number of elements in the RP matrix for \tau=1 is equal to N-m+1 and the recurrence rate is simply given as

    \begin{equation*} {\rm REC} = \frac{1}{(N-m+1)^2}\sum_{j,k = 1}^{N-m+1} {\rm RP}(j,k). \end{equation*}

The recurrence rate can also be calculated separately for each diagonal parallel to the line-of-identity (main diagonal). The trend of REC as a function of the time distance between these diagonals and the line-of-identity describes the fading of the recurrences for points further away.

The rest of the RP measures consider the lengths of the diagonal lines. A threshold l_{\rm min} = 2 is used for excluding the diagonal lines formed by tangential motion of the trajectory. The maximum line length is denoted l_{\rm max} and its inverse, the divergence,

    \begin{equation*} {\rm DIV} = \frac{1}{l_{\rm max}} \end{equation*}

has been shown to correlate with the largest positive Lyapunov exponentTrulla et al. 1996. The average diagonal line length, on the other hand, is obtained as

    \begin{equation*} l_{\rm mean} = \frac{\sum_{l=l_{\min}}^{l_{\rm max}} l N_{l}}{\sum_{l=l_{\rm min}}^{l_{\rm max}} N_{l}} \end{equation*}

where N_{l} is the number of length l lines. The determinism of the time series is measured by the variable

    \begin{equation*} {\rm DET} = \frac{\sum_{l=l_{\rm min}}^{l_{\rm max}} lN_{l}}{\sum_{j,k=1}^{N-m+1} {\rm RP}(j,k)}. \end{equation*}

Finally, the Shannon information entropy of the line length distribution is defined as

    \begin{equation*} {\rm ShanEn} = -\sum_{l=l_{\rm min}}^{l_{\rm max}} n_{l} \ln n_{l} \end{equation*}

where n_{l} is the number of length l lines divided by the total number of lines, that is

    \begin{equation*} n_{l} = \frac{N_{l}}{\sum_{l'=l_{\rm min}}^{l_{\rm max}} N_{l'}}. \end{equation*}

Table 3: Summary of nonlinear HRV parameters calculated by Kubios HRV software.
SD1 [ms] In Poincaré plot, the standard deviation perpendicular to the line-of-identity Brennan et al. 2001, Carrasco et al. 2001
SD2 [ms] In Poincaré plot, the standard deviation along the line-of-identity
SD2/SD1 Ratio between SD2 and SD1
ApEn Approximate entropyRichman & Moorman 2000, Fusheng et al. 2001
SampEn Sample entropyRichman & Moorman 2000
DFA, α1 In detrended fluctuation analysis, short term fluctuation slopePeng et al. 1995, Penzel et al. 2003
DFA, α2 In detrended fluctuation analysis, long term fluctuation slope
D2 Correlation dimensionGuzzetti et al. 1996, Henry et al. 2001
RPA Recurrence plot analysisWebber et al. 1994, Dabire et al. 1998, Zbilut et al. 2002
  Lmean [beats] Mean line length
  Lmax [beats] Maximum line length
  REC [%] Recurrence rate
  DET [%] Determinism
  ShanEn Shannon entropy
MSE Multiscale entropy for scale factor values τ=1,2,…,20 Costa et al. 2005


  1. R. M. Baevsky. Methodical recommendations use kardivar system for determination of the stress level and estimation of the body adaptability standards of measurements and physiological interpretation. 2009.
  2. A. Bauer, J.W. Kantelhardt, P. Barthel, R. Schneider, T. Mäkikallio, K. Ulm, K. Hnatkova, A. Schöming, H. Huikuri, A. Bundle, M. Malik, and G. Schmidt. Deceleration capacity of heart rate as a predictor of mortality after myocardial infarction: cohort study. Lancet, 367:1647–81, 2006.
  3. M. Brennan, M. Palaniswami, and P. Kamen. Do existing measures of Poincaré plot geometry reflect nonlinear features of heart rate variability. IEEE Trans Biomed Eng, 48(11):1342–1347, 2001.
  4. S. Carrasco, M.J. Caitán, R. González, and O. Yánez. Correlation among Poincaré plot indexes and time and frequency domain measures of heart rate variability. J Med Eng Technol, 25(6):240–248, November/December 2001.
  5. M. Costa, A.L. Goldberger, and C.-K. Peng. Multiscale entropy analysis of biological signals. Physical Rev E, 71:021906, 2005.
  6. H. Dabire, D. Mestivier, J. Jarnet, M.E. Safar, and N. Phong Chau. Quantification of sympathetic and parasympathetic tones by nonlinear indexes in normotensive rats. amj, 44:H1290–H1297, 1998.
  7. Y. Fusheng, H. Bo, and T. Qingyu. Approximate entropy and its application in biosignal analysis. In M. Akay, editor, Nonlinear Biomedical Signal Processing: Dynamic Analysis and Modeling, volume II, chapter 3, pages 72–91. IEEE Press, New York, 2001.
  8. P. Grassberger and I. Procaccia. Characterization of strange attractors. Phys Rev Lett, 50:346–349, 1983.
  9. S. Guzzetti, M.G. Signorini, C. Cogliati, S. Mezzetti, A. Porta, S. Cerutti, and A. Malliani. Non-linear dynamics and chaotic indices in heart rate variability of normal subjects and heart-transplanted patients. Cardiovascular Research, 31:441–446, 1996.
  10. B. Henry, N. Lovell, and F. Camacho. Nonlinear dynamics time series analysis. In M. Akay, editor, Nonlinear Biomedical Signal Processing: Dynamic Analysis and Modeling, volume II, chapter 1, pages 1–39. IEEE Press, New York, 2001.
  11. L.R. Keytel, J.H. Goedecke, T.D. Noakes, H. Hiilloskorpi, R. Laukkanen, L. Van Der Merwe, and E.V. Lambert. Prediction of energy expenditure from heart rate monitoring during submaximal exercise. J Sports Sci, 23(3):289–297, 2005.
  12. D.E. Lake, J.S. Richman, M.P. Griffin, and J.R. Moorman. Sample entropy analysis of neonatal heart rate variability. ajp, 283:R789–R797, September 2002.
  13. S.L. Marple. Digital Spectral Analysis. Prentice-Hall International, 1987.
  14. M.D. Mifflin, S.T. St Jeor, L.A. Hill, B.J. Scott, S.A. Daugherty, and Y.O. Koh. A new predictive equation for resting energy expenditure in healthy individuals. Am J Clin Nutr, 51:241-7, 1990.
  15. R.H. Morton, J.R. Fitz-Clarke and E.W. Banister. Modeling human performance in running. J Appl Physiol, 69(3):1171-7, 1990.
  16. O. Nasario-Junior, P.R. Benchimol-Barbosa, and J. Nadal. Refining the deceleration capacity index in phase-rectified signal averaging to assess physical conditioning level. J Electrocardiol, 47(3):306–310, 2014.
  17. C.-K. Peng, S. Havlin, H.E. Stanley, and A.L. Goldberger. Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series. Chaos, 5:82–87, 1995.
  18. T. Penzel, J.W. Kantelhardt, L. Grote, J.-H. Peter, and A. Bunde. Comparison of detrended fluctuation analysis and spectral analysis for heart rate variability in sleep and sleep apnea. IEEE Trans Biomed Eng, 50(10):1143–1151, October 2003.
  19. J.A. Richman and J.R. Moorman. Physiological time-series analysis using approximate entropy and sample entropy. Am J Physiol, 278:H2039–H2049, 2000.
  20. M.P. Tarvainen, P.O. Ranta-aho, and P.A. Karjalainen. An advanced detrending method with application to HRV analysis. IEEE Trans Biomed Eng, 49(2):172–175, February 2002.
  21. Task force of the European society of cardiology and the North American society of pacing and electrophysiology. Heart rate variability – standards of measurement, physiological interpretation, and clinical use. Circulation, 93(5):1043–1065, March 1996.
  22. L.L. Trulla, A. Giuliani, J.P. Zbilut, and C.L. Webber Jr. Recurrence quantification analysis of the logistic equation with transients. Phys Lett A, 223(4):255–260, 1996.
  23. H.P.A. Van Dongen, E. Olofsen, J.H. VanHartevelt, and E.W. Kruyt. Searching for biological rhythms: peak detection in the periodogram of unequally spaced data. J Bioloogical Rhythms, 14(6):617–620, 1999.
  24. C.L. Webber Jr. and J.P. Zbilut. Dynamical assessment of physiological systems and states using recurrence plot strategies. J Appl Physiol, 76:965–973, 1994.
  25. J.P. Zbilut, N. Thomasson, and C.L. Webber. Recurrence quantification analysis as a tool for the nonlinear exploration of nonstationary cardiac signals. Med Eng Phys, 24:53–60, 2002.