Analysis of heart rate variability based on quantitative approach

Heart rate variability (HRV) is the analysis of variations in the instantaneous heart rate in time series using beat-to-beat RR interval. HRV is considered as an indicator of the activity of autonomic regulation of circulatory function and as the one of the most significant methods of analyzing the activity of the autonomic nervous system (ANS). There are many methods that are used to detect heartbeats, such as ECG and blood pressure. Among them, ECG is considered as the best method to detect HRV signals since it provides a clear waveform, which makes it easier to exclude heartbeats not originating in the sinoatrial (SA) node.


Introduction
Heart rate variability (HRV) is the analysis of variations in the instantaneous heart rate in time series using beat-to-beat RR interval. HRV is considered as an indicator of the activity of autonomic regulation of circulatory function and as the one of the most significant methods of analyzing the activity of the autonomic nervous system (ANS). There are many methods that are used to detect heartbeats, such as ECG and blood pressure. Among them, ECG is considered as the best method to detect HRV signals since it provides a clear waveform, which makes it easier to exclude heartbeats not originating in the sinoatrial (SA) node.
The heart rate may be increased by slow acting sympathetic activity or decreased by fast acting parasympathetic (vagal) activity. The balance between the effects of the sympathetic and parasympathetic systems, the two opposite acting branches of the autonomic nervous system, is referred to as the sympathovagal balance and the reflected in the beat-to-beat changes of the cardiac cycle. Spectral analysis of the RR tachogram is typically used to estimate the effect of the sympathetic and parasympathetic modulation of the RR intervals.
The two main frequency bands of interest are referred to as the lowfrequency (LF) band (0.04 to 0.15 Hz) and the high-frequency (HF) band (0.15 to 0.4Hz).The ratio of the power contained in the LF and HF components has been used as a measure of the sympathovagal balance. The analysis of HRV signals is an important tool for studying the autonomic nervous system, as it allows the evaluation of the balance between the sympathetic and parasympathetic influences on heart rhythm.
Time-domain method is the simplest method to perform. Simple time-domain variables that can be calculated include the mean NN interval, the mean heart rate, the difference between the longest and shortest NN interval. In frequency domain methods, various spectral methods are used for the analysis of RR tachogram. Since RR intervals in time series are non-stationary and they are spaced unevenly, it needs to make re-sampling. The Lomb-Scargleperiodogram can overcome this weakness and it can estimate the PSD directly from irregularly sampled RR interval series without resampling. The most common classifications for HRV spectrum are: very low frequency (VLF: 0.003-0.04Hz), low frequency (LF: 0.04-0.15Hz) and high frequency (HF: 0.15-0.4Hz). The ratio of LF to HF appears to be a sensitive measure of the autonomic nervous system's (ANS) response to a sudden change in cardiovascular control and this ratio represents an evaluation of the ANS balance. [1][2][3][4][5][6][7][8][9][10] Overview of heart rate variability In this research, there are four major steps to implement the analysis of heart rate variability namely preprocessing, R-peak detection, and HRV analysis in time-domain and frequency-domain method.Raw ECG is taken from online MIT-BIH database, each data contains different ECG wave shape; each portion of the ECG waveform carries information that about the cardiac conditions. Preprocessing means the de-noising of raw ECG signal. During the recording of ECG signals, different types of noise from various sources can be superimposed to the original signal [12Nil]. In this reason, wavelet method is applied for the removal of baseline drift. Once the signal is free from noise the discrete wavelet transform is applied to extract the R-peak location and the amplitude in mV (millivolts) [16Sav,12Say]. The accuracy of the determined location of R-peak is essential for HRV analysis.
After detecting the R points, the RR interval analysis is performed. HRV metrics are calculated from an analysis of the RR tachogram, the time series of RR intervals. There are severalapproaches for HRV analysis, which can be subdivided into linear analysis (timedomain and frequency-domain) and non-linear analysis (wavelet). In this paper, the time-domain and frequency-domain method of linear analysis are used. The standard measurements of time-domain method are SDNN (ms), SDANN (ms), RMSSD (ms), NN50 (count), and PNN50 (%). These variables directly derived from the beat-to-beat intervals, such as the mean heart rate (HR) and the standard deviation (SD) for the entire record. The idea for the frequency domain methods is density (PSD) by using Lomb-Scargleperiodogram. After all, these results are verified with the normal values of standard measures of HRV. [11][12][13][14][15][16][17][18][19][20][21][22][23][24] Implementation of HRV analysis Autonomic nervous system (ANS) plays an important role in the regulation of the physiological processes of the human organism during normal and pathological conditions. Among the techniques used in its evaluation, the heart rate variability (HRV) has arising as a simple and non-invasive measure of the autonomic impulses, representing one of the most promising quantitative markers of the autonomic balance. The HRV describes the oscillations in the interval between consecutive heart beats (RR interval), as well as the oscillations between consecutive instantaneous heart rates. It is a measure that can be used to assess the ANS modulation under physiological conditions, such as wakefulness and sleep conditions, different body positions, physical training and also pathological conditions. Changes in the HRV patterns provide a sensible and advanced indicator of health involvements. Higher HRV is a signal of good adaptation and characterizes a health person with efficient autonomic mechanisms, while lower HRV is frequently an indicator of abnormal and insufficient adaptation of the ANS, provoking poor patient's physiological function. Because of its importance as a marker that reflects the autonomic nervous system activity on the sinus node and as a clinical instrument to assess and identify health involvements, this study reviews conceptual aspects of the HRV, measurement devices, filtering methods, indexes used in the HRV analyses [13 Hem], limitations in the use and clinical applications of the HRV.

