## 2 Materials and Methods
### Feature Extraction
#### 2.3.1 Empirical Mode Decomposition  
The EMD method is an adaptive data analysis method applicable to nonlinear and non-stationary signals [23]. The EMD model decomposes data into finite intrinsic mode functions (IMFs) based on directly extracting the energy associated with various intrinsic time scales. With the Hilbert transform, the IMFs yield instantaneous frequencies as functions of time that give sharp identifications of embedded structures. The final presentation of the results is an energy-frequency-time distribution, designated as the Hilbert spectrum, which accurately represents non-stationary and nonlinear signals [23].  
Unlike Fourier analysis or wavelet transforms, EMD breaks down a time signal into a set of base signals derived from the data itself; this unique approach allows EMD to preserve the full non-stationarity of the signal [24]. In addition to contrasting the Fourier analysis that produces a series of sine and cosine functions of fixed amplitudes to represent each frequency constituent in the signal, the IMFs are oscillatory modes whose amplitude and frequency vary over time [23].  
The EMD of a time-series signal x(t) can be represented as a sum of IMFs, \(\rm IMF_{i}(t)\), and a residue component, \(\rm r_{N}(t)\), as represented in Equation (1):  
\[\rm x(t)=\sum_{i=1}^{N}IMF_{i}(t)+\ r_{N}(t) \tag{1}\]  
The algorithm to extract the IMFs of a time-series signal using EMD follows an iterative process known as the sifting process, summarized as follows [23]:  
1. Identify the signal maxima and minima;
2. Compute the interpolated upper and lower envelopes and the instantaneous local mean of the envelopes;
3. Subtract the obtained local mean from the original signal x(t) to obtain the first component \(\rm IMF_{i}(t)\);
4. Check whether the component \(\rm IMF_{i}(t)\) satisfies the two basic conditions of the IMF: 1. The number of extrema--maxima and minima--and the number of zero-crossings in a signal should be either equal or differ by a maximum of one; 2. At any point, the mean value of two envelopes, one formed by connecting local maxima and the other by local minima, should be zero.
5. Repeat steps (i)-(iii) until it satisfies the conditions of the IMF (or by applying a stopping criterion such as the number of repetitions);
6. Repeat steps (i)-(iii) again for the calculation of the next IMFs, until no more components can be extracted (or by removing criteria such as the number of required IMFs).  
After the signal decomposition, the Hilbert transformation applied to the obtained IMFs provides an analytical signal representation (ASR) of IMFs (i.e., a complex-valued function with no negative frequency components). Any complex signal, z(t), can be considered the sum of its real part, \(\rm IMF_{i}(t)\), and its imaginary part, \(\rm Im(t)\), and rewritten in a polar coordinate system, as demonstrated in Equation (2). Equations (3) and (4) denote the instantaneous amplitude A(t) and phase \(\rm\theta(t)\) of the complex analytic signal, respectively:  
\[\rm z(t)=IMF_{i}(t)+jIm(t)=A(t)\ e^{i\rm\theta(t)} \tag{2}\]  
\[\rm A(t)=\sqrt{IMF_{i}(t)^{2}+\ Im(t)^{2}} \tag{3}\]  
\[\rm\theta(t)=tan^{-1}(\frac{Im(t)}{IMF_{i}(t)}) \tag{4}\]  
When applied to HRV analysis, the RR interval is decomposed into amplitude and frequency-modulated (AM-FM) signal components, the IMFs. The study in [41] demonstrated that the isolation of the first four components of the EMD is necessary to recognize the spectral bands of the autonomic modulation. Thus, this study limited the decomposition to four IMFs, considering they held the most significant signal variation.