Uncertainties in the Analysis of Heart Rate Variability: A Systematic Review

Heart rate variability (HRV) is an important metric with a variety of applications in clinical situations such as cardiovascular diseases, diabetes mellitus, and mental health. HRV data can be potentially obtained from electrocardiography and photoplethysmography signals, then computational techniques such as signal filtering and data segmentation are used to process the sampled data for calculating HRV measures. However, uncertainties arising from data acquisition, computational models, and physiological factors can lead to degraded signal quality and affect HRV analysis. Therefore, it is crucial to address these uncertainties and develop advanced models for HRV analysis. Although several reviews of HRV analysis exist, they primarily focus on clinical applications, trends in HRV methods, or specific aspects of uncertainties such as measurement noise. This paper provides a comprehensive review of uncertainties in HRV analysis, quantifies their impacts, and outlines potential solutions. To the best of our knowledge, this is the first study that presents a holistic review of uncertainties in HRV methods and quantifies their impacts on HRV measures from an engineer's perspective. This review is essential for developing robust and reliable models, and could serve as a valuable future reference in the field, particularly for dealing with uncertainties in HRV analysis.


I. INTRODUCTION
M ONITORING and analysing heart rate (HR) can provide valuable information about an individual's health status.The physiological basis and measurement of HR have been extensively studied in both healthy individuals and those with various diseases [1].Heart rate variability (HRV) refers to the variation in the time interval between consecutive heartbeats indicated by the oscillation in heart periods [2].It is considered a quantitative marker for assessing adequate cardiac regulation by the autonomic nervous system (ANS) in response to physical and psychological stimuli [3], [4].HRV is also affected by the activity of the sinoatrial node (SAN), the natural pacemaker of the heart [5].HRV is therefore widely used as a non-invasive parameter in healthcare applications, such as cardiovascular anomalies [6], critical care [7], and mental health disorders [8].
Generally, higher values of HRV parameters are associated with good health [9].For instance, a higher value of the root mean square of successive normal-to-normal intervals (RMSSD) has been linked to better self-rated health [10].Conversely, lower HRV values are often considered pathological and associated with reduced regulatory capacity [11].For example, the decline in frequency-domain HRV measures has been observed in older people with an increase in all-cause mortality [12].For these reasons, HRV has been used to monitor physical and psychological variables in various pathologies, such as critical care medicine [7], cancer patients [13], diabetes mellitus [14], bipolar disorder [15], irritable bowel syndrome [16], sleep apnea [17], myocardial infarction [6], silent stroke [18], depression [8], and psychotropic medication [19].It can be seen that numerous studies supported the use of HRV as an important tool for the diagnosis and monitoring a variety of clinical diseases.
The most common methods to obtain HRV data are through electrocardiography (ECG) and photoplethysmography (PPG) signals.However, ECG signals are often contaminated by noise and artifacts, which can subsequently affect the HRV analysis [20].On the other hand, PPG is an optical-based technique, which is similarly vulnerable to noise and artifacts, especially when data is collected from wrist-worn wearable devices [21].In addition, a variety of physiological factors can potentially affect HRV indices, such as age, weight, body mass index (BMI), body position, respiration rate, circadian rhythms, physiological states, and the effects of medications [13], [16], [22].Given these different types of uncertainties and their impact on HRV analysis, it is essential to quantify their effects to develop robust and reliable models for the analysis of HRV data.
The computation of HRV measures is a complex process that involves multiple computational procedures, including data denoising, outlier removal, cardiac cycle detection, missing data estimation, data segmentation, resampling, and spectral analysis [3].Each of these procedures can be implemented using different techniques, and their combination can lead to uncertainties in the analysis of HRV data.For instance, the HRV data can be analysed using short-term or long-term segments, but the two types of HRV measures are generally not interchangeable [23], [24].Given the impact of these computational techniques, there have been studies focusing on addressing specific types of uncertainties in the analysis of HRV data, such as noise and motion artifact removal [25], [26], [27], electromyogram interference [28], uncertainties in QRS complex detection [29], and the effect of limited sampling frequency [30].
Previous studies have developed a variety of techniques to simulate and analyse uncertainties in the analysis of HRV data, such as the spectro-temporal HRV analysis [31], Gaussian modelling [32], Monte Carlo simulation [33], Kalman filtering [34], non-linear noise reduction [35], and Bayesian deep learning [36].Given the different types of measurement uncertainties and computational techniques, it is important to provide a comprehensive review of these uncertainties and quantify their impacts on HRV indices.
This work aims to provide a comprehensive survey of uncertainties and computational techniques in HRV analysis, and quantify the impacts of these uncertainties on HRV measures from an engineer's perspective.This review informs researchers and practitioners about the various sources of uncertainty in HRV measurements and computational procedures, which would help facilitate the development of robust and reliable methods using HRV in healthcare applications.The review may also serve as a reference for the analysis of HRV data, providing potential solutions to mitigate the impact of uncertainties.

II. ACQUISITION AND MEASUREMENT OF HRV
HRV refers to the fluctuation in time intervals between adjacent heartbeats, which can be measured through different types of signal recordings [2].The standard approach to derive HRV data is through the use of ECG recordings, which utilise multiple electrodes placed on the skin to measure the electrical activities of the heart [37].Alternatively, the PPG technique uses an optical method to describe changes in pulse rate over time and is commonly used for measuring HRV data [38].Apart from the ECG and PPG signals, there are also other types of signals that can be used to derive HRV data, such as the ballistocardiograph (BCG) [39], heart sound signal [40], Doppler radar [41], and camera-based assessment, e.g., imaging-PPG and video-PPG [42].In the following sections, we describe two major types of signals (i.e., ECG and PPG) for HRV measurement.

A. HRV From ECG Signal
HRV can be determined by analysing HR data obtained from ECG signals, using specialized software for HRV calculation such as the PhysioNet's HRV Toolbox [43], Kubios HRV [44], and pyHRV [45].To derive HRV data from ECG signals, the R-peaks in the signal must be detected using appropriate QRS algorithms.For example, in Fig. 1(a), there are eight R-peaks in the ECG data; The beat-to-beat intervals can then be computed between the consecutive R-peaks; Next, HRV data can be obtained by processing the calculated RR intervals, using interpolation methods and data resampling to produce a consistently sampled RR tachogram for further analysis [46].