Procedure of HRV analysis
Analysis of Heart rate variability (HRV) provides a non-invasive method to assess the neuronal influences on the cardiaoregulatory function. Since as defined before HRV is the fluctuation of RR intervals, these physiological fluctuations reflect the nonlinear feedback control systems created by the interaction between sympathetic and parasympathetic activities. The HRV analysis can be processed as shown in Figure 1. It includes four major steps namely signal de-noising, R-peak detection and quantification of HRV.

Procedure for ECG signal de-noising
A signal is often corrupted by noise during its acquisition or transmission. The de-noising process is to remove the noise while retaining and not distorting the quality of the processed signal. The traditional way of signal de-noising is filtering. From few years a lot of research about non-linear method of signal de-noising has been developed. These methods are mainly based on thresholding the Discrete Wavelet transform (DWT) coefficients. Simple de-noising algorithm that used DWT consists of three steps as shown in Figure 2. Firstly, DWT is adopted to decompose the noisy signal and get the wavelet coefficients. And these wavelet coefficients are de-noised with wavelet threshold. Finally, inverse transform is applied to the modified coefficients and get de-noised signal. The equations of ECG signal de-noising with DWT method is as follows.
The DWT of a signal x(t) is given by: where , m k ψ is the wavelet function.
where j is the decomposition level.
The reconstruction phase begins with the oversampling data followed by two reconstruction filters and to obtain the original signal as expressed by the Equation 4.

Estimation of threshold value and noise level
In the Wavelet transform method, the estimation of threshold value and noise level includes in the essential rule. There are various rules for thresholding such as soft thresholding and hard thresholding. In the research, global threshold calculating rule is usedand soft thresholding is applied for de-noising. The mean from the absolute values of the first detail coefficients ( ) 1 cD k then the standard deviation can be calculated by the Equations 5 and 6.

(5)
Where N is the number of the detail coefficients ( ) Finally, the IDWT can be computed using the new coefficients to reconstruct the de-noised signal.

