Exploring the Intrinsic Features of EEG Signals via Empirical Mode Decomposition for Depression Recognition

Depression is a severe psychiatric illness that causes emotional and cognitive impairment and has a considerable impact on patients’ thoughts, behaviors, feelings and well-being. Moreover, methods for recognizing and treating depression are lacking in clinical practice. Electroencephalogram (EEG) signals, which objectively reflect the internal workings of the brain, is a promising and objective tool for recognizing and diagnosing of depression and enhancing clinical effects. However, previous EEG feature extraction methods have not performed well when exploring the intrinsic characteristics of highly complex and nonstationary EEG signals. To address this issue, we propose a regularization parameter-based improved intrinsic feature extraction method of EEG signals via empirical mode decomposition (EMD), which mines the intrinsic patterns in EEG signals, for depression recognition. Furthermore, our method can effectively solve the problem that EMD fails to extract intrinsic features. In this method, we first select an appropriate regularization parameter to generate the regularization matrix. Next, we calculate the sum of the matrix products of the IMFs and the regularization matrix and leverage the inverse of this matrix to extract the intrinsic features. The classification results of our method on four EEG datasets reached 0.8750, 0.8850, 0.8485 and 0.7768, respectively. In addition, compared with the iEMD method, our method requires less computational costs. These results support our claim that our method can effectively strengthen the depression recognition performance, and our method outperforms state-of-the-art feature extraction approaches.


I. INTRODUCTION
D EPRESSION is a severe psychiatric disorder that causes emotional and cognitive impairment and affects the thoughts, behaviors, feelings and well-being of patients, impairing their ability to live and work. According to recent reports from the World Health Organization (WHO), depression affects approximately 280 million people of all ages worldwide and cause more than 700,000 people die from suicide each year [1], [2]. Thus, it has become a serious threat to human life. Clinical diagnoses of depression, which are critical for receiving effective and timely treatment, rely on self-assessment scales and physician interviews, and depression is diagnosed by psychiatrists subjectively based on their experience. However, in clinical practice, experienced psychiatrists often have limited resources and the number of psychiatrists and patients is extremely unbalanced. Unfortunately, there are no objective criteria for assessing depression in clinical practice [3]. Therefore, it is crucial to develop an objective and valid method for recognizing and evaluating depression to understand the effect of different treatments.
Depression can be effectively identified with psychophysiological data-driven methods [4], [5], [6], [7], [8]. Specifically, physiological data can be used to develop objective and accurate tools for determining the physiological and psychiatric states of subjects. Electroencephalogram (EEG) is a widely recognized physiological signal that reflects the internal workings of the brain [9], [10]. EEG signals captured from the scalp are safe, noninvasive, low-consumption and easy to acquire. Thus, EEG signals can be used to develop a promising and objective tool for recognizing and diagnosing depression to enhance clinical outcomes. Several studies have confirmed that EEG signals can be used to effectively recognize depression [11], [12]. Furthermore, EEG signals in the frontal lobe and bilateral temporal region can be used as biomarkers for distinguishing depressed patients from healthy controls [13], [14], [15], [16].
Effective features are crucial in EEG recognition tasks. If a feature is invalid, a classifier can have difficulty identifying This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ group differences in an EEG signal. Existing EEG feature extraction methods in depression recognition research mainly utilize the following features. (1) Frequency domain feature methods use fast Fourier transformations (FFTs) [17] to convert EEG data from the time domain to the frequency domain to extract features such as the power spectral density, maximum power, and centroid frequency [18], [19], [20]. (2) Temporal domain features, which show changes in EEG signals over time, are implemented in the current study by using the autoregressive (AR) model [21] and recurrent neural networks (RNNs) [22], [23]. (3) Spatial domain feature methods utilize the original signals to generate a covariance matrix according to the EEG channel space and then construct spatial filters to learn the spatial patterns as features. The most impressive spatial domain feature method is the common spatial pattern (CSP) approach [24], [25], which extracts EEG features by exploring the spatial distribution components of different tasks. (4) Nonlinear features can be used to investigate the complexity and chaos of EEG signals. Traditional nonlinear features include the wavelet entropy [26], C0 complexity [27], Renyi entropy [28], permutation entropy [29], Kolmogorov entropy [30], and Lempel-Ziv (LZ) complexity [31]. (5) Image feature methods convert EEG signals into time-series images, spectral images, or topography maps and extract features from these images. Previous research has used convolutional neural networks (CNNs) [32], [33], [34] to mine the statistical characteristics of EEG signals in these images. Although these features are widely used in depression recognition research, they do not perform well when exploring the intrinsic characteristics of highly complex, nonstationary EEG signals.
Empirical mode decomposition (EMD) is an adaptive time-series signal decomposition method proposed by N.E. Huang et al. for complex and nonstationary signals [35]. EMD can decompose a complex signal into a finite number of simple signals known as intrinsic mode functions (IMFs), which reflect the essential physical characteristics of the signal nature. Therefore, EMD is an effective technical tool for investigating the intrinsic patterns and valid features of EEG signals [36], [37]. EMD-based feature extraction methods typically use IMFs to represent EEG signals and use the expansion coefficients as the features [4], [35]. However, due to mode mixing, the inverse of the matrix product of the IMFs and its transpose matrix does not exist. Thus, EMD cannot extract features of EEG signals effectively and accurately. To address this problem, a regularization parameter-based improved intrinsic feature extraction method of EEG signals via EMD is proposed in this work. In this method, we introduce a regularization parameter to effectively optimize the intrinsic features of EEG signals. First, the improved method multiplies the identity matrix by an appropriate regularization parameter to generate the regularization matrix. Next, we calculate the sum matrix by adding the matrix product of the IMFs and its transpose matrix to the regularization matrix. Then, we compute the inverse of this sum matrix. Since the regularization matrix is nonsingular, the sum matrix is also a nonsingular matrix. Thus, the inverse of the sum matrix can be used in place of the inverse of the matrix product of the IMFs in the subsequent EMD-based feature extraction to efficiently optimize the intrinsic features of EEG signals. Therefore, our improved intrinsic feature extraction method can effectively address the issue of IMF mode mixing in the traditional EMD method, which causes issues for EEG feature extraction. Consequently, our method can effectively extract the features by calculating the expansion coefficients of all IMFs. Furthermore, we demonstrate the performance of our EMD method in exploring EEG intrinsic characteristics and extracting effective features.
The main contributions of this work are as follows: 1) We provide an improved intrinsic feature extraction method that effectively enhances the robustness of EMD-based EEG feature extraction and supports effective depression recognition. 2) Our method effectively explores the intrinsic characteristics of EEG signals and improves the interpretability of these features according to physics and mathematics, thus enhancing the relationship between neuroscience and physiological electrical signal studies. 3) We illustrate a framework for intrinsic feature optimization and provide a novel idea and approach for recognizing depression, which will inspire innovative healthcare applications. The remainder of this paper is organized as follows. In Section II, the materials and theoretical basis are introduced. Then, our regularization parameter-based improved intrinsic feature extraction method of EEG signals via EMD for depression recognition is proposed in Section III. An analysis and comparison of the experimental results are demonstrated in Section IV. Finally, we discuss and conclude this work in Section V.

