A Novel Method for Baroreflex Sensitivity Estimation Using Modulated Gaussian Filter

The evaluation of baroreflex sensitivity (BRS) has proven to be critical for medical applications. The use of α indices by spectral methods has been the most popular approach to BRS estimation. Recently, an algorithm termed Gaussian average filtering decomposition (GAFD) has been proposed to serve the same purpose. GAFD adopts a three-layer tree structure similar to wavelet decomposition but is only constructed by Gaussian windows in different cutoff frequency. Its computation is more efficient than that of conventional spectral methods, and there is no need to specify any parameter. This research presents a novel approach, referred to as modulated Gaussian filter (modGauss) for BRS estimation. It has a more simplified structure than GAFD using only two bandpass filters of dedicated passbands, so that the three-level structure in GAFD is avoided. This strategy makes modGauss more efficient than GAFD in computation, while the advantages of GAFD are preserved. Both GAFD and modGauss are conducted extensively in the time domain, yet can achieve similar results to conventional spectral methods. In computational simulations, the EuroBavar dataset was used to assess the performance of the novel algorithm. The BRS values were calculated by four other methods (three spectral approaches and GAFD) for performance comparison. From a comparison using the Wilcoxon rank sum test, it was found that there was no statistically significant dissimilarity; instead, very good agreement using the intraclass correlation coefficient (ICC) was observed. The modGauss algorithm was also found to be the fastest in computation time and suitable for the long-term estimation of BRS. The novel algorithm, as described in this report, can be applied in medical equipment for real-time estimation of BRS in clinical settings.


Introduction
Baroreflex can be defined as a homeostatic process that controls blood pressure via the autonomic nervous system (ANS) negative feedback loop by modulating the total peripheral resistance and the heart rate (HR) [1]. The most common indicator used to assess the effect of baroreflex on heartbeat interval control is the baroreflex sensitivity (BRS) [2]. It has been suggested in the literature that BRS is a critical evaluation index in clinical settings. One study [3] showed that patients suffering from essential hypertension tended to have a lower BRS compared to normal subjects. A lower BRS was also found to be a direct indicator of the risk of sudden cardiac death after myocardial infarction [4]. In addition, BRS is considered to be a potential marker in the early detection of autonomic dysfunction in diabetic mellitus patients [5]. Patients with neurodegenerative disease are found to have computation is more efficient than for conventional spectral methods, and, most importantly, there is no need to specify any parameter. The present paper introduces a new method, called modulated Gaussian filter (abbreviated as modGauss in this article) for BRS estimation. It is constructed with a more simplified structure than GAFD using only two bandpass filters of dedicated passbands to separately decompose the evenly resampled SBP and IBI sequences into two dedicated bands. This novel strategy makes modGauss more efficient than GAFD in the computation, yet the advantages of GAFD are preserved.
The article is divided into the following sections: The proposed modGauss algorithm for BRS estimation is introduced in Section 2. This section also discusses the materials, conventional frequency-domain approaches, the GAFD algorithm for BRS estimation (in brief) and the statistical experiments performed in this study. The results of the analysis are provided in Section 3 along with corresponding discussion presented in Section 4. The conclusions of the article are summarized in Section 5.

Methods and Materials
This section is divided into five subsections. In the first subsection, the proposed algorithm is introduced. The following subsection briefly describes conventional spectral methods for BRS estimation. The estimation of BRS by the GAFD algorithm is concisely described in the third subsection. The fourth subsection describes the materials (including the signal resource and the computing environment) that were used for the computational experiments in this paper. The statistical testing methodologies used for comparing our developed algorithm with other methods are described in the final subsection.