Procedure for R-peak detection
For R-peak detection, the signal is decomposed into shifted and scaled versions of the original mother wavelet Db (6) and decomposed into level-8 because the wavelet Db6 closely matches with the shape of ECG QRS complex. An analysis of signal using ECG includes decomposition of the signal, thresholding of the wavelet coefficients and reconstruction of the signal using modified wavelet coefficients. And the following steps are managed for R-peak detection.
Step 1: ECG signal is read and the length is calculated.
Step 2: The signal is decomposed using db6 wavelet.
Step 3: 3 rd , 4 th and 5 th detail coefficients are selected, as most energy of the QRS complex is concentrated in these coefficients.
Step 4: The wave is reconstructed using detail coefficient 3, 4, 5 ( ) Step 5: A function d4*(d3+d5) /2^n is defined to reduce the oscillatory nature of the signal where d3, d4, d5 are the , , detail coefficients and n is the level of decomposition.
Step 6: Derivate up to level 5 is made using the transfer function Step 7: The differential equation: is applied to the signal, by using the transfer function and taking the amplitude response as Step 8: The signal is squared point by point using the equation (11) to emphasize R wave from the ECG signal.
Step 9: A moving window is integrated using the equation to obtain the waveform feature information.
Step 10: The threshold value is calculated corresponding to the product of max and mean of the signal to locate the end points of the moving window.
Step 11: The PQRST peaks are located based on the amplitude of the signal within each moving window.
Step12: The time intervals are calculated considering the positions of two consecutive same labeled peaks and stored.
Step 13: Diagnosis of various cardiac diseases is done by comparing ground truth conditions with the data.

HRV Metrics from the RR tachogram
The RR interval analysis consists in studying the Heart Rate Variability (HRV). The HRV is a measure of the variation of heart rate. It is usually calculated by analyzing the time series of beat to beat intervals from ECG (RR interval) or traces of blood pressure. It is obtained by measuring the time between RR intervals on the electrocardiogram. The values of RR intervalsare then plotted versus time, giving a curve called tachogram of HRV or RR tachogram as shown in Figure 3. This tachogram is a combination of sinusoidal waves of different frequencies.  HRV metrics are calculated from an analysis of the RR tachogram; the time series of RR intervals. It is important to note that this is an unusual time series in that both axes are time intervals, one being related to the other. Furthermore, since the variability in HR occurs on a beat-to-beat basis, the time series is inherently unevenly spaced along the horizontal axis. Figure 4 illustrates this concept; each star indicates the location of a beat in time (along the horizontal axis). The horizontal distance between each point (time stamp) is different for each adjacent pair, with the difference recorded on the vertical axis. The fact that the RR tachogram is inherently unevenly sampled leads to complications and errors in metrics that utilize interpolation.

HRV analysis in statistical indices of timedomain
Simple time-domain variables that can be calculated include the mean NN interval, the mean heart rate, the difference between the longest and shortest NN interval. For the analysis of RR intervals, we interest primarily in the extraction of beat frequency (F) in beats per minute (beats / min or rpm).
The first parameter is derived from the developed R peaks detection algorithm by the following formula: Beat frequency: Where L=total length of the recording (sample number), samp T =the sampling period(s), QRS N = number of detect QRS.
The root mean square of successive differences (RMSSD) is calculated for the purposes of HRV analysis by

HRV analysis by using lomb-scargle periodogram in frequency-domain
Since RR intervals in time series are non-stationary and they are spaced unevenly, it needs to make re-sampling. There are different interpolation methods for re-sampling. However, if the time series contains missing samples in heart rate time series, PSD estimates can be severely artificial and in such cases resampling is complicated by the need to infer probable values as replacements. The Lomb-Scargleperiodogram can overcome this weakness and it can estimate the PSD directly from irregularly sampled RR interval series without resampling.
This method is based on the definition of (DWT) as shown in Equation 16.
From this method, the power spectral density can be calculated several frequency bands of interest have been defined in humans.

Area calculation within frequency bands
The trapezoidal rule is a numerical method that approximates the value of a definite integral. We consider the definite integral as using the n + 1 point, The value of f(x) can be calculated at these points.
We approximate the integral by using n trapezoids formed by using straight line segments between the points ( )    x t

