A flexible and accurate method for electroencephalography rhythms extraction based on circulant singular spectrum analysis

Markov process amplitude EEG model

Simulated spontaneous EEG signals were used to verify the validity of the CiSSA method in alpha rhythm extraction. Spontaneous EEG signals were generated based on the Markov Process Amplitude (MPA) EEG model (Bai et al., 2001; Nishida, Nakamura & Shibasaki, 1986). The MPA EEG model is a powerful and widely used method to simulate and interpret EEG signals. With a few parameters, the model can represent the two major characteristics of EEG signals: rhythmic oscillation and randomness. Rhythmic oscillation is represented by sinusoidal waves, and randomness was represented by the stochastic process amplitude of the first-order Markov process. In recent years, the MPA EEG model has been applied in several studies analyzing spontaneous EEG signals, which have employed such techniques as feature expression, quantitative analysis (Nakamura et al., 1997), and algorithm verification (Xu et al., 2018).

In the MPA EEG model, EEG signals consist of several rhythmic oscillations expressed by a sinusoidal wave

(10) $$s\left(n\mathrm{\Delta }t\right)=\sum _{i=1}^K{a}_i\left(n\mathrm{\Delta }t\right)\sin \left(2\pi f_{i}n\mathrm{\Delta }t+{\theta }_i\right)$$

where n is the number of samples, Δt is the time interval, K is the number of rhythms, f is the dominant frequency of rhythm, θ is the initial phase (zero), and a is the rhythmic amplitude obtained from the following first-order Gauss–Markov process:

(11) $$a_i\left[\left(n+1\right)\mathrm{\Delta }t\right]={\gamma }_i{a}_i\left(n\mathrm{\Delta }t\right)+{\xi }_i\left(n\mathrm{\Delta }t\right)$$

where γ is the coefficient of the first-order Markov process, and ξ is a random increment of Gaussian distribution with mean zero and variance ${\sigma }^{\xi }$. Therefore, the rhythmic amplitude at the succeeding time $\left(n+1\right)\mathrm{\Delta }t$ depends only on the amplitude at time $\mathrm{\Delta }t$ and is determined only by two parameters: γ and ${\sigma }^{\xi }$. The parameters of the MPA EEG model are determined in the frequency domain to achieve the maximum likelihood with respect to the power spectrum of real EEG. H_i is defined as the amplitude, and B_i is the frequency width at half of H_i of the EEG power spectrum. Based on the literature (Bai et al., 2001), H_i, B_i can be described as