B. HRV From PPG Signal
In recent years, there has been a growing interest in developing alternative methods to measure HRV data using low-cost wearable devices [1].One of such methods is the use of PPG sensor, which consists of light sources and photodetectors to assess the amount of light absorbed or reflected by blood vessels in living tissue and detect volumetric changes in peripheral blood circulation [38].Then, the pulse rate variability (PRV) can be computed from the PPG signal, and previous studies have demonstrated the potential of using PRV as a surrogate for HRV data [1], [37].To derive the PRV data, as shown in Fig. 1(b), the pulse-to-pulse interval (PPI) or inter-beat interval (IBI) can be calculated from the PPG signal, and then it is straightforward to calculate the PRV data from the PPIs.The simplicity of this technique, cost-effectiveness, ease of signal acquisition, and remote monitoring capabilities are apparent advantages of PPG over the gold standard ECG for the estimation of HRV data.

C. Measurement of HRV
HRV analysis involves calculating meaningful variables from the RR or PP intervals in different time segments, such as long-term (24 hours or more), short-term (around 5 minutes), and ultra-short-term (less than 5 minutes) [37].However, longer recording periods are better suited for representing slower cardiac fluctuations and the response of the cardiovascular system to environmental stimuli and workloads [47].Therefore, shortterm and ultra-short-term HRV values are generally not interchangeable with long-term values [47].HRV variables can be evaluated using time-domain, frequency-domain, and non-linear analysis.We discuss some widely used HRV measures below, but for a more detailed description of these variables, we refer the reader to [3], [47].
Time-domain HRV indices quantify the variability in the time intervals between successive heartbeats.Widely used timedomain features are summarized in Table I, which include the mean value of normal-to-normal (NN) intervals (meanNN), the standard deviation of NN intervals over a 24-hour period (SDNN), the standard deviation of differences between adjacent NN intervals (SDSD), the standard deviation of the averages of NN intervals of all 5-minute segments (SDANN), the percentage of interval differences of successive NN intervals larger than 50 ms with respect to the total number of NN intervals (pNN50), and the HRV triangular index [47], [48], [49], [50].Because some of these time-domain measures are correlated, it is generally recommended to use specific indices for different HRV assessment purposes.For example, the SDNN and HRV triangular index can be used to estimate the overall HRV data, the SDANN can be used for the long-term components of HRV estimation, and the RMSSD can be used to estimate the short-term components of HRV data [3].
Frequency-domain analysis of HRV involves transforming the time-series data of RR intervals into the frequency domain, and the resulting power spectrum can be further categorized into several frequency bands as shown in Table II [47], [50].The power within each band reflects the activity of different physiological systems that influence heart rate, such as sympathetic and parasympathetic nervous systems [3].It is worth noting that frequency-domain HRV analysis can be performed across shortterm or long-term recordings.The recommended recording periods for frequency-domain HRV measures are as follows [47], [51]: the very-low-frequency (VLF) band (0.0033-0.04 Hz) requires a recording period of at least 5 minutes, but maybe over 24 hours; the low-frequency (LF) band (0.04-0.15 Hz) is typically recorded over a minimum period of 2 minutes; the high-frequency (HF) (0.15-0.40 Hz) is conventionally recorded over a minimum period of 1 minute.For the ultra-low-frequency (ULF) band, which is below 0.003 Hz, the index should be calculated for at least 24 hours or a longer period [51].
Non-linear HRV features provide insights into the complex non-linear dynamics of the cardiovascular system, which are influenced by various physiological variables, such as hemodynamic, electrophysiological, and humoral factors, as well as autonomic and central nervous regulations [3].Table III presents some widely used non-linear HRV measures [52], [53].Notably, some studies investigated the relationship between non-linear and spectral HRV analyses, and suggested that the coefficient α of the detrended fluctuation analysis (DFA) could be described by frequency-weighted spectral analysis, where the DFA method is a standard approach to assess long-range correlations embedded in non-stationary time series data [54].
We note that the same mathematical formula can be used to calculate both short-term and long-term heart rate variability (HRV) variables, but they have different physiological meanings and cannot be substituted for each other [55].For example, 24-hour data recordings capture slow fluctuations in circadian rhythms, core body temperature, metabolism, sleep cycle, and renin-angiotensin system contributions.As a result, HRV measures calculated from 24-hour recordings can guide biofeedback training for clinical and optimal performance [47].Longterm data recordings are typically used to assess autonomic nervous system (ANS) responses during normal daily activities in both healthy and unwell individuals, as well as in response to therapeutic interventions [56].In contrast, short-term HRV measures have practical advantages, such as easy application in out-of-hospital or laboratory settings and simplified data processing [57].As a result, short-term HRV measures are more commonly used as indicated in literature [47], [48], [52].Moreover, advances in wearable sensors and mobile monitoring systems have led to an increase in ambulatory acquired HRV data, resulting in a shift towards analysing ultra-short-term HRV features [57].This shift is largely due to real-life data constraints; even conventional 5-minute HRV recordings may be impractical to obtain in out-of-clinic environments [58].However, there is currently a lack of rigorous methods to assess the validity of short recordings and to identify reliable ultra-short HRV features in controlled settings.

III. UNCERTAINTIES IN HRV DATA ANALYSIS AND POTENTIAL SOLUTIONS
There are a variety of uncertainties that can affect the analysis of HRV data.For example, observational errors in the data acquisition, e.g., baseline wandering, motion artifacts, and electromyographic (EMG) interference; Meanwhile, there are also uncertainties arising from computational models in HRV analysis, e.g., QRS complex detection and data segmentation; In addition, some physiological factors such as age, gender, and BMI may also have impacts on the analysis of HRV data.
We note that there are different definitions of uncertainties.For example, the research in [59], [60] described uncertainty as the potential deficiency in any phase or activity of the process which can be characterized as not definite, not known, or not reliable.The uncertainties can be divided into model uncertainty (i.e., epistemic uncertainty), and data uncertainty (i.e., aleatoric uncertainty) [61].To be more specific, they can be summarised as uncertainties related to technology, incomplete knowledge, and clinical effectiveness [62].To this end, we sort out the different types of uncertainties in the analysis of HRV data, and divide them into three subsections, (i) uncertainties in data measurement, (ii) uncertainties in computational modelling, and (iii) influence on HRV indices from physiological factors.As shown in Fig. 2, we provide an overall diagram of reviewing the different types of uncertainties in the analysis of HRV data.In the following subsections, we review these uncertainties in details.