A. Data Acquisition and Preprocessing
The EEG signals used in this work were acquired from three spontaneous resting-state EEG datasets, including 128channel, 64-channel, and 3-channel EEG signals, and one auditory stimulus-evoked 3-channel EEG dataset. The 64-channel and 3-channel resting-state signals and the 3-channel auditory stimulus-evoked EEG signals were acquired independently in this study. The 128-channel resting-state EEG signals used in this paper were acquired from MODMA [38].
During the data acquisition process, there were no significant group differences in terms of the age and sex of the subjects in each dataset. All subjects were right-handed, had normal or corrected-to-normal hearing and vision, had no history of neurological disorders and had normal intelligence. The main inclusion criteria for the subjects were as follows: 1) no psychotropic medication treatment was taken at least two weeks before diagnosis and 2) no drug treatment was taken after diagnosis. To ensure the effectiveness and reliability of the EEG signals, the subjects were restricted to the age range of 18 to 55 years. Additionally, the subjects were required to have at least a primary school education level. To select subjects, psychiatrists used different scales, i.e., the Life Event Scale (LES), Mini-International Neuropsychiatric Interview (MINI), and Childhood Trauma Questionnaire (CTQ), to assess patient mental states based on the severity of depression, irritability, and stress levels. The subjects were informed of the purpose and protocol of the data acquisition We selected 15 depressed patients (female/male=9/6, 37.62±12.94 years old) and 20 healthy controls (female/ male=8/12, 31.33±7.48 years old) for the 64-channel EEG signal acquisition process. A 5-minute eyes-closed restingstate EEG signal was acquired for each subject using a 64-channel Brain Products (BP) electrode cap at a sampling rate of 1000 Hz with a 1-40 Hz analog bandpass filter. The impedance of each channel was maintained below 20 k . We set FCz and AFz as the reference and ground electrodes, respectively.
The 128-channel EEG signals in the MODMA dataset included 24 depressed patients (female/male=11/13, 30.88± 10.37 years old) and 29 healthy controls (female/male=9/20, 31.45±9.15 years old). The eyes-closed resting-state EEG signals of each subject were collected for 5 minutes using a 128-channel HydroCel Geodesic Sensor Net (Electrical Geodesics Inc.) at a sampling rate of 250 Hz.
The 3-channel EEG signals acquisition process included 81 depressed patients (female/male=50/31, 40.58±12.12 years old) and 89 healthy controls (female/male=48/41, 30.83± 10.32 years old). A 90-second eyes-closed resting-state EEG signal was acquired for each subject using a 3-channel wearable EEG collection device placed on the prefrontal lobe (Fp1, Fpz, and Fp2) at a sampling rate of 250 Hz [39]. We also used the device to acquire auditory stimulus-evoked EEG signals from 105 depressed patients (female/male=64/41, 41.02±12.00 years old) and 109 healthy controls (female/ male=60/49, 31.19±9.83 years old). In this work, we used six auditory stimuli with different emotional properties from the IADS-2 [4] to explore the differences between depressed patients and healthy controls, as illustrated in Table I. Each auditory stimulus was 6 seconds, and after each auditory stimulus, the patient was presented with 6 seconds of rest. Thus, 72 second EEG signals (6 auditory stimuli×6s+6 resting×6s=72 s) were acquired.
The raw EEG signals acquired with 64-channel and 128-channel electrode caps were preprocessed in MATLAB using the EEGLAB toolbox [40]. The sampling rate of the 64-channel raw EEG signal was downsampled to 250 Hz to reduce computational costs. The raw 128-channel and 3-channel EEG signals were filtered between 1 Hz and 40 Hz using a bandpass filter. Then, EEGLAB was applied to inspect all raw EEG signals for artifacts such as body movement, eye blinks, eye movement, and electromyogram (EMG).
All EEG signals were re-referenced using the common average reference (CAR) strategy. Then, the 64-channel and 128-channel EEG signals were decomposed to remove artifacts, such as electrooculogram (EOG), electrocardiogram (ECG), and EMG signals, using the independent component analysis (ICA) algorithm. Artifacts were removed from the 3-channel EEG signals with discrete wavelet transforms and Kalman filtering [41]. Due to noise, the EEG signals collected during and after the sixth auditory stimulus were removed. Finally, we used the processed 64-channel, 128-channel, and 3-channel resting-state EEG data of each subject and carefully selected four (4 × 10s = 40 s), eight (8 × 10s = 80 s), and four (4 × 10s = 40 s) 10-second valid epochs without artifacts, respectively, to build Dataset 1, Dataset 2 and Dataset 3. We divided the processed stimulus-evoked EEG data of each subject into ten 6-second (10 × 6s = 60 s) valid epochs to construct Dataset 4.