BRS Estimation by modGauss
The bandpass filters for the proposed method are constructed from the amplitude modulation on the Gaussian window. The Gaussian window of length (2M + 1) centralized at m = 0 is set as follows [34] in which β stands for the scalar inversely related to the standard deviation of the Gaussian distribution. Under the condition β > 2.5 [34], the truncated points of the Gaussian window are small enough to be omitted, and the associated discrete Fourier transform (DFT) is, in approximation, equal to From Equation (2), the spectrum of the Gaussian window substantially represents a non-causal lowpass filter (LPF) centered at ω = 0 rad with no phase offset. To ensure that the energy remains unchanged, the Gaussian window have to be normalized as follows.
where w[·] is specified in Equation (1). One important issue of concerned with respect to the Gaussian window is the selection of β. Assume the −6 dB cutoff (where the decay is one half of the highest amplitude) is adopted and the −6 dB cutoff frequencies are located at ω = ±ω d rad, then we have W(±ω d ) = 0.5 in Equation (2). By some algebraic manipulations, we have where ω d = 2π · f d / f s , f d is the −6 dB cutoff frequency and f s is the sampling frequency of the signal to be analyzed. Both f d and f s are in units of Hz. The bandpass filter in modGauss is acquired from the modulated Gaussian window w a [·]. Assume the resulting spectrum retains almost the same shape as the unmodulated spectrum. Let the carrier frequency in amplitude modulation be denoted as f c (in unit Hz), then the carrier signal can be generated by The modGauss filter can then be derived as The spectrum of w g [·] is the shifted version of W(ω) which is rescaled in amplitude and is now centered around ω = +2π · f c / f s and −2π · f c / f s rad. If the scale β is selected properly, the split spectra will have no evident overlap. In this case, the generated w g [·] shown in Equation (6) is, in essence, a non-causal bandpass filter with the passband ranging Assume the signal to be analyzed is denoted as x[n], with 0 ≤ n ≤ N − 1. Because the modGauss filter is non-causal, the analyzed signal x[·] requires to be extended outside the endpoints. In this study, the extension is performed by the following formula After the signal extension, the filtering operation is conducted as follows The derived y[·] is the bandpass-filtered components buried in signal x[·].
Another issue of concern for the modGauss filter is the determination of passband frequency range. In the spectral methods for human BRS estimation, the commonly adopted low-frequency (LF) band is within 0.04-0.15 Hz and the high-frequency (HF) band is located in the range 0.15-0.4 Hz [32]. In this condition, the center of the dedicated frequency range is selected to be the carrier frequency f c , which is 0.095 Hz for the LF band and 0.275 Hz for the HF band, respectively. The deviated frequency from the center was previously denoted as f d . It is assumed that the −6 dB cutoff frequencies are located at the endpoints of both bands in modGauss filters. In this case, the value of f d is equal to the deviation from f c to the end frequency. Therefore, it is 0.055 Hz for the LF band and 0.125 Hz for the HF band.
With selection of the f c and f d values referred to, the proposed method is constructed from two modGauss filters, one for the LF band and the other for the HF band. The filtering operation is based on Equation (8), which is conducted fully in the time domain but can derive the frequency-domain information in the same way as for other spectral methods.
The final issue for the modGauss filter requiring consideration is the filter length, which is determined by the parameter M. Because the components buried in the raw signal are closely related to the number of extremes (including both the local maxima and the local minima) appearing in the signal, the proposed method adopts an approach suggested by [35] according to the following formula where κ is a scalar parameter ranging from 1.1 to 3, N e denotes the number of extremes, N represents the signal length and · refers to the operator rounding the value to the closest integer toward −∞. To select the feasible value of κ, the estimation of the BRS value is performed with different κ values (1,1, 2 and 3); the results are shown in Figure 1. It can be observed that the results are identical for each of the values of κ. The value 2 is in the can be observed that the results are identical for each of the values of . The value 2 i the middle of the selection range and is an integer, which means it can save storage sp and offer fast computation in embedded applications. For these reasons, a constant va of 2 is selected for the parameter . . It can be observed that the BRS v is identical for each value of . In this study, is constantly set to 2. The data points beyond whiskers are marked by red + symbol.
In addition, the second derivative test [36] is chosen for detecting the extremes. T is, the points whose first derivative is zero, and for which the second derivative is than or greater than zero, are categorized as the local extremes. To prevent obtaining implausible value of , the limit of the extreme numbers is defined using the equa below in which ⌈⋅⌉ refers to the operator, rounding the value to the closest integer toward Signal processing is performed when the number of extremes detected is greater t . To facilitate readers' understanding of the proposed method, the modGauss f generation procedure for both the LF and HF bands is illustrated in Figure 2 from the t domain perspective. . It can be observed that the BRS value is identical for each value of κ. In this study, κ is constantly set to 2. The data points beyond the whiskers are marked by red + symbol.
In addition, the second derivative test [36] is chosen for detecting the extremes. That is, the points whose first derivative is zero, and for which the second derivative is less than or greater than zero, are categorized as the local extremes. To prevent obtaining an implausible value of M, the limit of the extreme numbers is defined using the equation below in which · refers to the operator, rounding the value to the closest integer toward +∞. Signal processing is performed when the number of extremes detected N e is greater than N e TH . To facilitate readers' understanding of the proposed method, the modGauss filter generation procedure for both the LF and HF bands is illustrated in Figure 2 from the time domain perspective.
In this figure, the generation of both modGauss filters is based on Equations (1), (3)-(6) with M = 50, f s = 4 Hz, f c = 0.095 Hz (for the LF band) and 0.275 Hz (for the HF band), and f d = 0.055 Hz (for the LF band) and 0.125 Hz (for the HF band). The corresponding modGauss filters from the frequency-domain perspective are shown in Figure 3, where the parameters are the same as those for Figure 2. It can be observed that the generated modGauss filters are bell-shaped bandpass filters centered at 0.095 and 0.275 Hz with the −6 dB cutoff frequencies corresponding to the frequency ranges of the LF and HF bands, respectively.   Figure  3, where the parameters are the same as those for Figure 2. It can be observed that the generated modGauss filters are bell-shaped bandpass filters centered at 0.095 and 0.275 Hz with the −6 dB cutoff frequencies corresponding to the frequency ranges of the LF and HF bands, respectively.
As shown in Figure 3, it can be observed that both the modGauss filters are symmetrical with a center at index 0, and the computation of Equation (8)     In this study, the SBP sequence expressed in mmHg, and its respective IBI sequence expressed in ms, are resampled using a 4 Hz sampling rate by cubic spline interpolation for BRS estimation. This is the primary reason why the parameter is set to 4 Hz in Fig-Figure 3. Generation procedure for the modGauss filters in LF and HF bands (from frequencydomain perspective). As shown in Figure 3, it can be observed that both the modGauss filters are symmetrical with a center at index 0, and the computation of Equation (8) is inherently the convolution sum of w g [·] and x e [·]. From the frequency-domain perspective, this is a direct multiplication of the frequency response of w g [·] and the spectrum of x e [·], and only the components within the dedicated frequency ranges are reserved.
In this study, the SBP sequence expressed in mmHg, and its respective IBI sequence expressed in ms, are resampled using a 4 Hz sampling rate by cubic spline interpolation for BRS estimation. This is the primary reason why the parameter f s is set to 4 Hz in Figures 2 and 3. Let the separated LF and HF component of the interpolated SBP sequence be designated as x SBP,LF [n] and x SBP,HF [n], for 0 ≤ n ≤ N − 1. The LF and HF band energies of SBP sequence are derived by Let x IBI,LF [n] and x IBI,HF [n], with 0 ≤ n ≤ N − 1, denote the extracted LF and HF components for the interpolated IBI sequence. The LF and HF band energies for IBI sequence are acquired by The BRS is then estimated using the so-called α index [37] in the LF band and the HF band, as well as in the averaged form, by where the expression unit of all indices is ms/mmHg. The computation procedure of our method is outlined in the following steps.
Step 1: The SBP sequence, expressed in mmHg, as well as its respective IBI sequence expressed in ms, are evenly resampled by f s Hz, which is a 4-Hz cubic spline interpolation in this study.
Step 2: Find the number of extremes for both the interpolated sequences. Then, determine the parameters M and β according to Equation (9) and Equation (4), respectively. In this study, the second derivative test [36] is applied for detection of extremes. In addition, the value of f d is 0.055 Hz for the LF band and 0.125 Hz for the HF band. If the number of extremes is less than N e TH , which is derived according to Equation (10), this implies the interpolated sequence is nearly a monotonic trend wave and the energies in both the LF and HF bands can be considered to be zero. The BRS estimation is only conducted once the number of detected extremes is greater than N e TH .
Step 3: Generate the carrier signals for both the LF and HF bands according to Equation (5), where the value of f c is 0.095 Hz for the LF band and for the HF band is 0.275 Hz.
Step 4: Extend the interpolated SBP and IBI sequences beyond the endpoints according to Equation (7), and then conduct the filtering process according to Equation (8) for both the LF and HF bands separately.
Step 5: Derive the energies of the interpolated SBP and IBI sequences in both the LF and HF bands according to Equation (11) and Equation (12), respectively. Finally, the BRS values are estimated according to Equation (13).

