A flexible and accurate method for electroencephalography rhythms extraction based on circulant singular spectrum analysis
- Published
- Accepted
- Received
- Academic Editor
- Luciano Fadiga
- Subject Areas
- Bioengineering, Bioinformatics, Computational Biology, Neuroscience, Radiology and Medical Imaging
- Keywords
- Rhythms extraction, Electroencephalography, Circulant singular spectrum analysis, Eyes-open and eyes-closed condition
- Copyright
- © 2022 Hu et al.
- Licence
- This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, reproduction and adaptation in any medium and for any purpose provided that it is properly attributed. For attribution, the original author(s), title, publication source (PeerJ) and either DOI or URL of the article must be cited.
- Cite this article
- 2022. A flexible and accurate method for electroencephalography rhythms extraction based on circulant singular spectrum analysis. PeerJ 10:e13096 https://doi.org/10.7717/peerj.13096
Abstract
Rhythms extraction from electroencephalography (EEG) signals can be used to monitor the physiological and pathological states of the brain and has attracted much attention in recent studies. A flexible and accurate method for EEG rhythms extraction was proposed by incorporating a novel circulant singular spectrum analysis (CiSSA). The EEG signals are decomposed into the sum of a set of orthogonal reconstructed components (RCs) at known frequencies. The frequency bandwidth of each RC is limited to a particular brain rhythm band, with no frequency mixing between different RCs. The RCs are then grouped flexibly to extract the desired EEG rhythms based on the known frequencies. The extracted brain rhythms are accurate and no mixed components of other rhythms or artifacts are included. Simulated EEG data based on the Markov Process Amplitude EEG model and experimental EEG data in the eyes-open and eyes-closed states were used to verify the CiSSA-based method. The results showed that the CiSSA-based method is flexible in alpha rhythms extraction and has a higher accuracy in distinguishing between the eyes-open and eyes-closed states, compared with the basic SSA method, the wavelet decomposition method, and the finite impulse response filtering method.
Introduction
Electroencephalograms (EEGs) are the electrical activity of the brain’s neurons recorded at the scalp surface (Henry, 2006). They consist of several rhythm bands: delta (1–4 Hz), theta (4–8 Hz), alpha (8–13 Hz), beta (13–30 Hz) and gamma (>30 Hz). Because the rhythms reflect different physiological and pathological information, EEG rhythms extraction has been widely applied in many areas. Examples include portable and wearable EEG devices (Hwang et al., 2018; Maskeliunas et al., 2016), mental fatigue assessment (Taran & Bajaj, 2017), disease diagnosis (Babiloni et al., 2016; Gupta & Pachori, 2019), and brain computer interface systems (Jeunet et al., 2019; Liu et al., 2020).
The accuracy of EEG rhythms extraction determines the physiological and pathological information it provides. Various methods have been proposed to extract the desired EEG rhythms. Filtering components have the ability to restrict a signal to a specific frequency band, and such bandpass filters were first used to extract EEG rhythms (Pfurtscheller et al., 1997). This method performed well in EEGs of high signal-to-noise ratio (SNR). Then, the wavelet transform (WT) method was used for EEG rhythms extraction (Duque-Muñoz, Pinzon-Morales & Castellanos-Dominguez, 2015). By estimating the rhythms with a customized wavelet, the WT method can extract time-varying EEG rhythms with changes in brain state. To facilitate the EEG rhythms extraction, the independent component analysis (ICA) method was then introduced (Kavuri, Veluvolu & Chai, 2018). By incorporating priori information about the desired rhythms as reference signals, the ICA method can extract EEG rhythms automatically. However, the extracted rhythms using the bandpass filter, WT, and ICA methods were contaminated by noise and artifacts overlapping in time–frequency space. In recent years, to improve the accuracy of EEG rhythms extraction, the singular spectrum analysis (SSA) method has been used (Akar et al., 2015; Mohammadi et al., 2016). This nonparametric method enables the separation of different sources even when they overlap in time–frequency space (Mohammadi et al., 2015).
The SSA method is a nonparametric signal processing method proposed to extract useful information from experimental data (Broomhead & King, 1986). SSA consists of two stages: decomposition and reconstruction. Decomposition involves time-delay embedding called Takens’ theorem (Takens, 1981), followed by singular value decomposition (SVD) (Mees, Rapp & Jennings, 1987). Reconstruction involves grouping and diagonal averaging (Vautard, Yiou & Ghil, 1992). SSA works well for single channel signals as well as multichannel signals. It allows the exploitation of SSA for a variety of applications including biomedical signal processing such as signal restoration, change detection, segmentation, anomaly detection and prediction (Sanei & Hassani, 2015; Xu et al., 2018). Compared with the narrowband filter methods which are not able to separate the component mixing overlapping in frequency space, the SSA method is utilized to find the structure of nonlinear and non-stationary signal and enables the separation of different sources even overlapping in the time-frequency space. Therefore, SSA has been successfully applied in the research of artifacts removal and rhythms extraction from the EEG signal (Maddirala & Shaik, 2016a, 2016b; Mohammadi et al., 2015; Xu et al., 2018).
In the basic SSA method, the grouping rule is important for SSA reconstruction. However, because of the lack of the information about the amplitude and frequency of the reconstructed components (RCs), there is no general grouping rule. Different grouping rules have been proposed depending on the research target, the types of signals, and noise. The conventional SSA grouping is performed according to the magnitudes of eigenvalues related to the power of each RC (De Carvalho & Rua, 2017; Teixeira et al., 2005; Yang et al., 2016). Mohammadi et al. (2016) proposed a new grouping rule based on eigenvalue pairs to extract the main rhythms from sleep EEG signals. Hai et al. (2017) proposed another efficient grouping rule based on the similarity between the eigenvalues and the peak frequency of RC, which makes SSA adaptive to EEG signals containing different levels of artifacts and rhythms. However, these grouping rules can only be applied to specific types of signals and must be incorporated with other methods (e.g., Fourier transform or wavelet decomposition) to pre-identify the frequencies of RCs, which is time-consuming and inflexible. Besides, the observed frequency mixing between different RCs leads to inaccurate EEG rhythms extraction (Xu et al., 2018).
In this paper, we introduce a novel circulant singular spectrum analysis (CiSSA) method (Bógalo, Poncela & Senra, 2021) on the basis of SSA to improve the flexibility and accuracy of EEG rhythms extraction. Compared with the basic SSA method, the CiSSA method has the advantage of avoiding the need for pre-identifying the frequencies of RCs. A set of orthogonal vectors are obtained by decomposing the circulant matrix, and the EEG signals can be decomposed into the sum of a set of orthogonal RCs of known frequencies. The RCs can be grouped automatically and flexibly to extract the specific EEG rhythms based on their frequencies. In addition, because the frequency bandwidth of each RC is limited to a particular band of the brain rhythm of interest and the artifacts are removed by the basic SSA, the extracted brain rhythms are accurate and no mixed components of other rhythms and artifacts are included.
Methods
The CiSSA method is a nonparametric signal extraction method proposed by Bógalo, Poncela & Senra (2021). CiSSA consists of four steps: embedding, decomposition, diagonal averaging, and grouping. As in the basic SSA method, in the time-delay embedding step, the single-channel EEG time series $\mathbf{s}=({s}_{1},{s}_{2},\dots ,{s}_{N}{)}^{T}$ (superscript T denotes the transpose of a vector) is mapped onto a multidimensional trajectory matrix X using a sliding window (Takens, 1981):
(1) $$\mathbf{X}=({\mathbf{S}}_{1},{\mathbf{S}}_{2},\dots ,{\mathbf{S}}_{K})=\left(\begin{array}{cc}\begin{array}{cc}{s}_{1}& {s}_{2}\\ {s}_{2}& {s}_{3}\end{array}& \begin{array}{cc}\cdots & {s}_{K}\\ \cdots & {s}_{K+1}\end{array}\\ \begin{array}{cc}\vdots & \vdots \\ {s}_{L}& {s}_{L+1}\end{array}& \begin{array}{cc}\ddots & \vdots \\ \cdots & {s}_{N}\end{array}\end{array}\right)$$
where L denotes the window length (or embedding dimension), K = N – L + 1, and ${\mathbf{S}}_{i}$ denotes the lagged vector.
In the decomposition step, the trajectory matrix is decomposed into elementary matrices of rank 1 that are associated with different frequencies. To do so, a related circulant matrix ${\mathbf{C}}_{L}$ is built based on the second order moments of the time series (Bógalo, Poncela & Senra, 2021):
(2) $${\mathbf{C}}_{L}(f)=\left(\begin{array}{cc}\begin{array}{ccc}{c}_{0}& {c}_{1}& {c}_{2}\\ {c}_{L-1}& {c}_{0}& {c}_{1}\end{array}& \begin{array}{cc}\cdots & {c}_{L-1}\\ \cdots & {c}_{L-2}\end{array}\\ \begin{array}{ccc}\vdots & \vdots & \vdots \\ {c}_{1}& {c}_{2}& {c}_{3}\end{array}& \begin{array}{cc}\ddots & \vdots \\ \cdots & {c}_{0}\end{array}\end{array}\right)$$
where
(3) $${c}_{m}=\frac{L-m}{L}{\gamma}_{m}+\frac{m}{L}{\gamma}_{L-m},\phantom{\rule{0.056em}{0ex}}{\gamma}_{m}=\frac{1}{N-m}\sum _{t=1}^{T-m}{s}_{t}{s}_{t+m},\phantom{\rule{1em}{0ex}}m=0,1,\dots ,L-1$$
The eigenvalues and eigenvectors of ${\mathbf{C}}_{L}$, respectively, are given (Gray, 2006)
(4) $$\begin{array}{c}{\lambda}_{k}=\sum _{m=0}^{L-1}{c}_{m}\mathrm{exp}\left(i2\pi m\frac{k-1}{L}\right)=f\left(\frac{k-1}{L}\right),\phantom{\rule{0.056em}{0ex}}k=1,2,\dots ,L\hfill \\ {\mathbf{u}}_{k}={L}^{-1/2}{\left({u}_{k,1},{u}_{k,2},\dots ,{u}_{k,L}\right)}^{H},{u}_{k,j}=\mathrm{exp}\left(-i2\pi \left(j-1\right)\frac{k-1}{L}\right),\phantom{\rule{0.056em}{0ex}}k=1,2,\dots ,L\hfill \end{array}$$
where $f(\cdot )$ denotes the power spectral density of the signal, and H indicates the conjugate transpose of a matrix. The k-th eigenvalue and the corresponding eigenvector are associated with the given frequencies by
(5) $${f}_{k}=\frac{k-1}{L}{f}_{s}$$
where ${f}_{s}$ is the sampling rate of the EEG signals. As a consequence, the diagonalization of ${\mathbf{C}}_{L}$ allows us to write X as the sum of the elementary matrices ${\mathbf{X}}_{k}$:
(6) $$\mathbf{X}=\sum _{k=1}^{L}{\mathbf{X}}_{k}=\sum _{k=1}^{L}{\mathbf{u}}_{k}{\mathbf{u}}_{k}^{H}\mathbf{X}$$
The symmetry of the power spectral density leads to ${\lambda}_{k}={\lambda}_{\mathrm{L}+2-k}$. The corresponding eigenvectors given by Eq. (4) are complex; therefore, they are paired with complex conjugates, ${\mathbf{u}}_{k}={\mathbf{u}}_{L+2-k}^{\ast}$ where ${\mathbf{u}}^{\ast}$ indicates the complex conjugate of a vector $\mathbf{u}$. Then, ${\mathbf{X}}_{k}$ and ${\mathbf{X}}_{L+2-k}$ correspond to the same frequency ${f}_{k}$.
To obtain the elementary matrices by frequency, we first form the groups of two elements ${B}_{k}=\left\{k,L+2-k\right\},k=2,3,\dots ,M,M=\lfloor \left(L+1\right)/2\rfloor $, $\lfloor \cdot \rfloor $ is the integer part operator, with ${B}_{1}=\left\{1\right\}$ and ${B}_{L/2+1}=\left\{L/2+1\right\}$ if L is even. Then, the real elementary matrix ${\mathbf{X}}_{{B}_{k}}$ is computed as the sum of the two elementary matrices ${\mathbf{X}}_{k}$ and ${\mathbf{X}}_{L+2-k}$, which are associated with eigenvalues ${\lambda}_{k}$, ${\lambda}_{\mathrm{L}+2-k}$ and frequency ${f}_{k}$, given by Eq. (5)
(7) $$\begin{array}{ll}{X}_{{B}_{k}}& ={X}_{k}+{X}_{L+2-k}={u}_{k}{u}_{k}^{H}X+{u}_{L+2-k}{u}_{L+2-k}^{H}X=({u}_{k}{u}_{k}^{H}+{u}_{k}^{*}{u}_{k}^{\prime})X\\ & =2({R}_{{u}_{k}}{R}_{{u}_{k}}^{\prime}+{I}_{{u}_{k}}{I}_{{u}_{k}}^{\prime})X\end{array}$$
where ${R}_{{\mathbf{u}}_{k}}$ and ${I}_{{\mathbf{u}}_{k}}$ denote the real and imaginary parts of ${\mathbf{u}}_{k}$, respectively, and the matrices ${\mathbf{X}}_{{B}_{k}}$ are real.
Then, in the diagonal averaging step (Vautard, Yiou & Ghil, 1992), several time series are reconstructed from the corresponding real elementary matrices ${\mathbf{X}}_{{B}_{k}}$. The reconstructed time series are generally called RCs. Theoretically, the frequencies of the RCs are given by Eq. (5). Finally, the alpha rhythm (8–13 Hz) can be extracted automatically by
(8) $${V}_{alpha}=\sum _{i=\lceil 1+8L/{f}_{s}\rceil}^{\lfloor 1+13L/{f}_{s}\rfloor}R{C}_{i}$$
The frequency bandwidth of each RC can be roughly expressed by Bozzo, Carniel & Fasino (2010) and Xu et al. (2018)
(9) $${f}_{b}={f}_{s}/L$$
As a consequence, the frequency bandwidth of each RC is limited to ${f}_{s}/L$. Considering the frequency of each RC given by Eq. (5), there is no frequency mixing between different RCs, and the extracted alpha rhythms do not contain mixed components of other rhythms. In order to remove the frequency components outside ${f}_{band}$, with ${f}_{band}$ denoting the bandwidth of the brain rhythm of interest (alpha rhythm), the frequency bandwidth ${f}_{b}$ is less than the bandwidth of the brain rhythm of interest ${f}_{band}$ ( ${f}_{b}\le {f}_{band}$). Therefore the embedding dimension L should satisfy the condition $L\ge {f}_{s}/{f}_{band}$.
The pseudo-code of the CiSSA method is shown in Algorithm 1. The Institutional Review Board of Tsinghua University granted ethical approval (20170010) to carry out the experiment within its facilities.
input s, L: single-channel EEG time series s and embedding dimension L |
output α: extracted alpha rhythm |
procedures |
(1) X: the trajectory matrix is constructed by Eq. (1) |
(2) C_{L}: the circulant matrix is built by Eqs. (2) and (3) |
(3) λ_{k}, u_{k}: the circulant matrix C_{L} is decomposed, and a set of eigenvalues and eigenvectors is derived by Eq. (4) |
(4) X_{Bk}: the real elementary matrices are derived by Eq. (7) |
(5) α: the RCs are grouped by Eq. (8) to obtain the alpha rhythm |
return α |
Simulation Results and Discussion
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)\mathrm{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}{\mathrm{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.
Symbol | Value | Comments |
---|---|---|
${f}_{1}$ (Hz) | 3.71 | Delta rhythm |
${\sigma}_{1}^{\xi}$ | 3.53 | |
${\gamma}_{1}$ | 0.98 | |
${f}_{2}$ (Hz) | 7.62 | Theta rhythm |
${\sigma}_{2}^{\xi}$ | 4.35 | |
${\gamma}_{2}$ | 0.95 | |
${f}_{3}$ (Hz) | 10.45 | Alpha rhythm |
${\sigma}_{3}^{\xi}$ | 1.65 | |
${\gamma}_{3}$ | 0.99 | |
${f}_{4}$ (Hz) | 15.43 | Beta rhythm |
${\sigma}_{4}^{\xi}$ | 0.24 | |
${\gamma}_{4}$ | 0.99 |
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\phantom{\rule{0.056em}{0ex}}\text{}\mathrm{H}\mathrm{z}$ and the bandwidth of alpha rhythm (8–13 Hz) ${f}_{band}=5\text{}\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\text{}\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\text{}\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\text{}\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\text{}\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\text{}\mu {V}^{2}/\mathrm{H}\mathrm{z}$ and $\epsilon =0.57\text{}\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\text{}\mu {V}^{2}/\mathrm{H}\mathrm{z}$ and $\epsilon =0.08\text{}\mu {V}^{2}/\mathrm{H}\mathrm{z}$, respectively, which were much lower than those obtained by the basic SSA, FIR and WDec methods.
Method | $\overline{\mathbf{\text{\epsilon}}}\left(\mathbf{\text{\mu}}{\mathit{V}}^{\mathbf{2}}/\mathbf{H}\mathbf{z}\right)$ |
---|---|
CiSSA (L = 80) | $0.29\pm 0.08$ |
basic SSA (L = 80) | $0.40\pm 0.18$ |
FIR | $0.49\pm 0.23$ |
WDec | $0.58\pm 0.21$ |
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.
L | $\overline{\mathbf{\text{\epsilon}}}\left(\mathbf{\text{\mu}}{\mathit{V}}^{\mathbf{2}}/\mathbf{H}\mathbf{z}\right)$ | L | $\overline{\mathbf{\text{\epsilon}}}\left(\mathbf{\text{\mu}}{\mathit{V}}^{\mathbf{2}}/\mathbf{H}\mathbf{z}\right)$ | L | $\overline{\mathbf{\text{\epsilon}}}\left(\mathbf{\text{\mu}}{\mathit{V}}^{\mathbf{2}}/\mathbf{H}\mathbf{z}\right)$ | L | $\overline{\mathbf{\text{\epsilon}}}\left(\mathbf{\text{\mu}}{\mathit{V}}^{\mathbf{2}}/\mathbf{H}\mathbf{z}\right)$ |
---|---|---|---|---|---|---|---|
20 | $0.82\pm 0.15$ | 70 | $0.56\pm 0.14$ | 120 | $0.50\pm 0.12$ | 170 | $0.50\pm 0.11$ |
30 | $0.89\pm 0.40$ | 80 | $0.29\pm 0.08$ | 130 | $0.33\pm 0.09$ | 180 | $0.38\pm 0.10$ |
40 | $0.39\pm 0.09$ | 90 | $0.44\pm 0.12$ | 140 | $0.45\pm 0.10$ | 190 | $0.45\pm 0.10$ |
50 | $0.76\pm 0.19$ | 100 | $0.59\pm 0.14$ | 150 | $0.54\pm 0.12$ | 200 | $0.54\pm 0.11$ |
60 | $0.33\pm 0.09$ | 110 | $0.38\pm 0.10$ | 160 | $0.40\pm 0.10$ |
Experimental Results and Discussion
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.
Subject # | CiSSA (L = 80) (%) | Basic SSA (L = 80) (%) | FIR (%) | WDec (%) |
---|---|---|---|---|
Subject 1 | 80.22 | 79.12 | 79.12 | 79.12 |
Subject 2 | 95.60 | 67.03 | 94.51 | 83.52 |
Subject 3 | 96.70 | 58.24 | 97.80 | 95.60 |
Subject 4 | 95.60 | 94.51 | 92.31 | 94.51 |
Subject 5 | 98.90 | 95.60 | 98.90 | 98.90 |
Subject 7 | 96.70 | 96.70 | 93.41 | 94.51 |
Subject 8 | 100 | 100 | 100 | 100 |
Subject 9 | 81.32 | 71.43 | 78.02 | 74.73 |
Subject 10 | 100 | 100.00 | 98.90 | 100 |
Subject 11 | 96.70 | 95.60 | 95.60 | 95.60 |
Subject 12 | 76.92 | 70.33 | 74.73 | 72.53 |
Subject 13 | 100 | 93.41 | 98.90 | 100 |
Subject 14 | 87.91 | 87.91 | 87.91 | 87.91 |
Subject 15 | 95.60 | 93.41 | 93.41 | 96.70 |
Subject 16 | 86.81 | 86.81 | 86.81 | 86.81 |
Subject 17 | 92.31 | 89.01 | 91.21 | 90.11 |
Subject 18 | 95.60 | 91.21 | 98.90 | 97.80 |
Subject 19 | 83.52 | 75.82 | 84.62 | 81.32 |
Subject 20 | 98.90 | 100.00 | 98.90 | 98.90 |
Subject 21 | 94.51 | 94.51 | 93.41 | 94.51 |
Subject 22 | 85.71 | 83.52 | 85.71 | 85.71 |
Average | 92.36 | 86.87 | 91.58 | 90.89 |
STD | 7.05 | 11.78 | 7.42 | 8.39 |
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.
Subject # | CiSSA (L = 80) (%) | Basic SSA (L = 80) (%) | FIR (%) | WDec (%) |
---|---|---|---|---|
Subject 1 | 91.23 | 87.72 | 89.47 | 88.60 |
Subject 2 | 92.98 | 93.86 | 95.61 | 95.61 |
Subject 3 | 95.61 | 91.23 | 95.61 | 95.61 |
Subject 4 | 93.86 | 80.7 | 88.6 | 84.21 |
Subject 5 | 92.11 | 84.21 | 84.21 | 96.49 |
Subject 6 | 86.84 | 86.84 | 84.21 | 89.47 |
Subject 7 | 96.49 | 85.09 | 93.86 | 94.74 |
Subject 8 | 97.37 | 98.25 | 99.12 | 99.12 |
Subject 9 | 80.7 | 77.19 | 79.82 | 66.67 |
Subject 10 | 93.86 | 89.47 | 91.23 | 77.19 |
Average | 92.11 | 87.46 | 90.17 | 88.77 |
STD | 4.74 | 5.86 | 5.79 | 9.69 |
Conclusions
In this paper, a flexible and accurate method based on CiSSA was proposed for alpha rhythm extraction from EEG signals. By decomposing the EEG signals into a set of orthogonal reconstructed components (RCs) at specific bandwidths of frequencies, the alpha rhythm can be extracted flexibly and accurately from EEG signals. The proposed method performed well on both simulated EEG data generated from the MPA EEG model and experimental EEG data, as well as the EEG data obtained from a public database. Features of the alpha rhythms extracted from experimental EEG signals were calculated to distinguish between the eyes-open and eyes-closed states. The CiSSA-based method showed higher classification accuracy and robustness than that of the basic SSA, FIR and WDec methods.
Supplemental Information
Experimental EEG data.
Experimental EEG data are in the .mat format (Matlab files). The link of the database EEG data is https://dataverse.tdl.org/dataset.xhtml?persistentId=doi:10.18738/T8/SS2NHB