An adaptive singular spectrum analysis method for extracting brain rhythms of electroencephalography

The adaptive grouping rule

With the SSA treatment the original EEG time series are decomposed into a set of RCs. The first several RCs dominate the trend of the EEG time series, which is represented by the large artifacts (Teixeira et al., 2006). Here, the first one or two RCs, determined by the EEG amplitude, are grouped as artifacts. When the amplitude is large, indicating a high level of artifacts, the first two RCs are grouped. Otherwise, the first RC is grouped: (3)${\displaystyle Artifacts=\left\{\begin{array}{c}\hfill R{C}_1+R{C}_{2}max\left(\mathbf{s}\right)>V_0\hfill \\ \hfill R{C}_{1}max\left(\mathbf{s}\right)<V_0\hfill \end{array}\right.}$ where max(s) denotes the EEG amplitude V₀ is a threshold. Since the amplitude of the spontaneous EEG without artifacts is usually below 100  µV (Ng & Chan, 2005), the threshold was set as V₀ = 200 µV.

After artifacts removal, the rhythms are then extracted from the EEG time series. Rhythms are the oscillatory components of EEG time series, including, but not limited to, the periodic components. In order to extract the EEG rhythms, the RCs are then divided into two groups: the periodic and non-periodic components. Firstly, the periodic components are extracted. Generally, a periodic time series will be factorised into some eigenvalue pairs with similar amplitude using SSA (Vautard, Yiou & Ghil, 1992). So, RCs with similar eigenvalues belong to the same periodic component. They are summed up as a periodic component (PC). The similarity of the eigenvalues for periodic component extraction is determined by the following criterion: (4)${\displaystyle |1-\frac{{\lambda }_j}{{\lambda }_i}|<0.05.}$ Then, these PCs and the remaining non-periodic RCs are grouped. Each PC and RC will fall into a narrow frequency band, when the window length L is large enough (Kouchaki, 2014; Mohammadi et al., 2015; Sanei, Ghodsi & Hassani, 2011). So, the peak frequency in the Fourier transform can be used to represent the frequency range of the PC and RC. (5)${\displaystyle f_{\max }=\underset{f}{argmax}\left\{abs\left[FFT\left(RC\right)\right]\right\}}$ where f_max is the peak frequency, FFT(RC) is the fast Fourier transform of the RC or PC. RCs and PCs with the peak frequencies in the same rhythm band are finally clustered into the same group, which constitutes the brain rhythm.

The pseudo-code of the adaptive SSA method is shown in Fig. 1.

Window length L is selected based on the lowest frequency of interest (L ≥ f_s∕f_l) to capture at least one period of the expected component (James & Lowe, 2003). In this paper, the window length L is set to be 40. After extracting the desired rhythms, the features of the rhythms are obtained for each EEG signal. The mean, standard deviation (SD), power and power ratio of different rhythms are usually selected as the features (Mohammadi et al., 2015). In this study, the power ($P=\sum V_{\alpha }^2∕N$) of the alpha rhythm were selected.

Simulated EEG data

Simulated EEG data consisted of two parts: the spontaneous EEG and the artifacts. The spontaneous EEG has two major characteristics: rhythmic oscillations and randomness. The artifacts mainly consist of EOG, baseline drift and white Gaussian noise.

The spontaneous EEG was simulated based on the Markov Process Amplitude (MPA) EEG model developed by Nishida, Nakamura & Shibasaki (1986). In the MPA EEG model, the EEG rhythmic oscillations were represented by sinusoidal waves, and the EEG randomness was represented by the stochastic process amplitude of the first-order Markov process (Al-Nashash et al., 2004; Bai et al., 2000; Bai et al., 2001). The spontaneous EEG was generated by a combination of K different oscillations, (6)${\displaystyle s\left(n\Delta t\right)=\sum _{i=1}^K{a}_i\left(n\Delta t\right)sin\left(2\pi m_{i}n\Delta t+{\theta }_i\right)}$ where n is the sample number, Δt is the time interval, K is the number of rhythms, m is the dominant frequency, θ is an arbitrary value representing the initial phase, a is the rhythmic amplitude obtained from a first-order Gauss–Markov process (7)${\displaystyle a_i\left[\left(n+1\right)\Delta t\right]={\gamma }_i{a}_i\left(n\Delta t\right)+{\xi }_i\left(n\Delta t\right)}$ where γ is the coefficient of the first-order Markov process, ξ is a random increment of Gaussian distribution with zero mean and variance σ.

Figure 2 shows the procedures of the spontaneous EEG simulation based on the MPA EEG model. Firstly, the power spectrum of a real EEG was calculated, as shown in Fig. 2A. In order to achieve the maximum likelihood with respect to the power spectrum of the real EEG, parameters of the MPA EEG model were determined in the frequency domain (shown in Table 1). Then, four oscillations representing the brain rhythms (delta (1–4 Hz), theta (4–8 Hz), alpha (8–13 Hz) and beta (13–30 Hz)) were generated based on the determined parameters, as shown in Fig. 2B. Finally, the spontaneous EEG was generated by a combination of the four rhythms as shown in Fig. 2C. The simulated spontaneous EEG lasted for 8 s with the interval Δt of 5 ms. The detailed mathematical description of the spontaneous EEG simulation is shown in the ‘Appendix’ (Markov Process Amplitude (MPA) EEG model).

The artifacts consist of EOG, baseline drift and white Gaussian noise. EOG is the main artifact and is caused by eye blinks and ocular movement. It is characterized by large amplitude, low-frequency electro-potential shift. Baseline drift originates from the head or body movement and is characterized by low-frequency electro-potential shift. White Gaussian noise was used to represent the measurement noise.

Generally, about 15–20 eye blinks, each lasting 0.3–0.4 s will be done in 60 s when the subject is relaxed. Based on the characters, a triangular pulse with the period of 3 s and the pulse width of 0.3 s was chosen to simulate the EOG artifact. A sinusoidal function with amplitude of 20 µV and frequency of 0.5 Hz was used to simulate baseline drift.

Parameters used in the EEG simulation are shown in Table 1. The data simulation and treatment were carried out using MATLAB version 7.9.0 on a personal computer with Intel(R) Core(TM) i5-4590 processor, 8 GB RAM and Windows 7 operating system.

Experimental EEG data

The experiments were approved with a protocol (NO. 20170010) by the Institutional Review Board of Tsinghua University. One male and two female aged 20–25 years participated in the experiments. These subjects were required to abstain from psychoactive substances for at least 4 h prior to experiments. Experiments were carried out with the subject sit on a comfortable chair in a room with normal lightness. The MP36 data acquisition and analysis system (BIOPAC Systems, Inc., Goleta, CA, USA) with a three electrodes system was used to acquire the EEG data. The Ag/AgCl electrode (Wuxi Sichiray Technology Co. Ltd, Wuxi, Jiangsu, China) flushed with conductive gel was attached to the scalp over the frontal region as the recording electrode. The other two electrodes were attached to the earlobe and mastoid, serving as a ground and a reference, respectively.

The experimental procedures were as follows. Initially, the subject closed eyes in relaxed state for 10 min. After that, the subject opened eyes and visually fixated on a small cross displayed on a computer screen in front of him/her. Meanwhile, the EEG data started to be recorded with the sampling rate of 200 Hz. Informed by the recorder, the subject then began 4 times of alternating periods of 60 s eyes open followed by 60 s eyes closed. The desired EEG segments were cut off from every period of the eyes-open and eyes-closed states. As shown in Fig. 3, each segment last for 8 s and was in the middle of each period. For each subject, the experiment was repeated three times in three separate days. Totally, 24 segments were obtained for each subject.

Simulated EEG data

Simulated EEG data based on the MPA EEG model was used to verify the artifacts removal using the adaptive SSA. Representative results are shown in Fig. 4. When the EEG was contaminated with large artifacts (the EOG amplitude was 400 µV ), as shown in Fig. 4A, the amplitude of the EEG signal was higher than the threshold (V₀ = 200 µV). The first two RCs were grouped as artifacts and removed. The corrected EEG and removed artifacts are shown in Figs. 4B and 4C, respectively. It can be seen that the amplitude of the corrected EEG was 54  µV, close to that of the simulated spontaneous EEG as shown in Fig. 2C. When the EEG was contaminated with small artifacts (the EOG amplitude was 100 µV), as shown in Fig. 4D, the first RC was grouped as artifacts and removed. The amplitude of the corrected EEG was 49 µV, close to that of the spontaneous EEG as well. Therefore, it could be concluded that the large or small artifacts could be removed adaptively by the adaptive SSA.

Experimental EEG data

The artifacts removal using the adaptive SSA was then further tested on experimental EEG data. Another recently reported SSA method (Maddirala & Shaik, 2016b), which will be called SSA 1# in this paper, was used for comparison. The SSA 1# used a novel grouping criterion based on the eigenvectors’ local mobility, which is a signal complexity measure, to remove motion artifacts. The comparison results are shown in Fig. 5. For nine of the total 24 segments, the corrected EEG using the adaptive SSA was the same as that using the SSA 1#. Representative results are shown in Figs. 5B and 5C. Furthermore, it can be seen from Fig. 5D that their power spectrums overlapped completely. For the other 15 of the total 24 segments, the amplitude of the corrected EEG using the adaptive SSA was higher than that of the corrected EEG using the SSA 1#, as shown in Figs. 5F and 5G, respectively. It can be seen from Fig. 5H that the SSA 1# removed more artifacts than the adaptive SSA. The excess removed artifacts were in the frequency of 3–5 Hz, which was within the EEG frequency band and should not be removed. Therefore, the corrected EEG using the adaptive SSA was more complete and accurate than that using the SSA 1#. It verifies that the adaptive SSA has a similar, even better, artifacts removal effect.

Simulated EEG data

Simulated EEG data based on the MPA model was then used to verify the rhythms extraction using the adaptive SSA. Results are shown in Fig. 6. The extracted alpha rhythm using the adaptive SSA is shown in Fig. 6A. It was similar with the simulated alpha rhythm (shown in Fig. 6B). To further examine the rhythms extraction, their power spectrums were calculated (shown in Fig. 6C). The power spectrum of the extracted alpha rhythm was in good consistence with that of the simulated alpha rhythm. It verifies the validity of the adaptive SSA in rhythms extraction.

Experimental EEG data

The rhythms extraction using the adaptive SSA was then further tested on experimental EEG data. After artifacts removal, the EEG alpha rhythms in the eyes-open and eyes-closed states were extracted. Figs. 7A and 7B shows the representative alpha rhythm extracted from EEG signals in the eyes-open and eyes-closed states, respectively. It is clear that the alpha rhythm in the eyes-open state was weaker than that in the eyes-closed state. Actually, the power value of the extracted alpha rhythm in the eyes-open state was 9.84  µV²; while in the eyes-closed state the power value was 50.52 µV². Figure 7C illustrates the alpha rhythm spectrogram, which is the square of the rhythm amplitude as a function of frequency. It presents an obvious difference between the eyes-open and eyes-closed states.

In order to verify the rhythms extraction using the adaptive SSA, another recently reported SSA method (Mohammadi et al., 2016) and wavelet decomposition method (Akar et al., 2015), which will be called SSA 2# and WDec respectively in this paper, were used for comparison. The comparison results are shown in Fig. 8. It can be seen from Figs. 8A and 8D that, in the eyes-open state, the extracted alpha rhythm using the adaptive SSA was of low amplitude and within the alpha band (8–13 Hz). Therefore, it could represent the real alpha rhythm of EEG. From Figs. 8B and 8D, it can be seen that the SSA 2# could not extract any alpha rhythm in the eyes-open state. It was because the SSA 2# could only extract the dominate component from EEG; while the alpha rhythm in the eyes-open state is weak and not the dominate component. Figures 8C and 7D both show that the amplitude of the extracted alpha rhythm using the WDec in the eyes-open state was even higher than that in the eyes-closed state. It was obviously inconsistent with reality. Besides, it contained a large number of components out of the alpha band. Therefore, in the eyes-open state, the adaptive SSA performed better than both the SSA 2# and WDec in rhythms extraction.

In the eyes-closed state, the extracted alpha rhythms using the adaptive SSA and SSA 2# were both of high amplitude and within the alpha band. However, the extracted alpha rhythm using the WDec contained components out of the alpha band. Therefore, in the eyes-closed state, the adaptive SSA performed as well as the SSA 2#, but better than the WDec, in rhythms extraction.