B. EEG Intrinsic Feature Extraction
EMD can be used to decompose complex signals into a finite number of IMFs, which reflect the intrinsic physical characteristics of the signal. The IMFs contain the time-frequency domain information of the signal and have the following properties: (1) The number of extremes (maximum and minimum values) is either equal to the number of zero crossings of the IMFs, or differs by one.
(2) The IMFs are locally symmetric, and the mean of the top and bottom envelopes of each IMF are zero.
(3) The IMFs contain at least one maximum value and one minimum value.
The decomposition process of the signal x(t) can be summarized as follows: (1) Find all extreme points in the signal x(t).
(2) Draw the envelopes e max (t) and e min (t) of the signal x(t) with cubic spline interpolation.
(3) Calculate the mean of the top and bottom envelopes m(t) = (e max (t )+e min (t )) 2 .
(4) Calculate the difference between the signal and the mean, i.e., (5) Determine whether there is a maximum value less than 0 or a minimum value greater than 0 in h(t). If so, repeat steps 1) to 4); if not, use c(t) = h(t) to represent an IMF. Then, repeat the above process with the residual signal until the signal x(t) cannot be decomposed further, and the remaining signal is a residue r (t).
After the decomposition process, the original signal can be represented as the sum of IMFs and the residue: where N is the number of IMFs decomposed by the EMD method.
To facilitate the subsequent operations in this paper, the time-series signals are represented as vectors with length t, i.e., the EEG signal is denoted by x, and the IMFs and residue are denoted by c 1 , c 2 , · · · , c N and r, respectively.
The relationship between IMFs and EEG signals illustrated in Eq. 1 can be used to determine the intrinsic features of the EEG signals. In the traditional intrinsic feature extraction method, two reference signals R 1 and R 2 are calculated by averaging the EEG signals of the depressed patients and the healthy controls in the training sets, respectively, as follows: where n 1 and n 2 are the numbers of EEG epochs with different labels in the training sets, and x 1i and x 2i are the i -th EEG epochs in the two groups. The reference signals are used to explore common information between different EEG epochs, which can be used to characterize the group properties of EEG signals.
Next, the reference signals are decomposed into IMFs: where N 1 and N 2 are the numbers of IMFs obtained by decomposing R 1 and R 2 using the EMD method, respectively, c 1 j and c 2 j are the IMFs, and r 1 and r 2 are the residues. Therefore, any EEG epoch x can be denoted by the IMFs and residues of R 1 and R 2 : where w 1 ∈ R 1×(N 1 +1 ) and w 2 ∈ R 1×(N 2 +1 ) are the expansion coefficients, and x 1 and x 2 are the EEG signals recovered using the IMFs of different reference signals. Let ) denote the matrices of the IMFs and residue acquired by decomposing the reference signals R 1 and R 2 , respectively. Therefore, Eqs. 4 and 5 can be rewritten in the form of matrix multiplication: Then the expansion coefficient can be calculated with the following equations: Therefore, all EEG epochs can use the above equations to derive the expansion coefficients w 1 and w 2 as intrinsic feature vectors. By extracting the features in this manner, we can explore the similarity between different EEG epochs and mental states to identify the intrinsic characteristics of the EEG signal. Theoretically, EEG signals that are characteristic of the depressed states should be more similar to the reference signals of depressed patients and less similar to the reference signals of healthy controls, and vice versa.