(12) $$\{\begin{array}{c}{H}_i=\frac{\mathrm{\Delta }t{\left({\sigma }_i^{\xi }\right)}^2}{4{\left(1-{\gamma }_i\right)}^2}\hfill \\ B_i=\frac{1}{\pi \mathrm{\Delta }t}{\cos }^{-1}\frac{4{\gamma }_i-1-{\left({\gamma }_i\right)}^2}{2{\gamma }_i}\hfill \end{array}$$

The simulation procedures of spontaneous EEG signals based on the MPA EEG model are shown in Fig. 1. First, the power spectrum of a real EEG signal with sampling rate f_s = 200 Hz is calculated in Fig. 1A. Then, f_i, H_i, and B_i, which represent the peak frequencies, amplitude, and the frequency width at half of amplitude of the EEG rhythms (delta, theta, alpha, and beta), respectively, are obtained according to the power spectrum. Based on Eq. (12), the parameters of the first-order Gauss–Markov process (γ and ${\sigma }^{\xi }$) are obtained. All parameters of the MPA EEG model are shown in Table 1. Then, the delta, theta, alpha, and beta rhythms are simulated based on the determined parameters, as shown in Fig. 1B. Finally, the simulated EEG signal (shown in Fig. 1C) is generated as the sum of the four rhythms. The simulated spontaneous EEG lasts for 8 s with a 5-ms Δt interval.

Circulant singular spectrum analysis of simulated EEG signals

The simulated EEG signal is processed by the CiSSA method. The embedding dimension L should satisfy the condition $L\ge f_s/f_{band}=40$, with the data sampling rate $f_s=200\mathrm{H}\mathrm{z}$ and the bandwidth of alpha rhythm (8–13 Hz) $f_{band}=5\mathrm{H}\mathrm{z}$. Moreover, the calculated frequency bandwidth $f_s/L$ is always less than the real frequency bandwidth of each reconstructed component. It is better to select the embedding length L twice times as long as $f_s/f_{band}$. Therefore, the embedding length L was set to 40 or 80.

Figure 2 shows the power spectrum density (PSD) of the first six reconstructed components (RCs) for a simulated EEG signal. When L = 40, as shown in Fig. 2A, every RC falls on the theoretical frequency derived by Eq. (5). Furthermore, the bandwidth of each RC is limited to $f_s/L$ = 5 Hz. Similarly, when L = 80, as shown in Fig. 2B, all RCs fall on the theoretical frequencies with the bandwidth limited to 2.5 Hz. However, the simulated EEG signal is processed by the basic SSA method with the embedding dimension set to L = 40. The PSD of the first six RCs is shown in Fig. 3. The frequency of each RC is unknown. To group the RCs by frequency, other algorithms like Fourier transform are introduced to calculate the frequency of RCs. Besides, according to the PSD of RC5 and RC6, some components fall outside of the bandwidth limit of $f_s/L$ = 5 Hz. This phenomenon is called component mixing.

Because the frequencies of the RCs processed by the CiSSA method are known, the alpha rhythm of the EEG signal can be extracted by Eq. (8). The error parameter for evaluating the performance of the alpha rhythm extraction is defined as Xu et al. (2018)

(13) $${\epsilon }_{ave}=\frac{1}{N}\sum _{i=1}^N|P_{\alpha }\left(i\right)-P_e\left(i\right)|$$

where ${\epsilon }_{ave}$ is the average error of the PSD between the simulated alpha rhythm and the extracted alpha rhythm, $P_{\alpha }\left(i\right)$ is the PSD of the simulated alpha rhythm, $P_e\left(i\right)$ is the PSD of the extracted alpha rhythm, and N is the length of the PSD.

Figure 4A shows the extracted alpha rhythm of the simulated EEG signal by the CiSSA method with the embedding dimension set to L = 40. Figure 2A shows that RC3 represents the alpha rhythm. The PSD of the simulated and extracted alpha rhythm by the CiSSA method when L = 40 is shown in Fig. 4B. There was component mixing (slash shadow), and the error of the extracted alpha rhythm was $\epsilon =0.43\mu V^2/\mathrm{H}\mathrm{z}$. Similarly, RC5 and RC6 represent the alpha rhythm from Fig. 2B with an embedding dimension of L = 80. The alpha rhythm extracted from the combination of RC5 and RC6 is shown in Fig. 4C. The PSD of the simulated and extracted alpha rhythm by the CSSA method when L = 80 is shown in Fig. 4D. There was less component mixing (slashed shadow) than that observed when L = 40, and the error of the extracted alpha rhythm was $\epsilon =0.27\mu V^2/\mathrm{H}\mathrm{z}$. In the basic SSA method, the alpha rhythm was extracted according to the adaptive grouping rule (Hai et al., 2017). The extracted alpha rhythm by the basic SSA method when L = 40 is the sum of RC3, RC4, RC5, and RC6, as shown in Fig. 4E. The PSD of the simulated and extracted alpha rhythms by basic SSA when L = 40, shown in Fig. 4F, illustrates the presence of more component mixing (slashed shadow) than that found with the CiSSA method when L = 40, and the error of the extracted alpha rhythm was $\epsilon =0.89\mu V^2/\mathrm{H}\mathrm{z}$. Similarly, the extracted alpha rhythm by the basic SSA method when L = 80 is the sum of RC4, RC5, RC6, and RC11, as shown in Fig. 4G. Figure 4H shows the PSD of the simulated and extracted alpha rhythms by basic SSA when L = 80. The error of the extracted alpha rhythm by the basic SSA method (L = 80) was $\epsilon =0.34\mu V^2/\mathrm{H}\mathrm{z}$, larger than that by the CiSSA method (L = 80). Besides, the extracted alpha rhythm by the basic SSA method (L = 40 or L = 80) contained the frequency components outside the alpha rhythm bands (8–13 Hz).

To compare the performance of the alpha rhythm extraction with that of other methods, the alpha rhythms were extracted by the finite impulse response (FIR) filtering methods and the wavelet decomposition (WDec) method. The FIR filter is a simple, stable and narrowband filter that has been utilized in many digital applications like image processing, wireless communication, biomedical etc. (Khurshid & Mir, 2017; Mahabub, 2019). The WDec can be performed by a series of high and low-pass filters (Mamun, Al-Kadi & Marufuzzaman, 2013). As a time-frequency representation of a signal, WDec has some advantages over other conventional techniques as optimal resolution in both the time and frequency domains and lack of the requirement of stationarity of the signal (Akar et al., 2015; Quiroga et al., 2001). The PSDs of the extracted alpha rhythms from a simulated EEG signal by the FIR with order 60 and WDec methods are shown in Figs. 5A and 5B, respectively. The PSDs of the extracted alpha rhythm by the FIR and WDec methods had a higher magnitude, and there were more component mixing than that by the CiSSA method. The errors of the extracted alpha rhythm by the FIR and WDec method were $\epsilon =0.51\mu V^2/\mathrm{H}\mathrm{z}$ and $\epsilon =0.57\mu V^2/\mathrm{H}\mathrm{z}$, respectively, which were higher than that obtained by the CiSSA method. This is because other rhythms and artifacts (electrooculogram (EOG), electromyography (EMG), baseline drift and stochastic noise) also have components in the alpha band (8–13 Hz). The FIR filter and the WDec methods are not able to separate the component mixing overlapping in frequency space. However, the CiSSA method can suppress a part of noise with overlapping frequencies.

In order to increase the robustness and validity of alpha rhythm extraction, 1,000 simulated EEG signals have been performed and the average error of the PSD between the simulated alpha rhythm and the extracted alpha rhythm was calculated as

(14) $$\overline{\epsilon }={\epsilon }_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}\pm {\epsilon }_{\mathrm{s}\mathrm{t}\mathrm{d}}$$