A. Uncertainties in Data Acquisition
As introduced earlier, HRV data can be derived from ECG or PPG signals.For both the two methods, the data acquisition process can be easily contaminated by noise and artifacts, e.g., motion artifacts.Although these contaminants can be reduced using hardware and software techniques during the experimental setup, it is unrealistic to remove all possible contaminants, and they would have effects on the subsequent HRV analysis.As summarised in Table IV, typical measurement uncertainties that occur during signal recording include baseline wandering, electromyogram (EMG) interference, electrode contact noise, lead placement, sampling frequency, and motion artifacts [25], [64], [65]; Meanwhile, there are measurement uncertainties in using different light sources and measurement locations for PPG signal acquisition [65], [66].These errors can lead to poor signal quality and affect the analysis of HRV data.
Uncertainties in ECG Measurement: 1) Baseline Wandering in ECG Measurement: The baseline wandering (BW), also known as drift, is a commonly seen noise during signal recording, which is usually caused by respiration or movement of the subject [67].For example, the BW in ECG measurement can be caused by the rotation of the heart's electrical axis during the respiratory cycle, which changes the thoracic impedance distribution and affects the ECG measurements [21].
Generally, the BW noise would affect the detection of peak values in the signal.For example, if the BW is large, it may lead to clipping of R peaks, resulting in immeasurable HRV [68].The BW noise is a low-frequency artifact at frequencies ranging from 0.15 to 0.3 Hz (i.e., <1 Hz) [69], and represented as a sinusoidal component added to the measured signal.Therefore, the BW noise can be processed using a high pass filter with a cutoff frequency greater than 1 Hz [70].There are also a variety of advanced filtering techniques that can be used to process the BW noise, such as adaptive filtering [71], and signal decomposition techniques [72].
2) Powerline Interference: The powerline interference in ECG measurements is mainly generated by the displacement current through a capacitive coupling between the human body and the power line [73].The interference is usually characterised as a high-frequency noise with the frequency of 50 ± 0.2 Hz (or 60 Hz), and an amplitude of up to 50% of peak-to-peak amplitudes in the measured signal [69].In a conventional ECG configuration with three or more electrodes, the powerline interference can be reduced using an additional reference electrode

TABLE IV MEASUREMENT UNCERTAINTIES IN HRV ANALYSIS
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
or the driven-right-leg (DRL) technique [74]; In the wearable ECG monitoring with two electrodes, the powerline interference may corrupt the quality of ECG recording severely [73], and the noise can be mitigated using some advanced hardware design, such as the ECG amplifier with a common-mode charge pump (CMCP) [75], and a biosignal read-out circuit with common-mode impedance compensation [76].Soft computational techniques have also been widely used to reduce powerline interference, such as digital filters (e.g., FIR and IIR filters) [77], adaptive filters [78], wavelet transform (WT) [79], and empirical mode decomposition (EMD) [72].
3) EMG Interference: The EMG interference is generated from muscle contractions, and it usually has an average amplitude of 10% full scale deflection (FSD) on the measurements, i.e., ECG data [25].The EMG signal has a broad range of frequency bands, and it sometimes overlaps the useful frequency range of ECG signal [80].As a result, a simple band-pass filter may be not adequate to remove the interference.Given that the EMG signal collected on the skin surface is localised in nature, the EMG interference in different ECG leads may be uncorrelated because of their different locations on the body.Therefore, the EMG interference can be suppressed by measuring multiple ECG leads combined with adaptive filtering techniques [81], [82].

4) Electrode Pop or Contact Noise:
The electrode contact noise is caused by the loss of contact between the electrode and skin, for example, the electrode being nearly or completely pulled off the skin [20].Electrode pops are common during body movements, where the signal manifests as sharp changes with periods of around 1 s.These errors could be eliminated by adding electrode gel or reapplication of the problem electrode when taking a recording.Afterwards, the artifact can generally be handled by switching to an alternative signal that does not come from the problem electrode [83].

5) Placement of ECG Leads:
The standard 12-lead ECG recordings are measured from different locations of the chest (precordial), wrists, and ankles on the body surface.Typically, the ECG leads I, II, and V2 are widely used for HRV analysis [85].There are studies investigating the effect of measurement locations on the QRS complex of ECG signals [95], [96], where the detection of QRS complex has substantial effects on the quality of HRV data.In particular, by analysing the clinical body surface potential map (BSPM) data with 120-lead ECGs, the study showed that the best locations for QRS complex are around the chest electrodes of the standard precordial V2, V3, and V4 ECG leads [95]; The results provide useful indications for the development of wireless ECG measurement systems.

