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
where is duration of exercise and 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.
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 successive beat intervals, i.e. . The mean RR interval () and the mean heart rate () are then defined as
where denotes the value of ‘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
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
is a measure of short-term (beat-by-beat) variability. For stationary RR series and SDSD equals the root mean square of successive differences (RMSSD) given by
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
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
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.
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).
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
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
Please note that DC and AC sums above are divided by 2 (instead of 4, as originally proposed in Bauer et al. 2006) in order to reflect the absolute change in RR interval (in msec) similar to the modified indices.
|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. 1990, Keytel 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).
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.
|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. 2001, Carrasco et al. 2001, approximate and sample entropy Richman & Moorman 2000, Fusheng et al. 2001, detrended fluctuation analysis Peng et al. 1995, Penzel et al. 2003, correlation dimension Guzzetti et al 1996, Henry et al. 2001, and recurrence plots Webber et al. 1994, Trulla et al. 1996, Zbilut 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.
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 as a function of . 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.
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 , respectively.
The ellipse is oriented according to the line-of-identity (). 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
where SDSD is defined as , 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
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 (). In addition, the lag can be bigger than 1, e.g., the plot ().
Approximate entropy (ApEn) measures the complexity or irregularity of the signal Richman & Moorman 2000, Fusheng 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 vectors is formed
where is called the embedding dimension and 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.
Next, for each the relative number of vectors for which is calculated. This index is denoted with and can be written in the form
Due to the normalization, the value of is always smaller or equal to 1. Note that the value is, however, at least since is also included in the count. Then, take the natural logarithm of each and average over to yield
Finally, the approximate entropy is obtained as
Thus, the value of the estimate ApEn depends on three parameters, the length of the vectors , the tolerance value , and the data length . In Kubios HRV software the default value of is set to be . The length of the data also affects ApEn. When is increased the ApEn approaches its asymptotic value. The tolerance 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 is , which is also the default value in Kubios HRV software.
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 for which also the vector itself is included. This ensures that 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 is eliminated by calculating as
Now the value of will be between 0 and 1. Next, the values of are averaged to yield
and the sample entropy is obtained as
The default values set for the embedding dimension and for the tolerance parameter 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 , having repeated itself within a tolerance for points, will also repeat itself for 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
- 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 , where represents the scale factor and is selected to range between . The length of each coarse-grained time series is , where is the number of RR intervals in the data. For scale , the coarse-grained time series is simply the original beat-to-beat RR interval time series.
- SampEn is calculated for each coarse-grained time series. SampEn as a function of the scale factor produces the MSE. MSE for scale factor 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
where is the average RR interval. Next, the integrated series is divided into segments of equal length . Within each segment, a least squares line is fitted into the data. Let denote these regression lines. Next the integrated series 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
This computation is repeated over different segment lengths to yield the index as a function of segment length . Typically 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 (the slope of the regression line relating to . Different values of indicate the following
- : Brown noise (integral of white noise)
- : Different kinds of noise
- : noise
- :Large values are likely to be followed by large value and vice versa
- : White noise
- : 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 obtained from the (,) graph within range (default values). Correspondingly, the slope obtained by default from the range characterizes long-term fluctuations, see Fig. 6.
Figure 6: Detrended fluctuation analysis. A double log plot of the index as a function of segment length . and are the short term and long term fluctuation slopes, respectively.
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 vectors
and calculate the number of vectors for which , that is
where the distance function is now defined as
Next, an average of the term is taken
which is the so-called correlation integral. The correlation dimension is defined as the limit value
In practice this limit value is approximated by the slope of the regression curve 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 when is increased. In the software, a default value of was selected for the embedding.
Figure 7: Approximation of the correlation dimension from the (, ) 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
where is the embedding dimension and the embedding lag. The vectors then represent the interval time series as a trajectory in dimensional space. A recurrence plot is a symmetrical matrix of zeros and ones. The element in the ‘th row and ‘th column of the RP matrix, i.e. RP(,), is 1 if the point on the trajectory is close to point . That is
where is the Euclidean distance (see above) and 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.
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 (default value) and (fixed), respectively. The threshold distance was selected to be (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 is equal to and the recurrence rate is simply given as
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 is used for excluding the diagonal lines formed by tangential motion of the trajectory. The maximum line length is denoted and its inverse, the divergence,
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
where is the number of length lines. The determinism of the time series is measured by the variable
Finally, the Shannon information entropy of the line length distribution is defined as
where is the number of length lines divided by the total number of lines, that is
|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|
|MSE||–||Multiscale entropy for scale factor values τ=1,2,…,20 Costa et al. 2005|
- 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.
- 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.
- 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.
- 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.
- M. Costa, A.L. Goldberger, and C.-K. Peng. Multiscale entropy analysis of biological signals. Physical Rev E, 71:021906, 2005.
- 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.
- 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.
- P. Grassberger and I. Procaccia. Characterization of strange attractors. Phys Rev Lett, 50:346–349, 1983.
- 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.
- 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.
- 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.
- 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.
- S.L. Marple. Digital Spectral Analysis. Prentice-Hall International, 1987.
- 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.
- R.H. Morton, J.R. Fitz-Clarke and E.W. Banister. Modeling human performance in running. J Appl Physiol, 69(3):1171-7, 1990.
- 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.
- 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.
- 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.
- J.A. Richman and J.R. Moorman. Physiological time-series analysis using approximate entropy and sample entropy. Am J Physiol, 278:H2039–H2049, 2000.
- 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.
- 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.
- 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.
- 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.
- 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.
- 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.