BRS Estimation with Spectral Methods
Three conventional spectral approaches are adopted for performance comparison in this study. One is the Welch's periodogram [38], whereas the others are the autoregressive (AR) spectrum estimation [39] and the wavelet [40]. The SBP, as well as the IBI sequences, should be resampled uniformly by f s Hz, which is a 4-Hz cubic spline interpolation in this study, before conducting the spectral analysis.
For the Welch's periodogram, the parameters of the 512-point Hanning window and 512-point fast Fourier transform (FFT), along with 50% overlap among adjacent windows, are selected for performance comparison. The energies of the LF band and the HF band are deduced by direct addition of the periodograms in the respective frequency range. Let the derived energies in LF and HF bands for the interpolated SBP and IBI sequences be represented by E SBP,LF , E SBP,HF , E IBI,LF and E IBI,HF , respectively. The BRS values are then estimated by Equation (13).
Another commonly used method for spectrum analysis is the AR model [39]. It consists of an all-pole approach whose transfer function is given by where σ is the scalar related to the estimation error, a i (1 ≤ i ≤ P) are the AR coefficients and p i (1 ≤ i ≤ P) denote the poles of the model. The order of the AR model is represented by the parameter P. The AR model of fixed order 25 is selected for both the IBI sequences and an interpolated SBP, and Burg's maximum extropy method is used to derive the estimation error σ as well as the AR coefficients. This method has the advantage of high resolution, even for short data records, and model stability [41]. The poles of the model, p i (1 ≤ i ≤ P), can be derived directly by factorization after the AR coefficients a i (1 ≤ i ≤ P) have been acquired. The resulting spectrum obtained from the acquired AR model is defined as Let P d ( f k ) denote the energy contributed from the dominant pole p k (the pole closest to the unit circle) corresponding to the frequency f k (in unit Hz) for a certain band of the AR spectrum, then P d ( f k ) can be derived as follow [42].
where f i represents the corresponding frequency (in unit Hz) of the pole p i (1 ≤ i ≤ P).  (13). The advantage of spectral analysis by the scalogram of complex Morlet wavelet is the ability to focus on analysis within the dedicated frequency range with a selected frequency resolution [40]. In this study, the frequency resolution is selected to be 0.005 Hz. After having acquired the scalogram under pre-determined scales (which are the scaled reciprocal of frequencies) in different time segment, the averaged spectrum is then time-averaged along each frequency. The energy within the HF band (0.15-0.4 Hz) and the LF band (0.04-0.15 Hz) for the SBP and IBI can then be estimated by Equation (13).