6) Sampling Frequency of ECG Signal:
The sampling frequency at which the ECG signal is digitised is important in deriving HRV indices [97].For example, low sampling frequencies could distort the R-peak waveform, and the distortion will then propagate during QRS detection and therefore affect the HRV measures.In particular, the pNN50 is demonstrated as a sensitive HRV measure for low sampling frequencies [23], [30].Some research showed that the uncertainty of spectral HRV indices is proportional to the inverse of the sampling frequency, and the bias is proportional to the inverse of the square sampling frequency [98].The Task Force of The European Society of Cardiology and The North American Society of Pacing and Electrophysiology indicated that a low sampling frequency might produce jitters in estimating the R-wave fiducial points [3], which will then considerably alter the HRV spectral estimation, and thus the optimal range of sampling frequency is recommended as 250 to 500 Hz or perhaps higher [3].
7) Motion Artifacts in ECG Measurement: Motion artifact (MA) refers to the uncertainty caused by the relative movement of the sensor measurement with respect to the skin [99], [100].MA has long been a problem in both the ECG and PPG measurements, which will in turn affect the HRV or PRV data that are derived from these signals.However, the ECG signal measures physiologic biopotentials with a frequency range from 0.05 (or 0.5) Hz to 100 (or 150) Hz [25], [101], and the PPG uses an optical approach to measure the volumetric variations of blood circulation with a frequency band from 0.5 to 4 Hz [88], [89].Therefore, the MA has different impacts on the two types of signals, and we discuss them in different subsections.
In particular, MA during the ECG measurement is introduced by the electrode motion, which causes deformations of the skin around the electrode site, and in turn causes changes in the electrical characteristics of the skin around the electrode [20].MA poses a major challenge in the long-term cardiac monitoring using wearables [102], and it is more abrupt and distinct in nature as opposed to the slow baseline wandering caused due to respiration [70].
ECG signal can be measured with wet or dry electrodes [87], where the wet electrode such as Ag/AgCl electrode is widely used in medical applications with excellent signal acquisition quality, it uses a conductive gel that acts as the electrolyte of the electrode and reduces the contact impedance; In comparison, the dry contact electrode is widely used in wearable device, it connected directly to the skin and therefore is more sensitive to MA [87].Previous study showed that the light abrasion of skin with fine sandpaper is an effective way to reduce MA [103], [104]; other techniques such as extracting electrode-tissue contact resistance and reactance can be used to monitor MA [87], [105], and accelerometers can be used as reference signals to reduce the MA [106].In addition, there are many advanced signal processing techniques that can be used to reduce the effect of MA, which will be discussed in the latter section.
Uncertainties in PPG Measurement: 1) Baseline Wandering in PPG Measurement: The BW noise in the PPG measurement is also caused by the respiratory rate, muscle tremor, and physiological changes.However, different from the chances of impedance distribution during ECG measurement, the BW in PPG measurement is due to changes in tissue blood volume, which can be caused by the changes in intrathoracic pressure transmitted through the arterial tree, and vasoconstriction of arteries during inhalation transferring blood to the veins [21], [107].Similarly, the above discussed techniques in processing BW of ECG measurements can also be applied to process the noise in PPG data [71], [72].
2) Measurement Locations of PPG Signal: In terms of PPG signal, the measurement locations will also have an influence on signal quality, and consequently affect the estimation of HRV data [38], [66].Previous studies investigated the influence from different measuring sites of PPG signals, such as the right thumb, right forefinger, right middle finger, left thumb, left forefinger, left middle finger, right wrist (posterior), right wrist (anterior), left wrist (posterior), left wrist (anterior), forehead, nasal bridge, right ear, left ear, right ankle, left ankle, manubrium, and xiphoid process [66], [108].Among these different locations for PPG measurement, it was found that the fingers, ears, and nose reveal higher perfusion than other sites.However, the selection of suitable measuring sites is related to the choice of light sources.For instance, the pulse amplitude of the green-light PPG signal at the upper arm was larger than the peripheral sites; while the infrared amplitude shows no significant differences across measurement sites in arms, wrists, and fingers [66].
3) Effects of Light Sources on PPG Signal: PPG sensors emit light onto the skin, and measure the pulse rate by detecting the amount of light transmitted or reflected.The quality of PPG signal has a strong relationship with the light sources, because it estimates changes in blood volume based on the light absorption characteristics of hemoglobin.Therefore, the variations of different wavelengths of light sources have effects on the HRV data that are derived on the PPG signal [109].For example, some studies investigated the effect of different light sources of PPG signals, including blue, green, red; and near-infrared light in the reflectance mode, and red and near-infrared light in the transmittance mode [65], [110].The results suggested that blue and green lights followed by near-infrared light in reflectance mode are the recommended settings to measure HR using PPG techniques [110].In addition, some research used a multi-wavelength PPG module by including different types of lights, i.e., blue, green, yellow and infrared; and it was demonstrated that the model achieved more accurate estimation than traditional PPG techniques [111].

4) Sampling Frequency of PPG Signal:
The sampling rate also has effects on the HRV data when using PPG signal.Previous studies investigated the impact of different sampling frequencies (i.e., 1,000, 500, 250, 125, 100, 50, and 25 Hz) on acquiring PPG signals, where the data were acquired from the forehead and finger using red and infrared lights [112].The study showed that the relative errors in the time and frequency HRV indices increase with the decrease of sampling rates [112].There are also studies indicating that for time domain analysis 100 Hz sampling rate is required for accurate PRV analysis without interpolation in normal variability series, and for frequency domain analysis 10 Hz can be sufficient with cubic spline interpolation [113].

5) Motion Artifacts in PPG Measurement:
The quality of PPG signal is highly influenced by MA, this is because MA generally has a frequency range of 0.01-10 Hz [114], which overlaps the frequency range for PPG signals of 0.5-4 Hz [88], [89], making it difficult to remove MA without affecting the PPG signal; The uncertainty from MA will in turn affect the subsequent PRV analysis.Different from MA in ECG signal, the interference in PGG signal is generated due to the movement between the photoelectric sensor and the skin, which can be caused by voluntary or involuntary subject movement, such as people's unconscious movement, respiration, and extrusion between the skin and the photoelectric sensor [115].Generally, MA in PPG data can be denoised using a reference signal collected from accelerometers [116], optical sensors [117], or impedance sensors [118]; and it also can be compensated using some advanced sensor design, such as the dual-channel organic photodiode (OPD) sensor as presented in [119].
Numerous signal processing techniques have been developed to reduce MA in both the ECG and PPG signals, such as filter-based techniques, time-frequency analysis, blind source separation, and machine learning-based algorithms [120], [121].As shown in Table V, we summarise these different types of techniques for MA removal.
(a) MA Removal Using Adaptive Filtering: Adaptive filtering techniques generally use a reference signal that has a high correlation with motion, and then it is possible to weaken the MA and obtain a clean signal.The reference signal may be obtained by accelerometers [116], skin stretch sensors [129], optical sensors [117], and impedance sensors [118] as mentioned above.Nevertheless, if the reference signal is not available from an auxiliary sensor, synthetic techniques can also be used to produce a reference signal [130].With the reference signal at hand, the adaptive least mean square (LMS) algorithm [122], recursive least square (RLS) adaptive filter [90], and adaptive tracking (AT) [123] can estimate the clean signal.The adaptive filter is easy to implement, but the major drawback of this technique is the need of a supplementary signal, which increases the complexity of the system [120].
(b) MA Removal Using Time-Frequency Techniques: Timefrequency techniques have been extensively used to remove artifacts in biomedical signals, as they can present both the time and frequency information of a signal.The widely used wavelet transform (WT) technique decomposes a signal into a series of sub-signals with multiple lower resolutions, and then removes the artifacts by processing the wavelet coefficients [25].The empirical mode decomposition (EMD) algorithm is another useful non-stationary signal processing technique.It decomposes the signal into multiple intrinsic mode functions with well-defined instantaneous frequencies, and enables estimating the HR data from the distorted signals with the decomposed components [91].Other techniques such as the short time Fourier transform (STFT) [131], and the variational mode decomposition (VMD) [132] has also been used for MA removal.
(c) MA Removal Using Blind Source Separation: The blind source separation (BSS) refers to the recovery of source signals from a set of observed mixtures with superimposed noise [133].The BSS model is generally expressed algebraically in terms of matrix factorisation, and the original signal sources and MA can be separated once an unmixing matrix is determined by the BSS model [120].Based on some assumptions on the source signals and the mixture matrices, different approaches of BSS can be implemented, including the independent component analysis (ICA) [124], principal component analysis (PCA) [92], canonical correlation analysis (CCA) [93], and non-negative matrix factorisation (NMF) [125].
(d) MA Removal Using Other Advanced Techniques: Machine learning approaches have been demonstrated as powerful techniques for MA removal [134].The support vector machine Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE V TECHNIQUES TO REDUCE MOTION ARTIFACTS IN HRV DATA ANALYSIS TABLE VI UNCERTAINTIES IN HRV ANALYSIS FROM COMPUTATIONAL MODELS
was applied to identify high confident heartbeats from MA corrupted ECG signal in wearable applications [127].A Bayesian learning-based probabilistic approach was developed to estimate HRV from artifacts distorted PPG signal recorded by wearable devices [94].A real-time modelling method was developed for MA detection and PPG signals reconstruction using real-time modelling based on multilayer perceptron (MLP), radial basis function (RBF) artificial neural networks (ANN), and adaptiveneuro fuzzy inference system (ANFIS) [126].A deep learning model with the convolutional neural networks (CNN) and long-short term memory (LSTM) was developed to estimate HR information and perform biometric identification from a singlechannel PPG signal collected in an ambulant environment [128].We note that these different types of models are usually not implemented independently for MA removal, instead, they can be implemented combining with other techniques (e.g., signal decomposition and adaptive filtering) to enhance the abilities of MA removal [91].