where ${\epsilon }_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}$ is mean of errors ( ${\epsilon }_{ave}$) in 1,000 simulations, ${\epsilon }_{\mathrm{s}\mathrm{t}\mathrm{d}}$ is the standard deviation (STD) of errors ( ${\epsilon }_{ave}$) in 1,000 simulations. The average errors of the extracted alpha rhythm by the CiSSA, basic SSA, FIR and WDec methods are shown in Table 2. The mean and STD of the errors of the PSD between the simulated alpha rhythm and the extracted alpha rhythm by the CiSSA method (L = 80) were $\epsilon =0.29\mu V^2/\mathrm{H}\mathrm{z}$ and $\epsilon =0.08\mu V^2/\mathrm{H}\mathrm{z}$, respectively, which were much lower than those obtained by the basic SSA, FIR and WDec methods.

We conclude that the RCs of simulated EEG signals processed by the CiSSA method fall on the theoretical frequencies limited to the selected bandwidth ranges. The alpha rhythm can be extracted automatically based on the frequency feature. The alpha rhythm extraction by the CiSSA method performed better than that of the basic SSA, FIR, and WDec methods. Therefore, the processing results of the simulated EEG verify the validity of the CiSSA method’s performance in alpha rhythm extraction. Furthermore, the error of the extracted alpha rhythm by the CiSSA method varied with the embedding dimension. The calculated average errors of the extracted alpha rhythm from 1,000 simulated EEG signals by the CiSSA method with different embedding dimensions are shown in Table 3. The average error attained a minimum value at L = 80. Therefore, the embedding dimension of the CiSSA method for alpha extraction was set to L = 80.

Results and discussion of database EEG signals