Tests and results
ECG de-noising is the initial procedure for heart rate variability analysis. Since the raw ECG loaded from the MIT-BIH database is normally corrupted with different types of noises. Therefore, it is needed to purify before R-peak detection stage. Figure 5 shows original ECG signal with base-line wander drift.

Wander baseline drift removal
Wander baseline drift and motion artifact belong to low frequency in which the wander baseline drift frequency is lower than 1 Hz. Since the ECG is a non-stationary signal, normal filters cannot be effective to remove the noise; so, several techniques are used to do so for such types of signals. Baseline wandering is a major source of noise which can be efficiently removed using wavelet based approach. In the original ECG signal, it needs to make de-trending and de-noising. Figure 6 shows the result of base line de-trending.

procedure with DWT
Generally, there are three steps for ECG signal de-noising by using Wavelet transform method. They are decomposition, thresholding detail coefficients and reconstruction. The decomposition of detail levels and the reconstruction of approximation and detail levels are shown in Figures 7-9 shows the base-line drift de-noised signal. The comparison of original and de-noised ECG signal is shown in Figure  10.

R-peak detection
In wavelet decomposition, the signal is decomposed into shifted and scaled versions of the original mother wavelet. An analysis of signal using ECG includes decomposition of the signal, thresholding of the wavelet coefficients and reconstruction of the signal using modified wavelet coefficients. The wavelet Db6 closely matches with the shape of ECG QRS complex. The ECG signals are decomposed into eight levels using Db6. Addition of third, fourth and fifth detail levels is taken as reconstructed signal as most energy of the QRS complex is concentrated in these coefficients. The peaks in reconstructed signal are observed in the vicinity of ECG R-peaks with maximum tolerance of six samples. All samples in reconstructed signal are scanned for peak. Figure 11 shows the R-waves that find in the restrictions of Wavelet Transform method and Figure 12 shows the detection of R-peak.

Overall performance of R peak detection
The performance of the R peak detection is essential for HRV analysis in detection of arrhythmia. Accuracies of R peak detection can be estimated according to the output results of 48 records. The detected R peaks are compared with MIT/BIH annotation records which represents the actual time intervals from patient. In matching with each annotation files, it is not needed that matching annotations have exactly equal times. While comparing with annotation, the performance of discrete wavelet transform (DWT) algorithm can be checked with the accuracy. From this result, we can see that DWT method for R peak detection is enough for HRV analysis.

Results for HRV analysis using time-domain method
Simple time-domain variables that can be calculated include the mean NN interval, the mean heart rate, the difference between the longest and shortest NN interval. The root mean square of successive differences (RMSSD) is calculated. Table 2 shows the results of statistical analysis in time-domain method.

RR intervals analysis
The RR interval analysis consists in studying the heart rate variability. It is usually calculated by analyzing the time series of beat to beat intervals from ECG (RR interval) or traces of blood pressure. It is obtained by measuring the time between RR intervals on the electrocardiogram. The values of RR intervals are then plotted versus time, giving a curve called tachogram of HRV or RR tachogram as shown in Figure 16. This tachogram is a combination of sinusoidal waves of different frequencies. From this RR-tachogram, this is an unusual time series in that both axes are time intervals, one being related to the other. Since the variability in HR occurs on a beat to beat basic, the time series is inherently unevenly spaced along the horizontal axis.

Conversion of RR iterval in time-series into PSD
Since RR intervals in time series are non-stationary and they are spaced unevenly, it needs to make re-sampling. There are different interpolation methods for resampling. However, if the time series contains missing samples in heart rate time series, PSD estimates can be severely artificial and in such cases resampling is complicated by the need to infer probable values as replacements. The Lomb-Scargleperiodogram can overcome this weakness and it can estimate the PSD directly from irregularly sampled RR interval series without resampling Figures 17-19 show the result of power spectrum for HRV analysis.