III. A REGULARIZATION PARAMETER-BASED IMPROVED INTRINSIC FEATURE EXTRACTION METHOD VIA EMD
According to N. E. Huang et al. [35], the orthogonality of the IMFs decomposed by the EMD method is local. For some data, the neighboring components of the IMFs may have the same frequency information, which leads to mode mixing. Therefore, the IMFs of the reference signals calculated by Eq. 3 could be linearly correlated, and Eqs. 8 and 9 may fail to effectively extract features. This result indicates that the traditional intrinsic feature extraction method has certain limitations and cannot be used with all types of data.
Since the IMFs could be linearly correlated, the determinant values of the matrix products A T 1 A 1 and A T 2 A 2 in Eqs. 8 and 9 are close to 0. In this case, the inverse of the matrix product does not exist. Therefore, the method of extracting the expansion coefficients as features is ineffective.
To address this problem, we propose an improved intrinsic feature extraction method that uses a regularization parameter to extract valid EEG features as accurately as possible. We let H 1 = A T 1 A 1 and H 2 = A T 2 A 2 denote the matrix products. Then, the rank of the matrix products H 1 and H 2 can be calculated as follows: If rank 1 is equal to N 1 +1, Eq. 8 can be utilized to effectively extract the intrinsic features of the EEG signal. However, if rank 1 is less than N 1 + 1, the following steps need to be taken to extract the features: Eq. 13 is equivalent to the following equation: Since H 1 is singular, the regularization parameter λ is introduced: where E denotes the identity matrix and λE denotes the regularization matrix. The regularization parameter λ is used to adjust the value of the identity matrix E to control its effect on the matrix product. The regularization parameter λ ensures that the matrix H 1 is as similar as possible to the matrix product H 1 . The value of the regularization parameter λ is set in the range of (0,1]. Since the matrix product H 1 can be decomposed into a nonnegative diagonal matrix and two unitary matrices by singular value decomposition (SVD) and λE is a nonsingular matrix, the matrix H 1 is a nonsingular matrix and its inverse exists. To verify that the matrix H 1 is nonsingular, we rewrite Eq. 15 as follows: Since 1 is a diagonal matrix with nonnegative elements and λE is a diagonal matrix with positive diagonal elements, the sum of 1 and λE is a nonsingular diagonal matrix with positive diagonal elements. Thus, the matrix H 1 is a nonsingular matrix with an existing inverse matrix, that can effectively and reliably be used in the following steps of the intrinsic feature extraction method. The inverse matrix can be expressed as follows: Thus, the intrinsic features w 1 can be calculated as follows: If rank 2 is equal to N 2 +1, w 2 can be calculated using Eq. 9. If rank 2 is less than N 2 + 1, we need to calculate w 2 as the intrinsic features of an EEG signal in the same manner as w 1 : It should be noted that although we demonstrated the regularization parameter-based improved intrinsic feature extraction method with only two groups of subjects, our method is also effective for multiple groups of subjects, which is discussed briefly below. First, the average EEG signals of different groups of subjects are used to acquire more reference signals with Eq. 2, such as R 3 , R 4 , · · · , R k . Then, the reference signals are decomposed into multiple IMFs using the EMD method, and the inverse matrices of the sum matrices of the reference signals are calculated using Eqs. 14 to 18. Finally, the intrinsic features of the EEG epochs are extracted by calculating the expansion coefficients from reference signals with different ground-truth labels with Eqs. 19 to 21.
The steps of the regularization parameter-based improved intrinsic feature extraction method are summarized in Algorithm 1. The improved method can effectively extract the intrinsic features of EEG signals by exploring the intrinsic differences between various groups. The method has a high generalizability and robustness, allowing it to extract intrinsic features from a wider range of time-series signals by neglecting their complex and nonstationary nature. For many types of time-series data, the improved intrinsic feature extraction method is as effective as traditional methods. However, in some situations, such as mode mixing, the improved method can effectively extract features that reflect the intrinsic characteristics of the signals, while traditional methods fail to extract the intrinsic features of time-series signals. Therefore, the improved method performs better than traditional methods. As illustrated in Algorithm 1, the computational complexity of the improved method is equivalent to that of traditional methods. Compared with traditional methods, the improved method requires slightly more time to calculate the ranks and the sum matrix. As a result, the improved method generally outperforms traditional methods. In addition, compared with other manual features, our method can more effectively explore the intrinsic differences between different groups of EEG signals and can infer differences in the neurophysiological mechanisms of various electrophysiological activities, resulting in an improved performance. Compared with the deep representations in end-to-end models, our method has more natural physics and mathematics interpretations, has a reduced computational cost and effectively prevents overfitting. The analysis of the proposed method shows that our method comprehensively considers intrinsic differences between groups, natural interpretations and the computational cost. Thus, our method can be used to effectively extract the intrinsic features of EEG signals.