BRS Estimation by GAFD
GAFD [33] is a novel approach for BRS estimation which adopts a three-level structure similar to wavelet decomposition but using Gaussian lowpass filters of three different cutoff frequencies (which are 0.04, 0.15 and 0.4 Hz, respectively). The SBP sequence in mmHg and its associated IBI sequences in ms were resampled evenly in 4 Hz by cubic spline interpolation. Before the signal decomposition, the Gaussian lowpass filters for GAFD are based on Equation (3), and the signal must also be extended out of the endpoints according to Equation (7). The parameters M and β are determined by Equations (9) and (4), with the values of f d being 0.04, 0.15, and 0.4 Hz, respectively, for GAFD. The filtering operation is conducted by Equation (8), with the filter w g [·] being replaced by w a [·], as shown in Equation (3). In the following content, the symbol x[·] is used to represent the resampled SBP or IBI sequences. The procedure for BRS estimation by GAFD is briefly summarized below. A detailed introduction is provided in [33].
Step 1: Perform the lowpass filtering to divide the sequences x[·] to components of above 0.4 Hz and below 0.4 Hz. Let the components below 0.4 Hz be referred to as x 04 [·].
Step 3: Perform the lowpass filtering on x 015 [·] to divide the sequences to components above 0.04 Hz and below 0.04 Hz. The components below 0.04 Hz represent the LF components for the estimation of BRS and are designated as x LF [·]. Let the LF components for IBI and SBP sequences be expressed by x IBI,LF [·] and x SBP,LF [·]. In a similar way to the HF band, the LF band energies can also be calculated by squared summation, as shown in Equations (11) and (12).
Step 4: Lastly, the values of BRS are approximated with the α index in the LF band and the HF band, along with the averaged form, independently, as shown in Equation (13).
The MATLAB codes for the BRS estimation by spectral approaches and GAFD referred to in this article can be found on the website: https://figshare.com/articles/software/ Sample_codes_for_the_assessment_of_baroreflex_sensitivity_BRS_/12820430 (accessed on 10 June 2022).