Results for HRV analysis using frequency-domain method
Since RR intervals in time series are non-stationary and they are spaced unevenly, it needs to make re-sampling. There are different interpolation methods for re-sampling. However, if the time series contains missing samples in heart rate time series, PSD estimates can be severely artificial and in such cases resampling is complicated by the need to infer probable values as replacements. The Lomb-Scargleperiodogram can overcome this weakness and it can estimate the PSD directly from irregularly sampled RR interval series without resampling. The results of spectral analysis in frequency-domain method are shown in Table 3.

Discussions on results
Heart rate variability (HRV) is a measure of variations of heart rate between two successive heart beats. It is usually calculated by analyzing the time series of beat-to-beat intervals from ECG signal. The detection of cardiac arrhythmias is a crucial point in the cardiac diseases diagnosis. An arrhythmia is characterized by the irregularity of the heart rate. A heart rate is regular if it is of the order of 60 beats per minute; otherwise; it's called bradycardia or tachycardia. The most commonly used modality for the arrhythmia diagnosis is the ECG. The detection of this cardiac irregularity is based on the R peaks detection and analysis of their regularity (RR intervals).
During the recording of ECG signals, different types of noise from various sources can be superimposed to the original signal. To purify the signal, many methods have been proposed: adaptive filtering, digital filtering and wavelet filtering. In this research, DWT algorithm is used for base-line wander drift de-noising. Using the knowledge that wavelet filtering is efficient and accurate in the compute of the R peaks positions without changing of the shape or position of the original signal, this technique is used to filter the ECG signal. For the choice of the mother function, a comparative study has achieved using various types of mother functions: coiflet, symlet and Daubechies. The best results are obtained by the use of the Db6 mother wavelet and decomposition level-8 is appropriate for de-noising. After denoising, the purified ECG signal is used as input for R-peak detection. For the R peak detection, DWT method is also used. Firstly, denoised ECG signal is read and length is calculated. ECG signals are decomposed into eight levels using Db6. After all, , and detail coefficients are selected and the wave is reconstructed using these coefficients. The peaks in reconstructed signal are observed in the vicinity of ECG R-peak with maximum tolerance of six samples. All samples in reconstructed signals are scanned for peaks.
The variation of heart rate is analyzed by using the linear method such as time-domain (statistical analysis) and frequency-domain (PSD) method. Simple time-domain variables that can be calculated include the mean NN interval, the mean heart rate, the difference between the longest and shortest NN interval. Since RR intervals in time series are non-stationary and they are spaced unevenly, it needs to make re-sampling. There are different interpolation methods for re-sampling. However, if the time series contains missing samples in heart rate time series, PSD estimates can be severely artificial and in such cases resampling is complicated by the need to infer probable values as replacements. The Lomb-Scargleperiodogram can overcome this weakness and it can estimate the PSD directly from irregularly sampled RR interval series without resampling. After all, the results of both methods are verified by normal value of standard measures of HRV. HRV plays an important role in an indication of cardiovascular health since it is considered as an indicator of the activity of autonomic regulation of circulatory function and as the one of the most significant method of analyzing the activity of the autonomic nervous system.

Conclusion
The quantification of HRV has been shown to afford a suggestion of cardiovascular health and can predict the hopeful different arrhythmia conditions such as Myocardial Infarction (MI), heart failure, angina, and sudden death. By applying the Lomb-Scargleperiodogram for HRV analysis, the variation of heart rate can be described more exactly. Moreover, the benefit of frequency-domain method (spectral analysis) is that it can make identification of different arrhythmia. The traditional time-domain method also yields the same results as abnormal. The performance of HRV analysis strongly depends upon the raw ECG signal de-noising and R-peak detection. By using DWT for these processes, the more accurate HRV analysis can be performed. And the HRV analysis also depends upon intrapatient basis (depending upon activity) and on an inter-patient basis (depending on their cardiovascular fitness).