B. Uncertainties in Computational Modelling
After data acquisition, the HRV data can be obtained using a variety of computational procedures, such as cardiac cycle detection, missing value imputation, data filtering, resampling, segmentation, and spectral estimation.However, these computational processes may produce potential uncertainties in the analysis of HRV data [135].For example, an unreliable detection of peak values from ECG or PPG signals would affect the HRV parameters, and the length of data segmentation would also have an impact on the time-or frequency-domain HRV indices [24], [136].As shown in Table VI, we summarise the uncertainties from computational modelling, and we discuss the details in the following subsections.

1) QRS Complex Detection of ECG Signal:
The accurate detection of QRS complex from the ECG signal is the basis for the efficient analysis of HRV data [151].Previous investigated the impact of the QRS detection process on timedomain, frequency-domain, and non-linear HRV parameters using the Monte Carlo simulation technique [135].The results indicated that the pNN50 and frequency-domain indices were the most sensitive indices; while the most robust indices were time-domain parameters as well as the short-term and long-term slopes of the DFA index [135].We note that many techniques have been developed for the detection of QRS complex, such as the amplitude threshold technique [137], the first derivative method [138], digital filtering [139], time-frequency representation [140], spectral analysis [152], geometry analysis [141], and learning-based techniques [153].For the different types of QRS detectors, it was suggested that the time-domain techniques are straightforward to compute but lack sufficient information, and time-frequency representation approaches could be considered a better option over geometry-based methods.In addition, previous study suggested that more reliable HRV results can be obtained by combining multiple QRS detection techniques [29].

2) Pulse-to-Pulse Interval Detection of PPG Signal:
Similar to the detection of QRS complex in ECG signal, the detection of PPIs is also important for the analysis of PPG-based HRV data [154], [155].Widely used detectors for the PPG signal include adaptive threshold approach [142], template matching technique [143], moving average filter [144], and waveform delineator [145].In particular, previous study compared three different techniques for pulse rate detection, and the results indicated that the traditional maximum peak detector had a high accuracy for the pulse rate estimation; the moving average filter can be used as an alternative algorithm; while the template matching technique had the worst performance [144].We note that many different fiducial points can be used for the detection of pulse-to-pulse intervals, such as the apex point, line-medium point, interpolate medium point, and the basal point [112], [146].It was demonstrated that when using the line-medium points as fiducial points for the finger PRV signal with a sampling rate of 25 Hz, the estimation of the power of LF index has a small estimation error [112].However, when using line-medium or medium interpolate points as fiducial points, there are no significant changes in HRV indices between data sampled from 50 Hz and 1,000 Hz for finger and forehead PPGs [112].Moreover, it was indicated that the middle-amplitude point, the apex point of the first derivative, and the tangent intersection point were the most suitable fiducial points for the PRV data analysis, which obtained low estimation errors between PRV and HRV indices [146].
3) Missing Data and Estimation Techniques: Missing data is a common type of uncertainty in the analysis of HRV data.Some research evaluated the effect of missing values on HRV data, and compared HRV measures derived from RR intervals with missing data and indices from a complete 5-min measurement [147]; and the results indicated that in terms of the root-mean-squared relative error (RMSRE) with a threshold of less than 10%, it was acceptable that the duration of missing data was approximate 140.6 s for the SDNN, 175.6 s for the RMSSD, 35.5 s for the pNN50, 285 s for the meanNN, and 0.9 s for almost all the frequency-domain HRV indices [147].
We note that there are a variety of techniques that can be used to impute missing values in the HRV data.For example, it would be straightforward to estimate missing data using interpolation methods, such as the nearest neighbour approach (NNR), linear, quadratic, spline cubic, and piecewise cubic Hermite (PCH) interpolation methods [33], [149].In terms of frequency-domain HRV features, previous research indicated that the spline and PCH interpolation methods have better performance than the NNR and linear interpolation [33]; In particular, the PCH interpolation has better estimation performance than the spline method on estimating the LF index; and the spline interpolation produces a low error than the PCH interpolation for a small duration of missing data on estimating the HF index.However, the errors caused by spline interpolation are larger than those caused by PCH method when the missing data duration is more than 60 s [33].Previous study investigated the impact of missing data on time-domain HRV features, such as the meanNN, SDNN, SDSD, RMSSD, and pNN50; The results indicated that the meanNN was the most robust parameter to missing data, and pNN50 was the most sensitive index to missing data [148].
In addition, some studies investigated the impact of missing values and interpretation techniques on estimating non-linear HRV measures, such as the Poincaré plot, detrended fluctuation, and entropy analysis [156].The results showed that among these non-linear HRV parameters, SD1 and SD2 obtained by Poincaré plot analysis were the most robust parameters to the missing data.However, comparing with time-and frequency-domain HRV parameters, the estimation errors of non-linear HRV parameters were relatively high.Therefore, apart from the SD1 and SD2 parameters derived from the Poincaré plot analysis, other parameters were not recommended for an accurate non-linear HRV analysis with missing data [156].
In terms of imputation strategies for missing data estimation, previous studies investigated two types of different approaches, i.e., interpolation of RR intervals on time (i.e., the timestamp when the heartbeats happen) and on the duration (i.e., the duration of the heartbeats) [149], and the results showed that the interpolation of missing data on time can produce more reliable HRV estimations than the interpolation on duration [149].Other than using interpolation techniques for missing data estimation, we note that machine learning methods such as deep recurrent neural networks, and mathematical tools have also been used to estimate missing values for the analysis of HRV data [157].