Materials
For the purpose of this research, both beat-by-beat SBP and IBI sequences were obtained from the EuroBavar dataset [43]. The dataset, which has been used in previous studies [33,44,45] to compare the performance of different methods used for BRS estimation, can be downloaded online. The EuroBavar dataset contains 46 recording files acquired from 21 subjects composed of 4 males and 17 females. The participating subjects are organized as follow: • 12 normotensive outpatients, including two patients who were treated for hypercholesterolemia, one diabetic patient with no cardiac neuropathy and one woman who was three months pregnant, • two patients who were treated for hypertension, • one non-treated hypertensive patient, • two patients suffering from cardiac autonomic failure (one with heart transplantation and one patient with diabetic neuropathy) • four healthy volunteers.
ECG signals (Cardiocap™ II, Datex-Ohmeda, Helsinki, Finland) and continuous noninvasive blood pressure (Finapres ® , Finapres Medical Systems, Enschede, The Netherlands) were recorded on all subjects at a sampling frequency of 500 Hz by 16-bit analog-to-digital conversion resolution in the standing position (files marked with the penultimate letter "S") and supine position (files marked with the penultimate letter "L") in a silent room according to a uniform procedure. Informed consent was obtained for all participants prior to collecting the EuroBavar dataset. Furthermore, the ECG and blood pressure signals from this dataset are available as waveforms (files with the ending letter "C") and beat-by-beat sequences (files with the ending letter "B"). For all the computer experiments of the study, the beat-by-beat SBP, along with the IBI sequences, were adopted.
All the statistical analyses and computational experiments in this study were performed on a MacBook Pro (mid 2015, with 2.5 GHz Intel Core i7 and 16 GB 1600 MHz DDR3 memory, Apple Inc., Cupertino, CA, USA) using MATLAB ® (2020b, MathWorks, Inc., Natick, MA, USA).

Statistical Analysis
The values of the BRS values were derived from Equation (13) by Welch's periodogram, then by GAFD and finally by modGauss. Various statistical tests were applied to assess the effectiveness of the suggested modGauss method. Firstly, the Shapiro-Wilk normality test [46] was performed to check if the BRS values in the whole dataset, as well as in the dataset for the supine and standing positions, were normally distributed or not. We then used non-parametric statistical methods for testing; the Wilcoxon signed-rank test [47] was used to assess if there was a statistical difference in the BRS values between the standing and supine positions for each approach. The Wilcoxon rank sum test [48] was employed to assess whether any statistical difference existed between the modGauss method and the other four methodologies. Lastly, to assess the level of correlation between modGauss and the others methodologies, the intraclass correlation coefficient (ICC) [49,50] was determined. A MATLAB algorithm, developed by Salarian [51] and based on [49], was adopted for the estimation of the 95% confident intervals (CI) and the ICC. Six different types of ICC computation can be used for the ICC analysis; however, in the context of this study, the ICC for a two-way mixed approach and single score (type "A-1" in the code) was chosen for the computation of the ICC; the lower and upper bounds of the ICC were calculated assuming a significance level of 0.05.