Algorithm 1 A Regularization Parameter-Based Improved
Intrinsic Feature Extraction Method via EMD Input: n segments of training EEG epochs x; m segments of testing EEG epochs x; Training label y Output: The intrinsic features w Method: 1: Calculate the reference signals R 1 and R 2 in different class in training signals, respectively; 2: Compute the IMFs and residues of R 1 and R 2 and let A 1 and A 2 be the set of IMFs and residues of R 1 and R 2 , respectively; 3: while i ≤ n + m do 4: Extract the features of i-th subjects: w 1i · A T 1 = x i and w 2i · A T 2 = x i ; 5: Compute H 1 = A T 1 A 1 and H 2 = A T 2 A 2 , and compute rank rank 1 and rank 2 of H 1 and H 2 , respectively; 6: if rank 1 = N 1 + 1 then 7: if rank 2 = N 2 + 1 then 12: The extracted intrinsic features of the i -th EEG epochs are w i = [w 1i , w 2i ]; 17: end while 18: return the intrinsic features w The framework of the regularization parameter-based improved intrinsic feature extraction method is shown in Fig. 1, which we can use to explore essential properties of the improved method. The reference signals in Fig. 1 are acquired by averaging EEG signals from channels with the same label. We use different colors to represent the different channels of the reference signal. In the feature extraction process, the reference signals are first decomposed using the EMD method to acquire the IMFs of different channels. Then, a matrix is used to denote the IMFs and residuals, and the matrix product is calculated. To calculate the inverse of the matrix product A T A as accurately as possible, this work introduces a regularization parameter λ to the intrinsic feature extraction method. The identity matrix E is adjusted by the regularization parameter λ to ensure that the sum of the matrix product and the regularization matrix is both similar to their product and invertible. Since the sum matrix is nonsingular, its inverse can be calculated. Finally, for any EEG epoch, this inverse matrix can be applied to extract the intrinsic features w, improving the effectiveness and reliability of our method. It should be noted that the extracted intrinsic features of the EEG signals are similar to the reference signals of different classes, and we can concatenate these features to reveal the intrinsic characteristics of different groups in various EEG epochs.
In this section, the regularization parameter-based improved intrinsic feature extraction method was presented in several ways, including a formal definition, an algorithmic description, and a graphical demonstration. In the improved method, we first introduce a regularization parameter λ to adjust the value of the identity matrix E, which ensures that the sum matrix and the matrix product of the IMFs are similar. Then, we calculate the inverse of the sum matrix, which can be used to approximately represent the inverse of the matrix product, to effectively extract the intrinsic features in the following steps of the method, thereby optimizing the intrinsic features of EEG signals.