4) Frequency Rate for Data Resampling:
The identified intervals between heartbeats are usually non-uniform sampled, and the sequence of the time intervals is generally resampled for spectral analysis.The resampling rate can range from 1 to 10 Hz, however, the 4 Hz remains the preferred rate for spectral HRV applications [150], [158].Previous studies compared the resampling rates of 1 Hz and 4 Hz for spectral HRV analysis, and showed that there were small estimation errors between the two resampling frequencies for spectral indices [158]; In particular, compared the rate of 1 to 4 Hz, the normalized LF (LFnu), normalized HF (HFnu), and the ratio of LFnu to HFnu obtained the mean relative errors (MRE) of 3.7%, 15.3% and 16.4%, respectively [158].In addition, by investigating a variety of resampling frequencies (e.g., 2, 4, 6, 8, and 10 Hz) for HRV data analysis [150]; the study suggested that the selection of resampling rate should consider the mean and minimum heartbeat intervals.For example, the resampling rate of 8 Hz and above would be appropriate for subjects with HR > 117 bpm, e.g., newborns; the rate of 6 Hz is suitable for subjects with 90 ≤ HR ≤ 117 bpm; and the rate of 4 Hz is sufficient for HR < 90 bpm [150].

5) Duration of Data Segmentation:
HRV can be analysed with data segments from short durations to very long recordings, depending on the features of interests.For example, the HRV task force suggested that the short-term 5-min recordings and nominal 24-hour long-term recordings would be appropriate for time-domain HRV measures [3].However, for the frequencydomain analysis, the VLF assessed from short-term recordings is dubious, and should be avoided when interpreting the power spectral density (PSD) of short-term ECGs [3].Some research indicated that long-term HRV analysis is more stable than shortterm analysis, but it is more expensive and time-consuming to obtain long-term recordings.Previous study suggested that the short-term HRV analysis can be used to track the changes of cardiovascular autonomic function across minutes, while Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
the long-term HRV analysis should be preferred for depicting changes across hours to even more extended periods [23].
There was also research investigating HRV parameters derived from ultra-short-term (UST) data segments, i.e., <5 mins [24].By comparing HRV parameters derived from data segments with 10, 20, 30, 60, 90, 120, 180, and 240 s recordings, the research suggested that the segment of 60 s achieved an acceptable estimation of SDNN, RMSSD, and pNN50; the segment of 90 s can be used to estimate the LF power, SD1, and SD2; and the segment of 120 s can be used to estimate the HRV triangular index and DFA α 1 [24], [136].However, it was indicated that factors such as age, health conditions, and techniques for artifact processing may have a greater impact on UST HRV measurements than on longer recordings [47]; and therefore the routine use of UST HRV measurements in medicine, performance, and personal fitness assessment need more detailed investigations [159].
6) Methods for Spectral Estimation: Generally, there are two types of approaches that can be used to estimate the spectral of HRV data, i.e., non-parametric and parametric methods [33], [160].In the non-parametric method, the PSD can be calculated directly from the signal data, and there are two widely used non-parametric methods; (i) The first type of non-parametric methods uses the fast Fourier transform (FFT) to calculate the spectral power; (ii) The second type of non-parametric methods is based on Welch's periodogram [33], the method uses a window function (e.g., Hamming window) to calculate periodogram for each data segment, and then computes the averaged periodogram for the data [160].In the parametric method, the PSD is calculated from the frequency response of the transfer function of a linear system [33], and there are also two types of parametric methods; (i) The first parametric approach is the Yule-Walker method, it uses the lagged-product autocorrelation to estimate the parameters of AR model [161]; (ii) The second parametric approach is the Burg method, which fits the AR model by minimizing the forward and backward prediction errors [33].It is noted that other than using the resampled time series data for HRV analysis, the Lomb-Scargle (LS) periodogram is also used for the estimation of HRV metrics, which requires no resampling of unevenly sampled signals [33].Both the non-parametric and parametric methods are widely used to calculate PSD for spectral analysis of HRV data.However, when comparing the two types of PSD estimation methods, previous research indicated that the parametric methods can provide smoother spectrum than the non-parametric methods [160].In terms of the two types of parametric methods, the Burg method was demonstrated with good properties, e.g., guaranteed stability, and therefore the method is often advised for the estimation of AR model in practice [161].

C. Influence From Physiological Factors
In addition to the uncertainties in data measurement and computational methods, as discussed in the previous sections, there are also physiological factors that would affect the analysis of HRV data, such as age, gender, and BMI [162].For example, it was reported that age plays an important role in women's cardiac autonomic modulation, followed by the BMI, and these modulations will have an influence on the analysis of HRV data [163].We provide a brief summary of the impacts of physiological factors in Table VII, and discuss the details in the following subsections.
1) Age: It is well established that the values of HRV indices decrease with a progressing age [164].For example, by investigating the impact of confounding factors on HRV indices in 653 subjects without heart diseases, and it indicated that most of time-and frequency-domain HRV parameters decrease with ageing [162]; In particular, the decreasing rate was 6.6 ms/decade for the SDNN, and 4.4 ms/decade for the SDANN [162], [165].Previous study also indicated that almost all short-term HRV indices were significantly age-dependent, and the greatest influence was observed in the population with ages ranging from 25 to 54 years [48].
We note that many existing studies focused on investigating HRV indices for adults.However, HRV parameters can be quite different for the pediatric population, in particular infants/neonates.This is because children usually exhibit higher HR and higher variations of HR than adults, which will have effects on the subsequent analysis of HRV data [166].Previous study investigated normal subjects with ages from 1 mo to 24 years in different conditions, such as awake, active, and quiet sleep, and it observed the age dependence of HRV; In particular, an increase in LF, HF, and total power from 0 -6 years, followed by a decrease to 24 years [167].
2) Gender: Gender dependencies were observed in several short-term HRV indices, and significant differences were observed between males and females for the resting HR, RRintervals, and frequency-domain HRV parameters [164].It was reported that females had higher total power and HF power, but lower LF power and LF/HF ratio than males; However, the gender difference would gradually disappear after the age of 40 -50 years old [23], [48].A meta-analysis examined 50 HRV measures for healthy participants, and showed that compared to men, women had higher mean HR, lower SDNN, and SDNN index values, particularly for 24-hour studies [169]; For the frequency-domain HRV indices, the results showed that women had lower total, VLF, and LF power, but greater HF power [47], [169].However, some research indicated no significant gender differences in time-domain HRV parameters [164].For example, the SDNN had a higher value in females than in males, but it did not reach a statistically significant level [168]; Meanwhile, it showed that compared to age, the gender difference was indicated with less impact on HRV indices [168].Therefore, it was suggested that we should keep in mind the effect of gender differences particularly in the spectral analysis of HRV data [23].
3) BMI: The BMI is considered as another factor of intervention in autonomic system regulation, which will clinically translate into HRV parameters.For example, previous research found an increase in the mean NN, SDNN, SDANN indices, and LF and HF components after weight loss [170].In particular, it was observed that the SDNN, pNN50, and RMSSD had lower values in women with BMI values less than 19 kg/m 2 or greater than 30 kg/m 2 [163].However, another study indicated that there was no significant difference in either the time-or frequencydomain HRV indices between the normal, overweight, and obese groups [162].Therefore, the impact of BMI on HRV parameters needs more detailed investigation.