Experimental Results in Time Domain and Frequency Domain
This article proposes a new method called modGauss for the estimation of BRS using α indices according to Equation (13). Figures 2 and 3 have previously demonstrated the generation procedure of modGauss filters from time-and frequency-domain perspectives, respectively. To demonstrate the effect of the modGauss filters, the decomposed results for one data record in the EuroBavar dataset (the record numbered B002LB) are shown in Figures 4 and 5 (for the resampled SBP sequence); Figure 4 presents the results from a time-domain perspective, while Figure 5 shows the corresponding results from a frequencydomain perspective.
The decomposed results for one data record of the EuroBavar dataset (B002LB) are presented in Figures 6 and 7 (for the resampled IBI sequence). Figure 6 shows the results in the time domain, while Figure 7 describes the corresponding results in the frequency domain.
The Welch's periodogram (with 512-point FFT, 512-point Hanning window and a 50% overlap between the neighboring windows) was used to acquire the spectra for the resampled sequences and the decomposed signals in Figures 5 and 7. For the figures in the time domain, it can be observed that the raw signals were separated into the lowerfrequency components and the higher-frequency components (refer to Figures 4 and 6). However, it is very evident that the spectra of the lower-frequency and higher-frequency components have been separated to the discrete ranges of 0.04-0.15 Hz and 0.15-0.4 Hz by the modGauss filters from the corresponding frequency-domain perspective (as shown in Figures 5 and 7). generation procedure of modGauss filters from time-and frequency-domain perspectives, respectively. To demonstrate the effect of the modGauss filters, the decomposed results for one data record in the EuroBavar dataset (the record numbered B002LB) are shown in Figures 4 and 5 (for the resampled SBP sequence); Figure 4 presents the results from a time-domain perspective, while Figure 5 shows the corresponding results from a frequency-domain perspective.  respectively. To demonstrate the effect of the modGauss filters, the decomposed results for one data record in the EuroBavar dataset (the record numbered B002LB) are shown in Figures 4 and 5 (for the resampled SBP sequence); Figure 4 presents the results from a time-domain perspective, while Figure 5 shows the corresponding results from a frequency-domain perspective.   The decomposed results for one data record of the EuroBavar dataset (B002LB) are presented in Figures 6 and 7 (for the resampled IBI sequence). Figure 6 shows the results in the time domain, while Figure 7 describes the corresponding results in the frequency domain.  The decomposed results for one data record of the EuroBavar dataset (B002LB) are presented in Figures 6 and 7 (for the resampled IBI sequence). Figure 6 shows the results in the time domain, while Figure 7 describes the corresponding results in the frequency domain.

Performance Evaluation of the modGauss
To evaluate the performance and effectiveness of our proposed method, four other approaches-GAFD, AR spectrum, Welch's periodogram and wavelet-were chosen for performance comparison. The present study used the EuroBavar dataset to confirm the applicability of the newly suggested approach. Although this dataset has been used in other studies [33,44,45], it is important to note that few studies have shown all the derived BRS values for each data record in the dataset. Figure 8 represents the BRS value for each of the methods discussed in this study. It can be observed that the median BRS value was almost identical for each of the five methods, which indicates that the proposed modGauss algorithm is comparable to the other methodologies for BRS estimation.
by the modGauss filters from the corresponding frequency-domain perspective (as sho in Figures 5 and 7).