IV. EXPERIMENTS AND RESULTS
In this section, we conduct a series of experiments with four EEG datasets to further demonstrate the performance of the proposed intrinsic feature extraction method at exploring intrinsic characteristics of EEG signals and evaluate the performance of the method for depression recognition. For notational convenience, we abbreviate the regularization parameter-based improved intrinsic feature extraction method via EMD as RIEMD in this section. We first analyze the effect of the selected regularization parameter on the RIEMD method. Then, we compare the performance of the RIEMD method with various feature extraction methods in terms of depression recognition.
In this work, we applied support vector machine (SVM), decision tree (DT), and k-nearest neighbor (kNN) models to perform the classification experiments. We used the subjectindependent 10-fold cross-validation strategy to identify suitable parameters for the models. An effective strategy for identifying the parameters is to exponentially grow sequences of σ and C, such as C = [2 15 , 2 14 , · · · , 2 −15 ] and σ = [2 15 , 2 14 , · · · , 2 −15 ] [42]. To construct an optimal DT model, we tune the depth d and the minimum number of samples at a leaf node s. We also linearly grow k to optimize the kNN model, i. e., k = [1, 2, · · · , k max ].
We applied the min-max normalization strategy to normalize the training features to range between 0 and 1. The testing features were normalized with the same strategy using the maximum and minimum values acquired with the training features. In this work, we utilized the accuracy, sensitivity, and specificity to evaluate the classification performance of the RIEMD method. The classification experiments were performed on the platform with a 3.2 GHz CPU and 16 GB of RAM and the Windows 10 operating system. We trained the SVM model using LIBSVM [42] tools in the Java version and implemented the DT and kNN models on the Python platform.

A. Performance Comparison of the Regularization Parameter Selection
We investigated the effect of the regularization parameter on the RIEMD method and selected five parameters {0.1, 0.3, 0.5, 0.7, 0.9} to extract features from four EEG datasets. In this comparison, the features were classified using the SVM model for depression recognition. For Dataset 4, we demonstrated the average results of all EEG segments using the same parameters.
As shown in Table II, the depression recognition results varied little as the regularization parameter changed. The same results were acquired with Datasets 1 and 2 for different parameters λ. This result indicates that the RIEMD method is effective for extracting intrinsic features of EEG signals. The difference between the best and worst classification results on Datasets 3 and 4 was less than 0.0060, verifying the robustness and effectiveness of the RIEMD method. This result is due to the fact that the RIEMD method changed only the values on the principal diagonal through the introduction of the regularization parameter when calculating the inverse of the matrix product. Since the value of λ is small, the values on the principal diagonal change little, while the values on the nonprincipal diagonal do not change. The change of the regularization parameter had little effect on extracting the intrinsic features of EEG signals, which is the main reason for the robust performance of the RIEMD method in depression recognition, as illustrated in Table II.
Although the RIEMD results varied slightly with different regularization parameters, the best results on the four datasets were acquired with a regularization parameter of 0.1, as shown by the results in Table II. The smaller the regularization parameter is, the smaller the change in the matrix product. However, the regularization parameter cannot be extremely small, as this would lead to singularization of the matrix product, and the corresponding inverse matrix cannot be solved effectively. While the change in the regularization parameter has little effect on the extracted features, we recommend a regularization parameter of 0.1 to ensure the effectiveness and generalizability of the RIEMD method.