A. Characterisation of Measurement Noises
Given the distinctive features and impacts of the variety of measurement noises, it is important to characterise the nature of different types of these interferences.Conventionally, except for the sudden abrupt interference, most noises can be assumed to be Gaussian (white) with flat frequency spectra [20].However, many real-world noises are non-stationary, and they have nonflat spectral density functions, e.g., coloured noise [171], [172].
The term of coloured noise refers to any non-white noise signal with non-constant spectra, and the spectral density function of the coloured noise is a function of the frequency [173].For example, the pink noise is known as flicker noise, which presents in almost all electronics [174].The spectra of pink noise are inversely proportional to the frequency, and it indicates that the pink noise is a significant noise source at low frequencies [175], [176].The log power spectrum provides an index to characterise the coloured noise with its slope such that [20], [177], where, S(f ) represents the PSD function, and f is the spectral frequency.The parameter β is the slope of the log power spectrum, which determines the colour of the signal.
r If β = 0, the signal is regarded as white (Gaussian) noise, and it has a flat spectrum; r If β = 1, it is known as pink noise (flicker noise), and is related to the observation noise; r If β = 2, it is the electrode movement noise with a Brow- nian motion-like form, known as brown noise.The white noise is a sequence of uncorrelated random variables with zero mean and finite variance, and it has an equal intensity at different frequencies and a constant PSD value; In contrast, the pink noise is the signal with a frequency spectrum inversely proportional to the frequency of the signal [178].During the data acquisition, the chopper stabilisation technique is widely used to reduce pink noise.Other techniques, such as auto-zeroing and correlated double sampling, can also be used to reduce pink noise [175].Other than the pink noise, the Brown noise relates to the baseline wander and electrode motion artifacts [171], [176], and the removal of MA has been extensively discussed in Section III-A.
An illustration of the impact of coloured noise is demonstrated in Fig. 3.The original ECG signal is shown in Fig. 3(a).We use the simulator as described in [173] to generate the coloured noise, and the value of signal-to-noise ratio (SNR) is set as 3 dB.Fig. 3(b) indicates the contaminated ECG signal by the pink noise, and Fig. 3(c) shows the contaminated signal by the brown noise.The PSD analysis of the two types of noises is demonstrated in Fig. 3(d).It can be seen from Fig. 3(d) that the slopes of the pink and brown noises match well with the characteristics as described in (1).When comparing Fig. 3(b) and (c), given the SNR of the two types of noises, it seems the ECG signal is more contaminated by the pink noise than the brown noise.This finding agrees with previous study that the whiter the noise, the more significant distortion for a given SNR value [20].
In addition, previous research indicated that the ultra-low and very low frequencies of the HRV spectra over 24 hours had a characteristic of 1/f shape that comes from the activity of the sinoatrial node [179]; While the ANS contribution might be modelled as a stochastic process of white noise with two periodic components at the characteristics baroreceptor and Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
respiratory peaks in the low-and high-frequency bands [5].However, it is noted that the shape of the spectrum may be altered due to measurement uncertainties, for example, the presence of missing data dampens low frequencies and enhances high frequencies [149].

B. Uncertainty Quantification and Sensitivity Analysis
Uncertainty quantification (UQ) and sensitivity analysis (SA) are techniques that can quantify how these uncertainties will affect the model performance, and therefore the techniques can help to develop trustworthy models [180].The UQ and SA have been investigated in many scientific areas, such as cardiovascular modelling [181], material science [182], mechanical systems [183], and astrophysics [184].However, there are limited studies focusing on quantifying uncertainties for the analysis of HRV data.For example, a previous study used artificially generated data to quantify the impact of resampling and beat replacement in HRV spectral analysis [185].This section outlines the concept, procedures, and techniques for the implementation of UQ and SA in various areas.
1) The Concept of UQ and SA: The UQ is a process that determines uncertainties in estimated parameters given uncertainties in the model formulation and experimental measurements, as well as estimates how uncertainties in model inputs affect the model output [186].The UQ can be generally described by statistical moments, percentiles, and confidence intervals, which can be calculated by the frequentist approach, Gaussian process, and Bayesian method [186], [187].While the SA includes the global and local assessment [180], and the aim of SA is to quantify the contribution to the output of particular uncertain inputs and their interactions.The SA can be described by the first and second-order sensitivity indices, and the total sensitivity indices to consider the interaction effects [188].
2) Procedures of UQ and SA: Some research outlined the procedures for the UQ and SA in system modelling [189].For example, a two-step method was suggested for the analysis.First, the SA is conducted to identify the key parameters whose tolerances contribute the most to the parametric uncertainty of the selected design variables; Second, quantifying the effect of increasing the model inputs on the estimated total uncertainty, and the predictive capability of the simulation models [189].In addition, a more detailed framework was introduced for the UQ and SA, the workflow consisted of the following steps [188], including (i) Identification of the output of interest; (ii) Identification and assessment of the distribution of the uncertain inputs; (iii) Sampling of the input space to acquire samples; (iv) Evaluation of the deterministic model; (v) Calculation of UQ and SA measures; and (vi) Assessment of convergence of UQ and SA measures.
3) Techniques for UQ and SA: Numerous techniques have been developed to qualify uncertainty and analyse sensitivity in system modelling.The most straightforward technique for UQ and SA is the stochastic simulation, the method has been applied to capture uncertainties in the quantities that determine model behaviour, and provide probabilistic representations of model parameters [190].The technique has been used to quantify the effect of inherent uncertainties from the cardiac output on the sensitivity of a human compliant arterial network response, and investigate the sensitivity of the pulse pressure and waves reflection magnitude over the arterial tree for different model uncertainties [191].
As an example, the patient-specific simulations were performed to investigate the UQ and SA of left ventricular (LV) function during the full cardiac cycle [192], which were implemented with the following procedures, (i) the research first developed the cardiac simulation with three models, i.e., constitutive model, active stress model, and circulatory model; (ii) then, the polynomial chaos expansion method was used to accelerate the UQ analysis [193], where the uncertainties were considered in the constitutive model.Statistical measures computed from the polynomial expansion were used to quantify the impact of the input uncertainties; and (iii) the Sobol sensitivity indices were used to perform the SA analysis [193], i.e., investigating the contribution of input parameters in the variability of model outputs.The implementation of UQ and SA in this research quantified the influence of uncertainties in the input parameters and geometry of a cardiac computational model on the prediction of clinical interests [192].
The Multilevel Monte Carlo method has also been applied for UQ in stochastic multi-scale systems.The technique has been used to evaluate the impact of computational uncertainties from a limited sampling rate and the process of QRS detection in HRV analysis [135].Previous research indicated that the simulation improved the precision of identifying polynomial chaos expansions versus the traditional heuristic approach, and lowered the computational cost of UQ [194].In another example, the Gaussian process-based Markov chain Monte Carlo approach was used to quantify uncertainties of model parameters in cardiac electrophysiology; the model provided posterior distributions of model parameters, and it can be used to identify potential factors contributing to parameter uncertainties [195].
In addition to the above simulation techniques, the Bayesian approach is also an efficient technique in performing the UQ and SA.The uncertainties can be first evaluated using statistical methods, e.g., moment estimation.Then the Bayesian approach can be used to compute the posterior probability for the quantification [196].The method was used for model optimisation with respect to the predictive mean and variance, which takes into account both the data density and measurement uncertainties [184].Other approaches, such as the fuzzy-set theory was applied to estimate the probabilities of uncertainties of the model input [183]; Mathematical techniques, e.g., the Lipschitz constant estimator, were used to calculate the bounds of parameter uncertainties [187].

