Abstract
Electroencephalogram (EEG) recordings are often contaminated with muscle artifacts. This disturbing muscular activity strongly affects the visual analysis of EEG and impairs the results of EEG signal processing such as brain connectivity analysis. If multichannel EEG recordings are available, then there exist a considerable range of methods which can remove or to some extent suppress the distorting effect of such artifacts. Yet to our knowledge, there is no existing means to remove muscle artifacts from single-channel EEG recordings. Moreover, considering the recently increasing need for biomedical signal processing in ambulatory situations, it is crucially important to develop single-channel techniques. In this work, we propose a simple, yet effective method to achieve the muscle artifact removal from single-channel EEG, by combining ensemble empirical mode decomposition (EEMD) with multiset canonical correlation analysis (MCCA). We demonstrate the performance of the proposed method through numerical simulations and application to real EEG recordings contaminated with muscle artifacts. The proposed method can successfully remove muscle artifacts without altering the recorded underlying EEG activity. It is a promising tool for real-world biomedical signal processing applications.
1. Introduction
The electroencephalogram (EEG) is frequently contaminated by various physiological activities of noninterest, such as electrocardiogram (ECG), electrooculogram (EOG), and electromyogram (EMG). These artifacts reduce the quality of the signal and blur features of interest. While ECG and EOG artifacts can be effectively removed by using adaptive filters and blind source separation (BSS) techniques [1], the perturbation induced by muscular activity (e.g., biting, chewing, and frowning) is particularly difficult to correct as recently reviewed in [2]. The main reason lies in the fact that EMG artifacts have higher amplitude (compared with the EEG signal), wide spectral distribution, and variable topographical distribution [2]. These muscle artifacts obscure EEG signals and make the interpretation of the EEG complicated or even unfeasible [3].
Low-pass filters are commonly employed to remove muscle artifacts. However, since the frequency spectrum of muscle artifacts significantly overlaps with that of brain signals, these filters not only suppress muscle artifacts but also interesting brain signals [4]. Recently, as one of the most popular BSS techniques, independent component analysis (ICA) has been extensively explored for this purpose [5–7]. ICA utilizes higher-order statistics and aims to separate the EEG recordings into statistically independent components (ICs). Clean EEG data can then be reconstructed by removing artifacts-related ICs from the raw EEG data. However, in some studies muscle artifacts seriously contaminate most ICs and crosstalk of brain and muscle activity can be observed [8, 9]. One possible reason is that ICA only exploits the spatial structure of source signals and the marginal distribution of the observations. Thus, it is suitable when source signals are temporally independent [10]. However, the artifacts typically have certain temporal structure which can be exploited for better source separation.
Second order blind identification (SOBI) takes temporal structure into consideration and simultaneously diagonalizes several covariance matrices at different time lags [10]. It has been shown that SOBI improved the performance significantly over ICA [11]. Yet, SOBI only considers stationary sources and it may suffer when there exist nonstationary ones, such as transient muscular activities [12]. More recently, a canonical correlation analysis (CCA) method has been proposed as a more suitable BSS approach for separating EMG artifacts from EEG [13]. Due to the broad frequency spectrum of the EMG artifacts, they resemble temporal white noise and thus have lower autocorrelation compared to EEG signals. The method exploits this characteristic for EMG removal and has been shown to outperform ICA on simulated data. Later, the result has been again documented by Gao et al. [14]. One possible reason for CCAs superior performance over ICA is due to the fact that muscle artifacts involve the movement of a group of muscles, which do not have a stereotyped topography [14]. Therefore, ICA does not function correctly here. Lately, a novel multiple time-lag CCA-based method was proposed to further improve the performance for removing muscular artifact in EEG [15]. However, it did not provide enough quantitative and comparative results based on either simulated or real data.
In recent years, biomedical signal measurement and processing are increasingly being deployed in ambulatory situations, particularly in healthcare applications. It is a trend to transit healthcare systems from hospital-centric toward ambulation-based, where minimal instrumentation and low computational complexity are required. The emerging wearable and portable wireless EEG system is a representative [16–18]. To reduce the complexity, quite a few ambulatory systems operate using only one single EEG channel, for example, [17, 18], in which case it is crucial to suppress muscle artifacts to extract as much valuable information as possible. However, almost all current methods for muscle artifact removal so far have been designed to handle multichannel/multidimensional datasets and will fail to isolate muscle activity in the current situation where only single-channel (unidimensional) EEG recordings are available.
In this study, we therefore propose a simple yet effective method to achieve the muscle artifact removal from single-channel EEG. Actually, it is a two-step modeling strategy. In the first step, unidimensional EEG is decomposed into multidimensional datasets. To implement this step, empirical mode decomposition (EMD) is a suitable option. EMD is a single-channel technique that decomposes nonstationary and nonlinear time series into a finite number of intrinsic mode functions (IMFs) [19]. Compared with other decomposition methods (e.g., wavelet transform), EMD is completely data-driven, meaning that it decomposes a signal in a natural way without requiring prior knowledge [20]. It has been shown to be efficient in many biomedical applications, for example, removing motion artifacts from functional near-infrared spectroscopy (fNIRS) data [21] and eliminating eye blink artifacts from EEG recordings [22]. However, the original EMD algorithm is highly sensitive to noise and may cause mode mixing. Recently, a noise assisted version of EMD, called ensemble EMD (EEMD), was proposed and has been demonstrated to be more robust in real-life applications [23]. In the second step, multiset canonical correlation analysis (MCCA) [24] is utilized as a BSS technique instead of the conventional CCA to the multidimensional datasets obtained from the first step. The advantage of MCCA over CCA will be discussed. The separation of muscle and brain activity components can then be achieved due to the relative low autocorrelation of muscle artifacts in comparison with brain activity. We denote the method as EEMD-MCCA by exploring the combination of EEMD and MCCA. The main contribution of the proposed method is to solve the practical muscle artifact removal problem from single-channel EEG, especially at the time when ambulatory healthcare is drawing continuously increasing attention.
We will examine the performance of the proposed EEMD-MCCA method on both synthetic and real datasets. We first validate it on simulated data by quantitative measures. We then apply it to real EEG recordings contaminated with muscle artifacts. We note that while EEMD-MCCA is proposed to remove muscle activity for the single-channel EEG case, it is generally applicable when one dataset contains relatively fewer channels (e.g., two or three) by first applying EEMD to each channel and then utilizing MCCA to the integrated signals after decomposition.
2. Materials and Methods
2.1. Methods
In this section, we first briefly introduce the techniques employed in this paper. Then we describe the new proposed EEMD-MCCA method.
Notations. Scalars are denoted by lowercase italic letters , vectors by lowercase boldface letters , matrices by boldface capitals , and the number of rows and columns by italic capitals . Matrix or vector transposition is denoted by an uppercase superscript (e.g., , ). The symbol (with size ) is used to represent the original single-channel signal. It can be also expressed like this , where is the value of the signal at the time point .
2.1.1. Ensemble Empirical Mode Decomposition
EMD is a single-channel decomposition method for nonstationary and nonlinear signals [19]. EMD decomposes a signal into a finite number of IMFs, which represent fast to slow oscillations. An IMF is a function that satisfies two following conditions [19]: (1) the number of extrema and the number of zero crossings must either be equal or differ at most by one; and (2) at any point, the mean value of the envelope defined by the local maxima and the envelope defined by the local minima are zero. To obtain an IMF from the original signal , a sifting process is performed [19]. First, all extrema of the original signal need to be identified. All local maximum points are connected by a cubic spline line as the upper envelope . Repeat the procedure for local minimum points to form the lower envelope . Their mean is calculated as The difference between the signal and the mean is defined as the first component as follows: In the second sifting process, is treated as the signal, and we can have Subsequently, we can repeat this sifting procedure times until is an IMF, with Therefore the first IMF component derived from the original signal is designated as
A criterion for stopping the sifting process to have an IMF has been established by limiting the size of the standard deviation (SD), calculated from the two consecutive sifting sequences as follows: A typical value for SD can be set between 0.2 and 0.3 [19].
To extract the 2nd IMF component, we remove from the original signal to have The residual is treated as a new signal and the same sifting process is applied to obtain the 2nd IMF component and the residual This procedure can be repeated on the subsequent residuals ’s until the final residual no longer contains any oscillation information: By summing up (7), (8), and (9), we can obtain Thus, we decompose the original signal into empirical modes ’s and a residue .
However, the original EMD algorithm is highly sensitive to noise. Recently, Wu and Huang introduced a new noise-assisted data analysis method, called EEMD [23]. The method defines the true IMF components as the mean of an ensemble of trials. Each trial consists of the signal plus an additive independent identically distributed white noise of the same standard deviation. In this case, although each individual trial may produce noisy results, it is canceled out in the ensemble mean of sufficient trials since the noise in each trial is assumed independently.
2.1.2. Canonical Correlation Analysis
Two zero-mean data sets are stored in two matrices, with size and with size , where means the number of observations and and indicate the numbers of variables in corresponding matrices. Conventional CCA is to find linear combinations of both and variables which have maximum correlation coefficient with each other [25]. This leads to the following objective function with constraints: where ’s are the weight vectors.
The solutions to this problem are the largest eigenvectors of the matrices and , respectively. The subsequent weights are the eigenvectors of the same matrix in the order of decreasing eigenvalues. The canonical variates () can be calculated directly from the original matrices ’s as . The corresponding rows between and are highly correlated, while the rows within each individual are uncorrelated with each other. The detailed derivation can be referred to [26].
Due to the aforementioned property, conventional CCA has been further extended to solve the BSS problem by assuming source components to be maximally autocorrelated and mutually uncorrelated in a functional magnetic resonance imaging (fMRI) study [27]. In the setting, let be the observed data matrix with mixtures and samples, and let be a temporally delayed version of the original data matrix . Thus, CCA can separate the recorded data into the self-correlated and mutually uncorrelated sources. As a potential alternative for the most widely used ICA method, CCA has been previously examined against a number of ICA algorithms using multichannel or single-channel recordings. The CCA-based methods were shown to outperform the ICA-based techniques for EEG/fNIRS artifact removal [13, 14, 21] and also demonstrated to be more computationally efficient when having similar qualitative results for EEG/fMRI source separation [27, 28] due to the usage of second-order statistics.
2.1.3. Multiset Canonical Correlation Analysis
MCCA extends the theory of CCA to more than two data sets to identify canonical variates that summarize the correlation structure among multiple random vectors by linear transformations. Unlike CCA where correlation between two canonical variates is maximized, MCCA aims to optimize an objective function of the correlation matrix of the canonical variates from multiple random vectors in order to make the canonical variates achieve the maximum overall correlation [24]. Suppose we have data sets with size . The aim of MCCA is to extract source components which are uncorrelated within each individual data set and meanwhile correlated well across the data sets. Analogously, it is straightforward to extend MCCA to handle the BSS problem under the similar assumption by letting
To demonstrate the advantage of MCCA over CCA, we will briefly discuss the source separability conditions. For more details, one can refer to [24]. When , the following condition must be satisfied to successfully recover the underlying sources by CCA: where represents the correlation coefficient between the th source from the 1st data set and the th source from the 2nd data set. When , the following requirement must be met to successfully recover the underlying sources from each of data sets correspondingly by MCCA: where represents the correlation coefficient between the th source from the th data set and the th source from the th data set.
It is important to note that condition (14) is more relaxed than (13), especially for our discussed BSS problem. More specifically, if two underlying sources have the same autocorrelation regarding to the time delay , condition (13) will not be met so that these two sources cannot be recovered successfully by CCA. However, as long as these two sources have different autocorrelations for one of the possible time delays (), they can be extracted completely by MCCA according to (14). The difference between (14) and (13) suggests that solving the BSS problem on a larger group of data sets is easier than doing that on a smaller group of datasets.
2.1.4. The Proposed EEMD-MCCA
To deal with the muscle artifact removal problem in single-channel EEG, we propose taking advantage of both EEMD and MCCA by exploring their combination and denote the proposed method as EEMD-MCCA. In fact, it is a two-step modeling strategy.
In the first step, EEMD is employed to decompose the single-channel EEG signal and derive a set of averaged IMFs. All the IMF components and the final residual are placed into a matrix . The size of is , where . Regarding the ensemble number , it is found that the performance of the technique becomes fairly consistent when using ten or more ensembles in our application. This is a proper number in practice considering the computational cost. The noise standard deviation has been suggested empirically to 0.2 times the standard deviation of the original signal [23].
In the second step, (12) is first used to generate temporally delayed versions of the matrix . Then MCCA is applied to the data sets and the underlying sources in are extracted and ordered in terms of autocorrelation from high to low. The sources with low autocorrelation correspond to muscle artifacts and can be removed by setting the corresponding row of the matrix to be zero. The source matrix is then passed through the mixing matrix to return the cleaned multichannel signals which are now, ideally, free of artifacts. The artifact-free single-channel recording can be determined by summing the recovered IMFs in the matrix . Regarding the parameter , we will discuss more in the simulation part. After these two steps, the muscle activity can almost be removed from single-channel EEG. The specific implementation procedure is shown in Algorithm 1.
| ||||||||||||||||||||||||||||||
2.2. Data Description
2.2.1. Synthetic Data
To demonstrate the performance of the proposed EEMD-MCCA method, in this section we will generate synthetic single-channel EEG with real-life muscle artifacts. Further, we employ some measures to test the performance since the ground truth is known.
Conventionally, the “ground truth” EEG signals without muscle artifacts were selected according to the visual inspection of experienced neurophysiologists. However, not only it is difficult to obtain clean EEG signals, but also there is no guarantee that the signals are completely free of muscle activity solely relying on visual inspection. Thus, in this study, we tend to use synthetic EEG data. A single-channel EEG series can be generated according to the phase-resetting theory [29, 30]. Similar to Mäkinen et al. [29], we generated our simulated data by summing 4 such sinusoids with frequencies chosen randomly from range 4–20 Hz. The sampling frequency was 250 Hz. Ten trials of EEG were generated and each trial was 1 second long. Then a 10-second series could be formed by concatenating the 10 trials, containing mainly theta, alpha, and beta activity. It should be noted that while each trial included 4 distinct frequencies, the frequencies chosen in different trials were also independent, which means that there was rich frequency information in the 10-second series.
To simulate real-life situations, obtaining pure muscle activity is quite necessary. It is insufficient to select muscle artifacts directly from the EEG as they contain both muscle and brain activity. To remove EEG activity and acquire muscle activity, ICA was utilized to decompose a real EEG data set with 21 channels. A neurophysiologist labeled the eye blink artifacts, eye movement artifacts, and muscle artifacts from all the decomposed ICs by inspecting some features such as the power spectral density and topography. It is important to note that a large number of ICs contained both EMG and ongoing EEG activity. Nevertheless, there existed one component containing pure EMG activity, denoted by . Since we focus on single-channel issues, it is unnecessary to reconstruct the component with the corresponding field distribution.
The EMG activity was superimposed on the EEG signal as follows: where represents the contribution of muscle activity. Figure 1 shows the original EEG signal and the EEG containing muscle artifacts (). The signal-to-noise ratio (SNR) can then be adjusted by changing the parameter : where the root mean squared (RMS) value is defined as The relative root-mean-squared error (RRMSE) is used as an evaluation measure to the effect of muscle artifact removal, which is defined as follows: where is the estimated EEG signal after muscle artifact removal. To further measure the capability of the proposed method for preserving the original EEG signal, correlation coefficients between the two waveforms and are also calculated. Hence, in this work, RRMSE and correlation coefficient (CC) serve as the main criteria for measuring the performance of muscle artifact removal.
(a)
(b)
2.2.2. Real Data
In the real data study, two datasets were employed to demonstrate the effect of the proposed method. One dataset was collected by us from eight health subjects during their stable cycling on an exercise bicycle. The study was approved by the University of British Columbia Ethics Board, and all subjects gave written, informed consent prior to participating. The EEG data were collected using an EEG cap (Quick-Cap, Compumedics, Texas, USA) with 9 electrodes based on the International 10–20 system, referenced to linked mastoids. The EEG data were sampled at 1000 Hz using SynAmps2 amplifiers (NeuroScan, Compumedics, Texas, USA). Data were later processed by a band-pass filter between 1~70 Hz. EEG recordings during exercise are easily contaminated with muscle artifacts and those artifacts can largely complicate the subsequent EEG signal processing such as brain connectivity analysis. Figure 2 displays one 10-second scalp EEG segment. All channels were more or less contaminated with muscle activity during the 10 seconds.
The other dataset is the public ictal EEG from the BioSource database established by Sabine Van Huffel (http://www.esat.kuleuven.be/stadius/members/biomed/biosource.htm). Ictal EEG is often severely contaminated with muscle artifacts, which make the determination and localization of the ictal onset complicated. Figure 3 shows the 10-second scalp EEG recordings with 21 channels from a long-term Epilepsy Monitoring Unit (OSG EEG recorders, Rumst, Belgium). Electrodes were placed according to the International 10–20 system with additional sphenoidal electrodes. The sampling frequency was 250 Hz. The EEG was digitally filtered by a band-pass filter between 0.3~35 Hz. A notch filter was applied to suppress the 50 Hz power-line interference. The seizure EEG was contaminated with muscle artifacts and eye blinks. Muscle artifacts can be observed between 0 s and 3.9 s on channels F7, T3, T5, C3, and T1 and between 5 s and 10 s on channels F8, T4, F4, C4, and P4.
Although the EEG recordings here are not based on single-channel, we can still apply the proposed EEMD-MCCA method to each channel individually and demonstrate its effectiveness for removing muscle artifacts from different regions in the brain.
3. Results and Discussion
3.1. The Synthetic Data Case
We applied the proposed EEMD-MCCA method to the synthetic single-channel data according to the procedure of Algorithm 1. As mentioned in Section 2.1.4, we tested the reliability of the method with different and SNR values in terms of RRMSE and CC as shown in Figure 4. Note that the number of data sets was chosen from 2 to 20 (i.e., ). When , it means that CCA was used. When , MCCA was employed. The parameter in (12) can be chosen empirically as it may highly depend on the data structure. In this simulation study, we found that similar results could be obtained by examining a set of values. One possible reason is that by applying MCCA to multiple time-delayed data sets it is sufficient to fully exploit the temporal structure of the original data set no matter what value is. Among the values, provided slightly better performance and was chosen here.
(a)
(b)
It can be seen from Figure 4 that MCCA can consistently achieve better performance than CCA, which is in accord with the theoretical analysis in Section 2.1.3. Another interesting observation is that when the performance becomes quite stable. This suggests that is a proper number of data sets in practical applications. To further demonstrate the practicability of the method, we tested the time cost when and , which were 2.61 seconds and 2.73 seconds averaged over 100 independent runs separately. The implementation was done in MATLAB (MathWorks Inc., Novi, MI, USA) and run under Microsoft Windows 8 x64 OS on the computer with Dual Intel(R) Core(TM) i-3427U 1.80 GHz CPU and 8.00 GB RAM. Considering the improved performance, the slightly increased time cost is well acceptable for removing artifacts from 10-second EEG data, especially for a number of ambulatory systems in which fast clean information and direct feedback are essentially important.
To see some details of the method, we also present the stepwise results in Figure 5. The IMF components extracted by EEMD were shown in Figure 5(a) from high to low frequencies. After applying MCCA, the uncorrelated sources were ordered in terms of their autocorrelations as displayed in Figure 5(b). The muscle activity was present in the last two components with lowest autocorrelations in the MCCA decomposition. Excluding the muscle artifact components in the reconstruction led to the cleaned EEG shown in Figure 5(c). To further illustrate the performance, an amplified version including both recovered and original EEG signals was present in Figure 5(d), from which we can see that the proposed method highly preserved the original brain activity.
(a)
(b)
(c)
(d)
In addition, we also implemented low-pass filter and ICA-based methods for a performance comparison study in terms of RRMSE and CC at various SNR values. Butterworth filter of order 8 was employed with three different cutoff frequencies equal to 10, 30, and 50 Hz. ICA was applied to the same IMF components extracted by EEMD as our proposed method did. It is termed EEMD-ICA. The joint approximate diagonalization of eigenmatrices (JADE) algorithm was adopted here for ICA implementation [31]. Muscle artifact components, after applying EEMD-ICA, had to be selected according to visual inspection. The cleaned signal was reconstructed by excluding the components related to the artifact. The results are shown in Figure 6, from which we can see that EEMD-MCCA has consistently better performance than the low-pass filters. The possible reason is due to the fact that the low-pass filters were insufficient to remove all artifacts without altering the underlying brain activity since the frequency spectrum of the muscle artifacts overlaps with that of the brain signal [13]. Another concern for the low-pass filters is the nonlinear phase-frequency response characteristic in practice, which will lead to signal distortion. It should be noted that we did not present the results of EEMD-ICA in Figure 6, because it was not capable of separating the muscle artifact from the EEG signal. Two examples of its decomposition are shown in Figure 7. We can easily observe that the EMG and EEG components were mixed together (e.g., IC9 in Figure 7(a) and IC1 in Figure 7(b)). It is quite difficult to determine which components should be excluded since we did not want to remove any brain activity.
(a)
(b)
(a)
(b)
3.2. The Real Data Case
In this case study, we applied the proposed EEMD-MCCA method to each individual channel of the EEG recordings as shown in Figures 2 and 3. We found that it was fairly easy to distinguish the muscle artifact components from the ones related to the brain activity. When processing each single-channel EEG recordings, muscle activity was almost present in the last several components in the MCCA decomposition. Excluding those components in the reconstruction of the EEG resulted in the cleaned EEG as shown in Figures 8 and 9. It can be seen that muscle artifacts have been sufficiently removed. In particular, for the ictal EEG, the ictal activity on the T2, F8, T4, and T6 electrodes was perfectly preserved. The ictal activity on F8 and T4, which originally was blurred by muscle artifacts, becomes visible by using the proposed EEMD-MCCA method. It should also be noted that there exist some obvious EOG artifacts in both real datasets, while their removal is beyond the scope of this paper. However, these EOG artifacts can help demonstrate the superior performance of our proposed method due to the fact that they were preserved with little distortion.
4. Conclusions
In this paper, we proposed a simple yet effective EEMD-MCCA method to realize muscle artifact removal from single-channel EEG. We illustrated the performances of the proposed method using both synthetic data and real-life data. We observed that the method is able to remove muscle activity effectively and efficiently and meanwhile preserve the brain activity very well. It is worth noting that while EEMD-MCCA is proposed to remove muscle activity for the single-channel EEG case, it is generally applicable when one dataset contains relatively fewer channels (e.g., two or three). The proposed method is a promising single-channel technique under the current situation that ambulatory healthcare systems are increasingly emerging.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
The work is supported by the National Natural Science Foundation of China (no. 61172037).