B. Performance Comparison of the Different Feature Extraction Methods
To investigate the effectiveness of the RIEMD method, we carefully chose three linear features in the EEG power spectrum, namely, the maximum frequency, mean frequency, and centroid frequency [43], [44], [45], and two nonlinear features, namely, the Kolmogorov entropy and LZ complexity [30], [46], [47], as our comparison baseline. In addition, we selected the EMD method proposed in [37], the differential entropy (DE) method proposed in [48], the CSP method used in [24] and the improved EMD (iEMD) feature extraction method proposed by Shen et al. [4] as the feature extraction methods for comparison. We also selected several state-of-the-art deep representations in end-to-end models as comparison methods, including the CNN model proposed by Zhang et al. [32], a deep neural network (DNN) model, a parallel model with a two-layer CNN and two-layer RNN (C_RNN) proposed by Saha et al. [49], and a covariance manifold-based long short-term memory (LSTM) model proposed by Zhang et al. [23]. These methods have been widely applied in EEG recognition tasks and have been demonstrated to be effective for depression recognition [46], [50]. The selected methods cover the main EEG feature extraction methods used in depression recognition research, as mentioned in Section I. In this paper, this comparison can be used to synthetically validate the superiority of the RIEMD method. The results of the experiments are demonstrated in Tables III and IV, where the results on Dataset 4 are the average of all EEG segments.
As shown in Table III, for EEG Datasets 1 and 2, the proposed RIEMD method achieved the best results, outperforming the other methods. Most of the depression recognition results obtained using the different feature extraction methods performed above the chance level of 0.5 on both datasets, except for the DT model for DE features on Dataset 1 and the kNN model for the baseline on Dataset 2, where the classification results were only 0.4188 and 0.4650, respectively. Therefore, these models are not effective and generalizable for depression recognition tasks. Furthermore, the depression recognition accuracy using both the DE and baseline methods was less than 0.6500, while the accuracies of the EMD, iEMD and our RIEMD methods was greater than 0.7700 on both datasets. This result indicates that these methods perform better than the DE and baseline methods in exploring the intrinsic characteristics of EEG signals. Moreover, this result demonstrates that using intrinsic features to represent EEG signals is more effective than using frequency domain features, temporal domain features and nonlinear features for depression recognition. Although the classification results achieved with the CSP method were better than those obtained with the DE and baseline methods, the results of CSP method acquired with the kNN model on Dataset 1 were considerably worse than those obtained using other methods. Furthermore, the sensitivity and specificity metrics achieved with the CSP method using the kNN and DT models differed significantly on Dataset 1 in depression recognition tasks, indicating a high risk of imbalance between the models. Therefore, the advantage of the CSP method over DE and baseline methods may not be stable enough.
On Datasets 3 and 4, the RIEMD, iEMD and CSP methods achieved one, three and three of the best metrics in the depression recognition results, respectively. Although the CSP method achieved the three best results, these results varied greatly from the remaining results achieved by the CSP method, while the best results achieved with the RIEMD and iEMD methods do not vary considerably from the remaining results. This result indicates that the RIEMD and iEMD methods are more stable than the CSP method. The depression recognition results obtained on both datasets were better than 0.5, and most of the results were greater than 0.6 (only the accuracy of the baseline with the DT model on Dataset 4 was 0.5945), which suggests that the effect of the different models on Datasets 3 and 4 reaches an acceptable level. In addition, Table III shows that the stability of the classification results on Datasets 3 and 4 was better than the stability on Datasets 1 and 2 because the results with different features and classifiers were similar. Among all the features, the DE and baseline methods have the greatest differences in overall stability, which may be caused by two reasons. First, these feature extraction methods do not optimize the spatial information of multichannel EEG data, and sufficient EEG spatial information cannot be learned by concatenating features extracted from different channels. The second reason is the lack of training EEG epochs in Datasets 1 and 2. Datasets 1 and 2 included only 53 and 35 subjects, respectively, while there were 170 and 214 subjects in Datasets 3 and 4, respectively. Additional subjects would be beneficial for the DE and baseline methods to investigate interclass similarities and differences. In addition, the RIEMD method demonstrated stable recognition results with high reliability and robustness on the four datasets, illustrating the effectiveness of optimizing intrinsic features of EEG signals for depression recognition.
As illustrated in Table IV, we investigated the performance of the proposed RIEMD method and the selected deep representations in end-to-end models. In the depression recognition tasks, the RIEMD method outperformed the CNN, DNN, C_RNN and LSTM models. The accuracies on the four EEG datasets were improved by approximately 0.0425 to 0.2000, 0.1100 to 0.2850, 0.1367 to 0.2132 and 0.0950 to 0.1404, respectively. As end-to-end models, the CNN and LSTM achieved better classification results than the DNN and C_RNN models. This result is mainly because the DNN and C_RNN models do not comprehensively consider the temporal and spatial domain characteristics of EEG signals and represent EEG data from only a single domain, while the CNN and LSTM models explore the characteristics of EEG signals in multiple domains. Therefore, classification results of the DNN and C_RNN models are slightly reduced. Compared with deep representations, the RIEMD method can effectively explore the intrinsic characteristics of EEG signals for depression recognition tasks. Moreover, the RIEMD method has more natural interpretations in physics and mathematics, requires less computational time and prevent overfitting.
To further demonstrate the effectiveness of the RIEMD method, the Friedman test [51] and Bonferroni-Dunn test [52] were used to compare the performance of different feature extraction methods, as shown in Table III. The CNN, DNN, C_RNN and LSTM models were excluded from the test due to obtaining only one result for each dataset. Since there were six methods and twelve classification results for each method in the test, F 5,55 = 40.42 was calculated based on the Friedman test, which is much higher than the critical value of 2.38 for F(5,55) with α = 0.05. This test revealed significant ( p < 0.05) differences between the classification results. The post-hoc two-tailed Bonferroni-Dunn test was used to compare the various methods. In depression recognition tasks, the RIEMD method significantly ( p < 0.05) outperformed the DE and baseline methods, the iEMD method performed significantly ( p < 0.05) better than the CSP, DE and baseline methods, and the EMD method performed significantly ( p < 0.05) better than the DE and baseline methods. Thus, the proposed RIEMD method is highly effective for depression recognition and could optimize the intrinsic features of EEG signals. Although the classification results of our RIEMD method are not significantly better than those of the iEMD method, the results are comparable. Furthermore, the RIEMD method requires less computational time than the iEMD method. To guarantee the effectiveness of solving the inverse matrix of the matrix product of the IMFs, the RIEMD method must perform n more additions than the EMD method, while the iEMD must implement SVD with a computational complexity of O(n 3 ), where n denotes the size of the matrix product.
To further demonstrate the computational differences between the RIEMD and iEMD methods, we conducted two experiments. The first experiment involved the artificial matrix products of the IMFs, which were randomly generated for different sample sizes. In this work, we increased the number of input samples from 1 to 10000 and compared the computational cost. The comparison results are illustrated in Fig. 2.
According to the experimental results, the computational costs of the RIEMD and iEMD methods are consistent with the above analysis. As the sample sizes increased, the computational costs of the RIEMD and the iEMD methods increased. Moreover, the computational cost of the iEMD method increased more quickly than that of our RIEMD method.
The second experiment involved the four EEG datasets. The computational costs of the RIEMD and iEMD methods with these EEG datasets are shown in Table V, where RIEMD_INV and iEMD_INV indicate that we solved the inverse of the matrix product using the RIEMD and iEMD methods, respectively. The computational time of the RIEMD and iEMD methods was the average time to extract the features of one EEG epoch with the RIEMD or iEMD method, respectively. We repeated the computation of each EEG epoch 100 times to guarantee stable results. As shown in Table V, the advantage of the RIEMD method in terms of the computational time is consistent with the theoretical analysis and the results illustrated in Fig. 2. However, this advantage is not clear during the feature extraction stage because decomposing EEG signals using the EMD method requires more computational time than solving the inverse of the matrix product. Therefore, in consideration of the classification results and the computational time, the performance of our RIEMD method is better than that of the iEMD method. In addition, the RIEMD method provides a novel idea and approach for effectively solving the issue of inaccurate and impractical extraction of intrinsic features with the EMD method due to mode mixing and may inspire more innovative and intelligent applications in the field of depression recognition.

V. CONCLUSION
This study presented a regularization parameter-based improved intrinsic feature extraction method via EMD and its application in depression recognition. In this method, we introduced a regularization parameter to effectively optimize the intrinsic features of EEG signals with more natural physics and mathematics interpretations. The experimental results on four EEG datasets demonstrate that the RIEMD method can efficiently capture relationships between the intrinsic characteristics of EEG signals and depression, achieving good depression recognition performance. One of the limitations of our work is that the dataset is small. Besides, the channel space of EEG signals should contain more useful information for depression recognition. In future studies, we will collect more EEG signals to improve the generalizability of our method and further optimize the spatial information from EEG signals.