V. DISCUSSION
HRV has emerged as an important indicator in healthcare, and there are many existing reviews focusing on summarising clinical applications of HRV measures [7], [13], [15], [16].It is understood that HRV indices are derived from biomedical signals using a variety of computational techniques.However, uncertainties from the signal measurement and computational Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
approaches will have substantial effects on HRV analysis.We therefore summarise these different types of uncertainties in HRV, and provide the review as a reference to engineers when processing HRV data.To the best of our knowledge, this would be the first comprehensive review to cover different sources of uncertainties that may affect the analysis of HRV data, quantify their impact, and outline the potential solutions.
There are a variety of signals that can be used to calculate HRV data, and we limited the discussion of measurement uncertainties from ECG and PPG signals in this review, which are two major sources to derive HRV data.It is understood that ECG and PPG sensors use distinct principles for data measurement, and therefore they can be affected by different types of measurement uncertainties.For example, the electrical interference may produce undesirable effects on ECG signals, but it is not applicable to PPG data.Meanwhile, some research investigated the use of PRV as an alternative to HRV, and indicated that it would be sufficiently accurate to use PRV as an estimation of HRV for healthy people (and mostly younger) at rest condition [37]; However, physical activity and mental stressors may impair the agreement of PRV and HRV measures.
It is a complex process to calculate HRV parameters using a variety of procedures, such as the computation of beat-to-beat intervals, missing data estimation, data resampling, segmentation, and spectral estimation.In addition, each of these computational procedures can be implemented using different approaches, and therefore the combination of these techniques produces uncertainties in the results of HRV analysis.We note that the selection of optimal techniques for the computation of HRV has a strong relation to the data source.For example, ECG signals acquired in primary care usually have short-time durations, therefore, they are suitable to calculate short-term or ultra-short-term HRV parameters; In comparison, when using wearables for continuous monitoring, the standard 5-min data segments or long-term HRV can be derived from the measurements.
We note that other than age, gender, and BMI, there are also many factors that could potentially affect the analysis of HRV data, such as diet and environmental conditions [197].Meanwhile, it is known that the consumption of caffeinated drinks stimulates the ANS system and impacts cardiovascular activities [198].In terms of the impact of smoking, previous research indicated that acute, chronic active, and passive smoking may generate marked disruptions in the normal ANS functioning, which can be characterised by an increased sympathetic drive, reduced HRV, and parasympathetic modulation [199], [200].There are also studies investigating the use of HRV parameters to evaluate the physical training of athletes [201], and the results indicated that time-domain HRV measures are more consistent than frequency-domain in describing the chronic cardiovascular autonomic adaptations; In particular, the mean and standard deviation of RR intervals are the most consistent measures [201].
We highlight that the different types of uncertainties have correlations in the analysis of HRV data.For example, MA may produce missing values in the data measurements, and they also would affect the accurate detection of the QRS complex.It is therefore proper guidelines should be proposed to investigate and quantify their impacts on the results of HRV analysis.Our discussed techniques for uncertainty quantification and sensitivity analysis provide promising approaches to identify the important factors and quantify their impact on HRV data.For instance, with a human-based electromechanical model [202], it is possible to identify important HRV measures and quantify uncertainties (e.g., coming from computational approaches) in the healthcare management of cardiac diseases.

VI. CONCLUSION
HRV has emerged as an important tool in healthcare, including the diagnosis of cardiovascular diseases, screening of diabetes mellitus, and monitoring of psychiatric disorders.This survey provides a comprehensive review of uncertainties that can impact HRV analysis and affect healthcare outcomes, including measurement noises in HRV data analysis, computational modelling of HRV data, and physiological factors.Then, we evaluate the technical solutions to address these uncertainties and provide potential suggestions accordingly.Finally, we discuss the possible approaches to quantify the impacts of these uncertainties.The review could be used as a reference for researchers and practitioners, providing insights into future research directions in HRV data analysis.

Fig. 1 .
Fig. 1.Illustration of deriving HRV data from ECG and PPG signals, (a) RR intervals in the ECG signal, and (b) PP intervals in the PPG signal.N.U.: Normalised units.

Fig. 2 .
Fig. 2. The overall diagram of reviewing uncertainties in the analysis of HRV data.

Fig. 3 .
Fig. 3. Coloured noise in ECG signals.(a) Original ECG signal, (b) ECG signal contaminated by pink noise, (c) ECG signal contaminated by brown noise, (d) PSD comparison of the two types of noises.

TABLE VII PHYSIOLOGICAL
FACTORS IN HRV ANALYSIS