Determination of maximum noise level in an ECG channel under SURE Wavelet filtering for HRV extraction

Heart Rate Variability (HRV) is the measure of variation between R-R interbeats, it has been demonstrated to be a good representation of physiological features, especially to the alterations in the Autonomic Nervous System (ANS). Considering the values that compose a HRV distribution are extracted from electrocardiography (ECG), many of the electrical disturbances that affect ECG-based diagnosis can also interfere with the results of the HRV analysis. This paper uses a 30-minute portion of a healthy patient (no arrhythmias detected or annotated) from the MIT-BIH ECG database to analyze the effectiveness of the SURE Wavelet denoising method for extracting the HRV from a progressively noisier ECG channel. Results show that the minimum SNR for reliable HRV extraction under these conditions is approximately 5dB and outlines the exponential behavior of HRV extraction for escalating noise levels in the ECG signal.


FIGURE 1.
A pair of heartbeat cycles as observed by the electrocardiograph, indicating the major 4 points of reference in electrocardiographic analysis: the T and P waves, as well as the QRS complex. The period between two R-peaks is also represented in the picture and is of great importance for HRV acquisition [6] .

INTRODUCTION
Heart Rate Variability (HRV) is currently the most common digital application of analysis of the ECG signal trace. It refers to a precise time measurement (usually in the order of milliseconds) of the intervals between each heartbeat, represented in Figure 1, labeled as R-R interval [1] .
The object of interest when analyzing variations in heart rate is autonomic nervous system (ANS), which is a division of the peripheral nervous system that acts in the regulation of bodily functions such as heart rate, respiratory rate and digestion [2] .
It is largely regarded as the primary control mechanism of the fight-or-flight reaction (or hyper-arousal) [3] .
In the past few decades, researchers have discovered correlation of the HRV analysis with many different health related issues, ranging from the autonomic function of insulin-dependent diabetic patients [4] to the aftereffects of myocardial infarction [5] .
The most usual method of extracting HRV data is the detection of R-wave peaks in an ECG trace plot.
Therefore, when dealing with acquired signals, the HRV can be compromised by the occurrence of artifacts originated by electrode motion, respiration motion (baseline wander) or even irradiated electromagnetic noise [7] . Considering that the relevant information to the HRV analysis is the reliable detection of the R peak, the design of filtering methods for this purpose is slightly different from the popular methods used for ECG filtering [8] .
The signal to noise ratio (SNR) is a dimensionless ratio comparing the signal and noise power in a given channel [9] . (1) This ratio is described in Equation 1 in its most usual form, where P s stands for the power of the signal being analyzed and P N stands for the power of the noise source. The result is displayed in decibels (dB) [10] .
The electrical engineering academic consensus is that 3dB in signal-to-noise ratio is the minimum value for the system to be considered acceptable in terms of the fidelity of signal output [11] , This threshold was obtained using Equation 1.
The Wavelet Transform is represented in its most general form in Equation 2 for a given function (x), having ψ as a continuous function for the mother Wavelet, accepting two parameters: a for scaling and b as the translational value. It can be used as a denoising method and is especially useful dealing with sharp transients in the signal.
Different methods using Wavelet have been developed over the years. Wavelet Shrinkage provides a great filter with low computational power [12] . Some other well-known Wavelets denoising methods are known as Bayes, BlockJS, Minimax, SURE and Universal Threshold [13] , all of which have been gaining on popularity since being incorporated into programming libraries on different platforms such as MATLAB and Python.
This paper proposes to analyze the process of extracting the HRV from a progressively noisier ECG channel in order to find out what is the minimum SNR with which SURE Wavelet denoising is still sufficient. The aim is to determine the limitation of this method of filtering specifically for HRV analysis.

MATERIALS AND METHODS
The ECG signal used for the HRV analysis was obtained from the MIT-BIH Arrhythmia ECG database [14] . This database was selected due to the long monitoring duration of each single patient (6-hour periods) available, enabling the collection of a sizeable sample of R-R intervals. From the records available, a 30-minute ECG record from a healthy patient (i.e. without observed arrhythmias) was chosen.
A short sample of the signal used is represented in Figure 2. It is important to consider that this signal has not yet went through the R-R identifying algorithm.  BioSigKit V1.0 [16] which was used for the extraction of the R-peaks through the Pan-Tompkins algorithm.
Signal was sampled at 360 Hz and filtered with the SURE Wavelet Transform as a denoising method with 5 levels of decomposition. This particular method of denoising was chosen because of its availability in most Wavelet toolboxes, its low demand in processing power and its good performance for biophysiological signal processing (especially for ECG and EEG applications) [17] .
The R-peak detection algorithm (specifically the Pan-Tompkins implementation [16] ) is implemented onto the signal represented in Figure 3 after it is filtered.
Points marked with red in Figure 3 were used to create the HRV vector. The process is a simple subtraction between the time values associated with two consecutive R-peak points, repeating the process all along the record. This is done once in order to obtain the real mean HRV (which consists of an arithmetic mean of all of the R-R values contained in the HRV vector) observed in the original 30-minute ECG record. The ECG signal is then summed with a matrix with random noise, this process is repeated in order to generate an ECG matrix with escalating noise level (each time noise power is incremented by 0.1).
Every newly generated noisy signal then undergoes the same SURE Wavelet method of denoising presented before, which in turns created the final matrix. points in a determinate process [18] .
From this final matrix, HRV is then extracted, also it is labeled with its SNR level using Equation 1 which will reflect the white noise added to it. This whole process is repeated 21 times, creating the matrixes in different noise-levels that will be needed for the analysis.
A possible way to represent HRV extracted from the progressively nosier channel is utilizing a Poincaré plot [18] . Which is a recursive plot since its axis represent a current measurement (called RR(n)) and a future measurement (RR(n+1)) repeated for many instances [18] . In a completely stable phenomenon, all points of the plot should be grouped together as tight as possible. When the subject is the human heart, that doesn't always happen for different reasons.
If the Poincaré in Figure 4 were data gathered from one subjects heart-rate variability, points observed in the extremities of the plot could either be the representation of some kind of intermittent cardiac function abnormality or the malfunctioning of the R-peak detection algorithm induced by noise in the signal.

RESULTS AND DISCUSSION
The chosen portion of the signal is filtered through the SURE Wavelet denoising tool and a small part of is represented in Figure 5. The Pan-Tompkins R-peak detection algorithm shown in Figure 3 is implemented on the filtered signal in order to extract the locations of every R-peak detected throughout the recording.
A Poincaré plot is used to represent the extracted HRV distribution in Figure 6. It useful here because it makes possible the visual analysis of increasing noise.
As expected from a filtered low-noise channel, the higher density of points in the center of the plot is visually recognizable.
Such artifacts are expected as HRV analysis usually focuses more on the total distribution of the plot and less on the outlier points [19] . The effects of the noisy ECG can be observed in the Poincaré plot of the HRV as represented in Figure 6.
The scattered points are more frequent and can be found further away from the centroids location.
In this case, calculated mean HRV is 274.298ms which is not a physiologically meaningful value, since it exceeds even the highest possible expected heart rate variation.  Experimental values as described in Table I are represented graphically in Figure 7 and a trend curve is fitted using the Smoothing Spline method.

AUTHOR CONTRIBUTIONS
R.N.C. conceived and designed the work, he participated in the data acquisition, initial analyses, data interpretation and first drafting of the writing. P.B.F.
Contributed substantially to the manuscript, reviewed and approved the final version, reviewed all aspects of the work, ensuring that questions related to the accuracy or integrity of the project were appropriately investigated and resolved. Both authors have equally contributed to the work in each stage of its development.