The database EEG signals reported in literature (Trujillo, Stanfield & Vela, 2017) were used. A total of 22 subjects (11 male and 11 female) aged from 18–26 participated in the experiment. The experiments were carried out with the subjects sitting quietly on a comfortable chair in a darkened room. Every subject underwent 8 min of eyes-open and eyes-closed resting state, with 4 min eyes open and 4 min eyes closed interleaved in 1-min intervals. The EEG signals of one subject (subject #6) were removed because of a technical recording error. A total of 72 channels of EEG signals were recorded by a BioSemi Active II amplifier system. The EEG signals were downsampled to 256 Hz.

Channel Fpz was selected for EEG analysis. The EEG data of the Fpz channel were divided into 8-s (2,048 samples) epochs with 50% overlap, thus producing 91 epochs of eyes-open and eyes-closed conditions for each subject. This was done because artifacts, including those resulting from electrooculogram, electromyography, baseline drift, and stochastic noise, interfere with the rhythm extraction. The adaptive SSA method (Hai et al., 2017) was used to remove the artifacts, and the results are shown in Fig. 6. Figures 6A and 6B show EEG epochs of the eyes-open and eyes-closed conditions of subject #17, respectively. Figures 6C and 6D show the corrected EEG signals after artifact removal. The electrooculogram artifacts were removed from the EEG signals of the eyes-open condition, and the spontaneous EEG signals were preserved in both conditions.

The EEG signals after artifact removal were processed by the CiSSA method with the embedding dimension set to L = 80. The first six RCs of the EEG signals in the eyes-open and eyes-closed conditions are shown in Figs. 7A and 7B, respectively. Each RC falls on the theoretical frequency derived by Eq. (5), and the bandwidth of the RCs is limited to $f_s/L$ = 3.2 Hz, which is agreement with the simulation results. Figures 7A and 7B show that RC4 and RC5 represent the alpha rhythm in both the eyes-open and eyes-closed conditions, respectively. Thus, the alpha rhythms of the EEG signals can be extracted automatically as the sum of RC4 and RC5, which is agreement with Eq. (8). Figures 8A and 8B show the extracted alpha rhythms of the EEG signals in the eyes-open and eyes-closed conditions, respectively. The amplitude of the alpha rhythm in the eyes-open condition was lower than that in the eyes-closed condition. This was consistent with the results of previous studies, in which the alpha rhythm in the resting state in the eyes-open condition with visual stimulation was much weaker than that in the eyes-closed condition (Barry et al., 2007). Figure 8C illustrates the spectrogram of alpha rhythms in the eyes-open and eyes-closed conditions, which is the square of the rhythm’s amplitude as a function of time and frequency. It illustrates a significant difference between the eyes-open and eyes-closed states.

The performance of alpha rhythm extraction by the CiSSA method was compared with that of three other methods: the basic SSA method, the WDec method, and the FIR method. The alpha rhythms under the eyes-open and eyes-closed conditions were extracted using the CiSSA, basic SSA, WDec, and FIR methods. The PSD of the extracted alpha rhythms by the four methods under the eyes-closed and eyes-open conditions are shown in Fig. 9. Figures 9A and 9D show that the extracted alpha rhythms using the CiSSA method were within the alpha band (8–13 Hz) under both the eyes-open and eyes-closed conditions. In addition, the power of the extracted alpha rhythm under the eyes-open condition was lower than that under the eyes-closed condition. Therefore, the alpha rhythm extracted using the CiSSA method could represent the real EEG alpha rhythm. However, the alpha rhythm extracted by the basic SSA method under both the eyes-open and eyes-closed conditions contained frequency components less than 8 Hz, which were outside the alpha band. Similarly, the alpha rhythms extracted using the WDec method under eyes-closed condition contained components outside the alpha band (component mixing), as shown in Fig. 9E. Figures 9B, 9C and 9F show the alpha rhythms extracted using the WDec method under eyes-open condition, using the FIR method under eyes-open and eyes-closed conditions, respectively. The extracted alpha rhythms fell into the alpha band, and that the power of the extracted alpha rhythm was stronger than that extracted using the CiSSA method. This is in agreement with the simulation results shown in Fig. 5 because the WDec and FIR methods are unable to remove artifacts and noise from the alpha rhythm with an overlapping frequency spectrum. Therefore, the CiSSA method performed better than the basic SSA, WDec, and FIR methods at alpha rhythm extraction.

To further verify the CiSSA method’s performance, the extracted alpha rhythms were used to distinguish between the eyes-open and eyes-closed states, and the classification results produced by the CiSSA method were compared with those by the basic SSA, FIR, and WDec methods. In this study, the power $\left(P={\sum }_{i=1}^N{V}_i^2/N\right)$ and the mean of the absolute value $\left(\overline{V}={\sum }_{i=1}^N|V_i|/N\right)$ were selected as the features of the alpha rhythm (Mohammadi et al., 2015), where $V_i$ represents the amplitude of the extracted alpha rhythm, and N represents the number of samples. Figure 10A shows the values of the power and the mean of the absolute value of the extracted alpha rhythm by the CiSSA method for subject #17. The power value and the mean of the absolute value under the eyes-open condition were lower than those under the eyes-closed condition. Then, the support vector machine method was used to classify the features under the eyes-open and eyes-closed conditions. The classification accuracy was 92.31%. Figures 10B–10D show the power values and mean absolute values of the extracted alpha rhythms by the basic SSA, FIR, and WDec methods, respectively. Similar to the results obtained by the CiSSA method, the power and mean of the absolute value of the extracted alpha signal under the eyes-open condition were generally lower than those under the eyes-closed condition. The classification accuracy of feature extraction by the basic SSA, FIR, and WDec methods was 89.01%, 91.21%, and 90.11%, respectively, and these values were lower than the accuracy obtained by the CiSSA method.

We calculated the power values and means of the absolute value of the extracted alpha rhythms of all 21 subjects, and the classification results are shown in Table 4. The classification accuracy varied greatly between different subjects because of individual differences in EEG signals. The mean and standard deviation of the classification accuracy was calculated for all subjects to compare classification performance between the CiSSA, basic SSA, FIR, and WDec methods. The mean value of the classification accuracy for all subjects by the CiSSA method was 92.36%, which was higher than those obtained by the basic SSA (86.87%), FIR (91.58%), and WDec (90.89%) methods. The standard deviation of the classification accuracy across all subjects by the CiSSA method was 7.05%, which was lower than those obtained by the basic SSA (11.78%), FIR (7.42%), and WDec (8.39%) methods. Therefore, the CiSSA method’s classification performance was better and more robust than that by the basic SSA, FIR, and WDec methods.

Results and discussion of experimental EEG signals

Additional experimental EEG signals were recorded and used to further verify the validity of the CiSSA method. The Institutional Review Board of Tsinghua University granted ethical approval (20170010) to carry out the experiment within its facilities. Ten subjects (8 male and 2 female) aged 20–29 participated in the experiments and the written informed consent was obtained from the subjects. The experiments were carried out with the subjects sitting on a comfortable chair in a room with normal lightness. The experimental EEG signals were recorded using the MP160 data acquisition and analysis system (BIOPAC Systems, Inc., Goleta, CA, USA) with the sampling rate of 200 Hz. The active Ag/AgCl electrode flushed with conductive gel was used as the recording electrode and attached to the frontal region of the subject’s scalp. Another two electrodes, which served as ground and reference, were attached to the earlobe and mastoid, respectively, as shown in Fig. 11B. Every subject underwent 30 s of eyes-open resting state, followed by 30 s of eyes-closed resting state, and repeated this procedure 57 times. Segments lasting 8 s were extracted from the middle of each period, as shown in Fig. 11A. Consequently, 57 segments each of the eyes-open and eyes-closed states for each subject were obtained. Artifacts were removed from each EEG segment using the adaptive SSA method (Hai et al., 2017).

The experimental EEG signals were processed by the CiSSA method, and the alpha rhythms recorded under the eyes-open and eyes-closed conditions were extracted. The power and mean absolute values of the extracted alpha rhythms were calculated as features for classification using the support vector machine method. The classification accuracy by the CiSSA method for subject 1# was 91.23%, which was higher than that obtained by the basic SSA (87.72%), FIR (89.47%), and WDec (88.60%) methods, as shown in Fig 12. The classification results of all 10 subjects are shown in Table 5. The mean value of the classification accuracy for all subjects by the CiSSA method was 92.11%, which was higher than those obtained by the basic SSA (87.46%), FIR (90.17%), and WDec (88.77%) methods. The standard deviation of the classification accuracy across all subjects by the CiSSA method was 4.74%, which was lower than those obtained by the basic SSA (5.86%), FIR (5.79%), and WDec (9.69%) methods. We therefore concluded that the CiSSA method’s classification performance was better than that by the basic SSA, FIR, and WDec methods.