Performance Evaluation of the modGauss
To evaluate the performance and effectiveness of our proposed method, four o approaches-GAFD, AR spectrum, Welch's periodogram and wavelet-were chosen performance comparison. The present study used the EuroBavar dataset to confirm applicability of the newly suggested approach. Although this dataset has been used other studies [33,44,45], it is important to note that few studies have shown all the deri BRS values for each data record in the dataset. Figure 8 represents the BRS value for e of the methods discussed in this study. It can be observed that the median BRS value almost identical for each of the five methods, which indicates that the proposed m Gauss algorithm is comparable to the other methodologies for BRS estimation. After estimation of the BRS values for each of the methods, detailed statistical va were obtained using different tests, including the Shapiro-Wilk normality test, the W coxon signed-rank test, the Wilcoxon rank sum test and the intraclass correlation co cient (ICC). First, the Shapiro-Wilk normal test [46] was applied to assess if the BRS va for the whole dataset, as well as for the datasets for supine and standing position, w normally distributed. In Table 1 below, the test results for the five methods, includ After estimation of the BRS values for each of the methods, detailed statistical values were obtained using different tests, including the Shapiro-Wilk normality test, the Wilcoxon signed-rank test, the Wilcoxon rank sum test and the intraclass correlation coefficient (ICC). First, the Shapiro-Wilk normal test [46] was applied to assess if the BRS values for the whole dataset, as well as for the datasets for supine and standing position, were normally distributed. In Table 1 below, the test results for the five methods, including modGauss, GAFD, AR, Welch and Wavelet in different conditions (whole dataset, supine and standing position) are presented.
It was observed that all the p values were smaller than 0.05, indicating that there was no Gaussian distribution of BRS values for each test. Based on this, the non-parametric Wilcoxon signed-rank test [47] was utilized to assess if there was a statistical difference in the BRS values between the standing and supine positions for each of the four methods. The experimental results for this test are shown in Table 2. "Whole" refers to the analysis of the entire dataset; "Supine" stands for the analysis of the data in supine position; "Standing" stands for the analysis of the data in standing position.  Since all the p values were lower than 0.05, this implies that there were statistical differences in the BRS values for the supine and standing positions for all the five approaches. Therefore, a further test was needed to determine if there was any statistically significant difference between the proposed modGauss and the other four methods for BRS estimation. The Wilcoxon rank sum test [48] was then to this end; Table 3 shows the comparison results under each condition (supine and standing position, whole dataset). the p value of BRS α index in HF band under different "Test Condition"; Test Condition: "Whole" refers to the analysis of the entire dataset; "Supine" stands for the analysis of the data in supine position; "Standing" stands for the analysis of the data in standing position.
It was observed that all the p values were higher than 0.05, which implies that there was no substantial difference between the modGauss and the other four BRS estimation approaches. The ICC [49,50] was used to evaluate the level of correlation between the BRS values deriving from modGauss and the other four methodologies. The results of the ICC analysis comparing modGauss with the other approaches under different conditions (supine and standing position, whole dataset) are shown in Table 4. The 95% confidence interval (CI) is also provided in this table for each comparison, along with the ICC value. Even considering the lowest CI (which is 0.91 for α LF compared between the modGauss and GAFD under the supine condition), the obtained results indicated "excellent" agreement (based on the guideline described in [52]) for the estimation of BRS between the modGauss and the other four approaches. Whole" refers to the analysis of the entire dataset; "Supine" stands for the analysis of the data in supine position; "Standing" stands for the analysis of the data in standing position.

