## 2 Materials and Methods
### Measurement and Signal Processing  
An AESCULON ICG device (Osypka Medical, Berlin, Germany) [34] was used for simultaneous HR and SV measurements. The time variate thoracic bioimpedance was acquired using a set of four electrodes (each) on the thorax and neck (sampling frequency = 200 Hz, sampling period = 5 ms). The erythrocyte orientation changes during opening and closing of the aortic valves; the corresponding change in impedance with each beat enables SV estimation. The AESCULON ICG device uses one-lead electrocardiography (ECG). QRS waves in the ECG were identified and the time-domain index of the R-peak was recorded as index [n].  
\[\text{HR [n]}=\text{round}\{\ \frac{60\times 200}{\text{Index [n]}-\text{Index [n -1]}},\ 0\ \} \tag{1}\]  
The function round [X,0] represents where X is rounded to the specified number of digits (0). The number will round to the nearest integer. For instance, round[98.98, 0] = 99.  
Missing RR interval and SV series points occur during riding stages due to motion artifacts during impendence acquisition. If HR[n] > 200 or HR[n] < 2 ppm, an autoregression extrapolation using the first 10 points of n were used to calculate a reasonable R peak. Similarly, if SV[n] > 200 or SV[n] < 2 ppm, the first 10 points of n were selected for autoregression and extrapolation of the point to calculate a reasonable SV value. The calculation equation for each beat is SV \(\times\) HR = CO, and the standard deviation (SD) and mean of CO, HR and SV were also calculated for every stage. The maximum SV value indicates the ACP stage; SV begins to decrease from this stage, which is also the limit of left ventricular volume diastolic expansion. The results of time-domain were used to compare physiological responses between rest and incremental exercise. The phase of SVV might affect the correlation between HRV and SVV, using cross-correlation to test the correlation in SV and RR intervals. Fast Fourier transform analysis was used for frequency domain analysis of HR and SV, and a power spectrum was generated for analyzing the spectrum. The distribution of components in the power spectrum of HRV power is divided into LF (0.04-0.15 Hz) and HF (0.15-0.4 Hz); both are absolute values (ms2). On the basis of the results of our previous study, the distribution of components in the power spectrum of SVV power is divided into LF (0.035-0.13 Hz) and HF (0.13-0.28); both are absolute values too (ms2). HRV and SVV are normalized from 0 to 1 \(\times\) 100%, as is conventional in HRV and SVV analysis, as follows:  
\[\text{nHF}=\frac{\text{HF}}{\text{LF}+\text{HF}}\times 100\% \tag{2}\]  
\[\text{nLF}=\frac{\text{LF}}{\text{LF}+\text{HF}}\times 100\% \tag{3}\]  
The LabVIEW (LabVIEW 2020, National Instruments Corp., Austin, TX, USA) platform was used to develop programs for processing RIP signals. The original raw RIP signal has substantial noise; thus, complementary ensemble empirical mode decomposition (CEEMD) was used to decompose the RIP signal into intrinsic mode functions (IMFs) to filter this noise [35]. An IMF is defined as follows: (1) the number of local maxima and local minima must be equal to the number of zero-crossings or dissimilar by at most one, and (2) at any point, the value of the average envelope must be approximately equal to zero. The main IMF among various IMFs was found in order to calculate the BR. CEEMD was used to identify the IMFs from the raw data, and the main signal was extracted on the basis of the correlation between the main IMF and raw data.  
The next step is to time-shift the main IMF component in order to identify the peak of the RIP signal; these acquired peaks are used to calculate BR. The time difference between adjacent peaks was calculated as well as breath-to-interval for subsequent analysis. A flowchart of the RIP signal processing method is displayed in Figure 2.