Discussion
The various experiments and tests performed using the proposed modGauss algorithm have shown that the algorithm is highly comparable to other conventional BRS estimation methods. More importantly, a major advantage of our proposed modGauss method is its simplicity for performing complex computation. As its highly simplified structure uses only two dedicated bandpass filters and bypasses the three-level structure of the GAFD algorithm [33], this strategy is more computationally efficient than the GAFD method while retaining all the advantages of the latter. Another major advantage of the proposed modGauss method is its fast computation time. A fast-computational time is an important criterion in defining the performance of computational algorithms. In this study, a comparative analysis of the proposed modGauss method was performed to evaluate the computational time of the algorithm using the EuroBavar dataset. Four different methods (GAFD, AR, Welch and Wavelet) were chosen for comparison and the results are presented in Table 5. It can be observed from this table that the modGauss algorithm had a shorter computational time in the 0.05 quartile, median and 0.95 quartile. The p values of the Wilcoxon rank sum one-tailed tests for modGauss versus the other four methods were all less than 0.05. This means that the modGauss algorithm was faster in computation, which makes it more attractive for performing rapid computation for BRS estimation. Although the elapsed time for each method was of the order of milli-seconds, it should be noted that the values were evaluated under the platform of a mid-2015 MacBook Pro. If the computation were conducted in the embedded micro-controller of a medical device, the elapsed time might be very different, especially for the spectral methods because of their relatively high demand for computational resources. When selecting a computational algorithm for BRS estimation, a key consideration is the ability of the algorithm to perform long-term analysis. In this regard, a comparative analysis was performed for long-term BRS estimation between the proposed modGauss and sleep stage, with the results shown in Figure 9. The signal slp45 from the MIT-BIH Polysomnographic Database [53] of PhysioNet [54] was adopted for the long-term BRS analysis. In this figure, the proposed modGauss algorithm was used to perform an analysis of the BRS estimate for a duration of 6 h. From Figure 9, it can be seen that the BRS values generally increased after the awake stages, which implies arousal of the baroreflex function when a person is awoken from the sleep stage. With respect to sleep, the increase in BRS values and the delay in BRS elevation after waking up are topics which may warrant further study.
Moreover, when compared to the other four methods reviewed in this study, the long-term analysis of the modGauss algorithm suggest that it performed as well as the other methods, as shown in Figure 10. This, combined with the previous observations on elapsed time (Table 5), leads to the conclusion that, overall, the modGauss algorithm is better and more suitable for the long-term estimation of BRS.
Although the new modGauss method proposed in this paper has various advantages, limitations may include the double reflection that is required at the start and end points of the filter to perform signal filtering. In practical use, a short waiting time is required to assure the signal entirely entering the modGauss filtering window. However, if the results obtained in this waiting window can be omitted without influencing the clinical decision, this limitation does not affect the performance of the algorithm because the modGauss algorithm only uses convolution, which results in faster computation time, as shown in Table 5. The methodologies, such as AR, Welch, and wavelet require more complex computation for model estimation, complex multiplication or to be spread out in the time and frequency domains, which causes them to have longer computation times because of the relatively higher demand on computing resources. The novel modGauss algorithm as discussed in this paper, provides new insights into computational methods for BRS estimation. Future development of this work could lead to the incorporation of the modGauss algorithm into embedded medical devices for real-time BRS analysis. Real-time estimation of BRS is important in clinical application as it provides critical information for continuous control of patient blood pressure [55,56]. Despite the existence of other methods for real-time BRS estimation [57,58], these methods typically use advanced microcontrollers for computation, making them more demanding on resources and more expensive. The modGauss, which is a based on a linear convolutional approach, does not require significant computational resources or advanced microcontroller systems. values generally increased after the awake stages, which implies arousal of the baroreflex function when a person is awoken from the sleep stage. With respect to sleep, the increase in BRS values and the delay in BRS elevation after waking up are topics which may warrant further study. Moreover, when compared to the other four methods reviewed in this study, the long-term analysis of the modGauss algorithm suggest that it performed as well as the other methods, as shown in Figure 10. This, combined with the previous observations on elapsed time (Table 5), leads to the conclusion that, overall, the modGauss algorithm is better and more suitable for the long-term estimation of BRS. Although the new modGauss method proposed in this paper has various advantages, limitations may include the double reflection that is required at the start and end points of the filter to perform signal filtering. In practical use, a short waiting time is For BRS estimation, the phase spectrum is sometimes adopted in spectral approaches to depict the timing difference from the sympathetic and vagal nerves in different conditions [29]. In addition, some spectral approaches utilize the threshold of coherence (e.g., greater than 0.5) between SBP and IBI spectral estimates to determine the bands used for the computation of BRS [45]. The lack of phase and coherence information is another limitation of the presented modGauss algorithm because its computation is conducted thoroughly in the time domain. For this issue, correlation analysis between the separated SBP and IBI sequences in the LF and HF bands may be helpful to compensate for such constraints. Due to the differences in the selected band and in the spectrum analysis method, there have been eleven different spectral approaches for BRS estimation proposed [45]. Because there has been no gold standard method for BRS estimation until now, this study selected three representative spectral analysis methods for performance comparison.
BRS estimation has been shown to be crucial in clinical settings and the evaluation of BRS has been applied in cardiovascular disease (CVD) [3,4,9,10,18], diabetes [5], neurodegenerative disease [6], diagnosis of brain stem death [7], OSAS [8], surgery [16], physiological function evaluation while conducting physical and psychological tasks [11], CKD [17], and pain perception [19,20]. No dedicated medical device has been available for BRS estimation until now. One attractive approach is to embed the BRS estimation algorithm in contemporary medical equipment. Accumulated observations from BRS estimation in various situations may provide new insights into the role that baroreflex function plays in clinical settings. The proposed algorithm is a strong candidate for practical implementation because of its computational simplicity and efficiency.

Conclusions
A new algorithm named modGauss, developed for the estimation of BRS values, is described in this paper. The proposed method utilizes two bandpass filters of dedicated passbands to separately decompose the evenly resampled SBP and IBI sequences into LF and HF bands. Convolution and squared summation are used without the need to specify a parameter, thus rendering the developed algorithm highly efficient in computing. The EuroBavar dataset was used to perform computational experiments and the BRS estimation results were compared with four different methods, namely the Welch's periodogram, GAFD, the AR spectrum and wavelet. It was found that the BRS values under different conditions (whole dataset, supine position and standing position) obtained through modGauss showed "excellent" concordance compared to the other four methodologies. In addition, further tests showed that the proposed modGauss algorithm, which has a very simplified structure, had faster computation times compared with the other methodologies (Table 5) and was also suitable for the long-term analysis of SRB estimation (Figures 9 and 10).