Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Automatic Classification of Artifactual ICA-Components for Artifact Removal in EEG Signals

Automatic Classification of Artifactual ICA-Components for Artifact Removal in EEG Signals Background: Artifacts contained in EEG recordings hamper both, the visual interpretation by experts as well as the algorithmic processing and analysis (e.g. for Brain-Computer Interfaces (BCI) or for Mental State Monitoring). While hand-optimized selection of source components derived from Independent Component Analysis (ICA) to clean EEG data is widespread, the field could greatly profit from automated solutions based on Machine Learning methods. Existing ICA-based removal strategies depend on explicit recordings of an individual’s artifacts or have not been shown to reliably identify muscle artifacts. Methods: We propose an automatic method for the classification of general artifactual source components. They are estimated by TDSEP, an ICA method that takes temporal correlations into account. The linear classifier is based on an optimized feature subset determined by a Linear Programming Machine (LPM). The subset is composed of features from the frequency-, the spatial- and temporal domain. A subject independent classifier was trained on 640 TDSEP components (reaction time (RT) study, n = 12) that were hand labeled by experts as artifactual or brain sources and tested on 1080 new components of RT data of the same study. Generalization was tested on new data from two studies (auditory Event Related Potential (ERP) paradigm, n = 18; motor imagery BCI paradigm, n = 80) that used data with different channel setups and from new subjects. Results: Based on six features only, the optimized linear classifier performed on level with the inter-expert disagreement (<10% Mean Squared Error (MSE)) on the RT data. On data of the auditory ERP study, the same pre- calculated classifier generalized well and achieved 15% MSE. On data of the motor imagery paradigm, we demonstrate that the discriminant information used for BCI is preserved when removing up to 60% of the most artifactual source components. Conclusions: We propose a universal and efficient classifier of ICA components for the subject independent removal of artifacts from EEG data. Based on linear methods, it is applicable for different electrode placements and supports the introspection of results. Trained on expert ratings of large data sets, it is not restricted to the detection of eye- and muscle artifacts. Its performance and generalization ability is demonstrated on data of different EEG studies. Background artifactual sources. Typical artifacts of the EEG are Signals of the electroencephalogram (EEG) can reflect caused either by the non-neural physiological activities the electrical background activity of the brain as well as of the subject or by external technical sources. Eye the activity which is specific for a cognitive task during blinks, eye movements, muscle activity in the vicinity of the head (e.g. face muscles, jaws, tongue, neck), heart an experiment. As the electrical field generated by neural activity is very small, it can only be recognized beat, pulse and Mayer waves are examples for physiolo- by EEG if large assemblies of neurons show a similar gical artifact sources, while swaying cables in the mag- behavior. Resulting neural EEG signals are in the range netic field of the earth, line humming, power supplies or of micro volts only and can easily be masked by transformers can be the cause of technical artifacts. Brain-Computer Interfaces (BCI) are based on the sin- * Correspondence: irene.winkler@tu-berlin.de gle trial classification of the ongoing EEG signal and can Machine Learning Laboratory, Berlin Institute of Technology, Franklinstr. 28/ improve the life quality of disabled individuals especially 29, 10587 Berlin, Germany © 2011 Winkler et al; licensee BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 2 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 in combination with other assitive technology [1]. The independent components, stationarity of the sources and exclusion of artifacts is of special interest for BCI appli- the mixture, and prior knowledge about the number of cations, as the intended or unconscious use of artifacts components), their application usually leads to a good for BCI control are usually not desirable when the BCI separation, with only a small number of hybrid compo- system is tested on healthy subjects. Furthermore, as nents that contain both, artifacts and neural signals averaging methods have to be avoided, these real-time [9-12]. Existing methods for artifact rejection can be sepa- systems BCIs rely on relatively clean EEG signals. The rated into hand-optimized, semi-automatic and fully same holds true for other Mental State Monitoring automatic approaches. Semi-automatic approaches applications, that monitor a subject’s mental state con- tinuouslyand onafinegranulartimeresolutionto require user interaction for ambiguous or outlier com- detect changes e.g. of wakefulness, responsiveness or ponents [13,14]. While fully automated methods were mental workload as early as possible [2]. proposed for the classification of eye artifacts [15,16], The two physiological artifacts most problematic for these methods do not easily generalize to non-eye arti- BCI applications are ocular (EOG) and muscle (EMG) facts or even require the additional recording of the artifacts. EOG activity is either caused by rolling of the EOG [17,18]. Viola et al. and Mognon et al. [19,20] both eyes or by eye blinks which occur approx. 20 times per developed an EEGLAB plug-in which finds artifactual minute [3]. Both result in a low-frequency activity most independent components. Both plug-ins have a fully prominent over the anterior head regions, with maximal automatic mode that has been shown to recognize and frequencies below 4 Hz. In contrast, EMG activity reject major artifacts like eye blinks, eye movements and (caused by chewing, swallowing, head or tongue move- heart beats, while the detection of muscular or more ments) is usually a high-frequency activity (>20 Hz) subtle artifacts has not been reported. The plug-in which ranges from rather small to very large amplitudes developed by Viola et al. relies on a user-defined tem- [4]. plate, while Mognon’s approach does not require user For an extensive review of artifact reduction techni- interaction. ques in the context of BCI-systems, the reader can refer Existing more flexible approaches for the general clas- to Fatourechi et al. [5]. Since the rejection of artifactual sification of different artifact types were reported for trials amounts to a considerable loss of data, a method EEG data of epileptic patients [21], where the authors that removes the artifacts while preserving the underly- report a Mean Squared Error (MSE) of approx. 20% for their system based on a Bayesian classifier. Halder et al. ing neural activity is needed. For example, linear filter- [22] report a classification error below 10% for their ing is a simple and effective method if artifactual and neural activity are located in non-overlapping frequency Support Vector Machine (SVM) based system for a bands. Unfortunately, artifacts and the brain signal of fixed number of electrodes if dedicated artifact record- interest do usually overlap. Nevertheless, ocular activity ings are available for the classifier training. But even if can be partially removed by regression-based methods, such optimized conditions are present, difficulties of which subtract a part of the activity measured at addi- separating muscle artifact components from neural com- tional electrooculogram (EOG) channels from the EEG ponents are common [22]. (see [6] for a review). Regression-based methods require The review of the existing literature did not reveal a the reliable recording of additional EOG channels and systematic screening of potentially discriminant features are limited by the fact that the EOG is contaminated by for the general task of artifact detection/removal. More- brain activity which is removed as well. Furthermore, over, most approaches restrict themselves to part of the they cannot eliminate non-eye activity. available information, e.g. rely on spatial patterns only If artifactual signal components and neural activity of [19], or spatial patterns and spectral features [22], or interest are not systematically co-activated due to a dis- spatial pattern and temporal features [20]. advantageous experimental design, methods of Blind Our proposed solution for a general artifact detection Source Separation (BSS) like Independent Component method is motivated by the needs of EEG practitioners. Analysis (ICA) are promising approaches for their First, it is desirable that a method efficiently and reliably separation [7,8]. A common approach is the transforma- detects all classes of artifacts,e.g.isnotrestrictedto tion of the EEG signals into a space of independent eye-, heart beat-, or muscle artifacts. Second, a practical source components, the hand-selection of non-artifac- method must be applicable post-hoc, i.e. without the need of dedicated artifact recordings at the time of the tual neural sources and the reconstruction of the EEG experiment. Third, it is difficult to convince EEG practi- without the artifactual components (for an example of tioners to use a method of artifact rejection if it is a independent source components, see Figure 1). While assumptions for the application of ICA methods are black box and refuses introspection. As the goal must only approximately met in practice (linear mixture of be to develop a method, that delivers interpretable and Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 3 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 (a) (b) (c) sec Hz Filter P atte rn Figure 1 Three example independent source components. Time series (first column), spectrum (second column), filter (third column) and pattern (fourth column) of three components. The first row (a) shows an alpha generator in the occipital lobe. The second row (b) shows a rare muscle artifact component with an increased spectrum in higher frequencies. The third line (c) shows an eye artifact component that appears regularly, has an increases spectrum in lower frequencies and a typical front-back distribution in the pattern. easy to understand results, we decided for a linear clas- signal pre-processing methods (including a temporal sification method. Luckily, linear methods have proven a variant of ICA) are introduced, we describe 38 features high performance for a number of classification tasks in that are candidates for the artifact discrimination task. the field of EEG-based BCI systems. However, to be able BasedonlabelsprovidedbyEEGexperts, athorough to estimate the performance loss compared to a poten- feature selection procedure is described, that is used to tially better, but difficult to interpret, non-linear classifi- condense the 38 features to a small subset. Furthermore, cation method, the results of a Gaussian SVM are classification methods are introduced. The methods sec- reported in parallel. tion ends with a description of two other EEG para- We decided to use a sparse approach (sparsity in the digms (auditory Event Related Potential (ERP) and features) although it is a mixed blessing. It leads to a motor imagery for BCI), that will be used to validate the generalization approach of the proposed artifact classi- trade-off between efficiency and interpretability, as redundant but slightly less discriminative features are fier. In the results section, the outcome of the feature removed with high probability from the overall set of selection procedure is given, together with the artifact features. This has to be kept in mind during the analysis classification performance on unseen data of the RT of results. To reach the goal of a sparse method that paradigm, data of a unseen auditory ERP paradigm. delivers physiologically interpretable results, we decided Finally the method is applied in the context of a motor to incorporate a thorough feature selection procedure in imagery BCI setup, before the paper closes with a combination with a linear classification method that is discussion. based on features of all three available information domains of EEG data: the spatial domain (e.g. patterns Methods of independent components), the frequency domain and In the following subsections, we will describe how the the temporal domain. proposed new artifact classification method is set up. The paper is organized as follows: In the methods sec- Then we will introduce two further studies that are uti- tion, a reaction time (RT) paradigm is introduced, as lized to test the classifier’s generalizability. data from this study forms the basis for the construction Participants of the studies described below provided of the proposed artifact detection method. After the verbal and written informedconsent andwere freeto dB Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 4 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 stop their participation at any time. All collected data of a given multivariate time series X =(x ,..., x ) into 1 k was anonymized before any subsequent analysis or pre- unknown, assumed mutually independent source com- sentation took place. ponents S =(s ,..., s ) . Note that both the demixing W 1 k and the source components S are unknown, and that Classifier Construction using a RT study BSS algorithms differ in the definition of independence The artifact classifier is set up based on labeled indepen- between components. While ICA algorithms exploit dent components gained from a reaction time (RT) higher order statistics, TDSEP relies on second-order study. statistics by taking the temporal structure of the time Experimental Setup series into account. TDSEP amounts to finding a demix- Data from 12 healthy right-handed male subjects were ing W which leads to minimal cross-covariances over used to train and to test the proposed automated com- several time-lags between all pairs of components of S. ponent classification method. Every subject participated For a mathematical discussion, let in one EEG recording session of approx. 5 hours dura- be the cross-covariance (τ):= E(X (t)X (t − τ )) tion. EEG was recorded from 121 approx. equidistant matrix of the whitened data X at time-lag τ,where the sensors and high pass filtered at 2 Hz. During this ses- whitening transformation linearly decorrellates and sion, 4 repeated blocks of 3 different conditions (C0, C1, scales the data such that Σ(0) = I. Consider now that (1) C2) were performed. Each block lasted approx. 45 min- Whitening reduces the BSS problem to finding an utes. During all three conditions, subjects performed a ˜ ˜ orthogonal demixing matrix ;(2) W(τ )W equals forced-choice left or right key press reaction time task the cross-covariance matrix of the source components S upon two auditory stimuli in an oddball paradigm. The at time-lag τ; and (3) The independence assumption key press actions were performed with micro switches yields that the cross-covariance matrix of the source attached to the index fingers. During condition C0 sub- components S at time-lag τ is a diagonal matrix. TDSEP jects had to gaze at a fixation cross without any further thus computes ˜ as the matrix that jointly diagonalizes visual task. Condition C1 introduced an additional dis- a set of whitened cross-covariances Σ(τ). Here we use τ traction, as a video of a driving scene had to be watched = 1,..., 99. passively on a screen. Condition C2 introduced an addi- In the context of EEG signals, TDSEP finds k indepen- tional second task: subjects infrequently had to follow dent components contributing to the scalp EEG. They simple lane change instructions and control a steering are now characterized by their time course, a spatial wheel. By design, EEG recordings under condition C2 pattern given by the respective column of the mixing -1 were inevitably more prone to muscle and eye artifacts, matrix A := W , and a spatial filter given by the respec- whileC1 possiblystimulated eye movement artifacts, tive row of the demixing matrix W.The pattern con- but not muscle artifacts. However, all subjects had been tains the projection strengths of the respective instructed during all conditions to avoid producing component onto the scalp electrodes, whereas the filter artifacts. gives the projection strength of the scalp sensors onto Unmixing and data split the source component (see, e.g [24]). All resulting To avoid the artificial split of signal components due to source components were hand labeled into artifactual the high dimensionality of the data, the separation of and non-artifactual components by two experts who the EEG signals by an ICA method was preceded by a each labeled one half of the ICA components based on dimensionality reduction by Principal Component Ana- four plots per component, namely the time series, the lysis (PCA) from 121 EEG channels in the sensor space frequency spectrum and one scalp plot of the compo- into k = 30 PCA components. This choice of k was nent’s filter and one of its pattern. Not all components based on previous experience, but was probably not the were unambigious but instead contained a mixture of optimal choice. The TDSEP algorithm (Temporal Dec- neural and artifactual activity. Discarding all those com- orrelation source SEParation) [23] was used to trans- ponents which contain traces of artifacts would remove form the 30 PCA components into 30 independent too much of the relevant neural activity. Therefore, only source components. PCA and TDSEP were applied in a those mixed components were labeled as artifacts, that subject specific way, i.e. PCA and TDSEP matrices were revealed a relatively small amount of neural activity calculated seperately for each subject. compared to the strength of the artifact contained. TDSEP is a BSS algorithm to estimate a linear demix- For the training of the proposed automated classifica- ing tion method, 23 EEG recordings of 10 minutes duration were taken from the first experimental block only, lead- (1) WX = S ing to 690 labeled source components. Neural Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 5 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 components and artifact components were approx. Features derived from a component’s spectrum equally distributed (46% vs. 54%). Figure 1 shows typical 1. k , l, k and Fit Error describe the deviation of a 1 2 examples of two artifacts and one neural component. component’s spectrum from a prototypical 1/fre- quency curve and its shape. The parameters k , l, k The trained classifier was tested on 36 unseen EEG 1 2 >0 of the curve recordings from the third experimental blocks. Among these 1080 source components were 47% neuronal com- ponents and 53% artifact components. f → − k (2) Feature Extraction In order to provide substantial information to an auto- mated classification method, we construct an initial fea- ture set that contains 13 features from a component’s are determined by three points of the log spectrum: time series, 9 features from its spectrum and 16 from its (1) value at 2 Hz, (2) local minimum in the band 5- pattern. Based on this collection of 38 features a subset 13 Hz, (3) local minimum in the band 33-39 Hz. of the most discriminative features is determined in a The logarithm of k , l, k and of the mean squared 1 2 feature selection procedure. error of the approximation to the real spectrum are Features derived from a component’s time series used as features. 1. Variance of a component’s time series. It is not The spectrum of muscle artifacts, characterized by possible to determine the variances of the indepen- unusual high values in the 20-50 Hz range, are thus -1 dent components, as both S and A := W are approximated by a comparatively steep curve with unknown, and the solution is thus undetermined up high l and low k . to scaling. We estimate the impact one independent 2. 0-3 Hz, 4-7 Hz, 8-13 Hz, 14-30 Hz, 31-45 Hz, component s has on the original EEG by calculating the average log band power of the δ (0-3 Hz), θ (4-7 Var(std(A )· s )where A denotes the respective pat- i i i Hz), a (8-13 Hz), b (14-30 Hz) and g (31-45 Hz) tern. The idea here is to calculate the standard band. deviation of one independent component when its corresponding pattern has unit variance. Features derived from a component’s pattern 2. Maximum Amplitude 1. Range Within Pattern, logarithm of the differ- 3. Range of the signal amplitude ence between the minimal and the maximal activa- 4. Max First Derivative, approximated for the dis- tion in a pattern s(t )−s(t ) i+1 i crete signal s(t)in t by s (t ) ≈ i 2. Spatial Distance of Extrema, logarithm of the 5. Kurtosis Euclidean norm of the 2D-coordinates of the mini- 6. Shannon Entropy mal and maximal activation in a pattern 7. Deterministic Entropy, a computationally tract- 3. Spatial Mean Activation Left, Left Frontal, able measure related to the Kolmogorov complexity Frontal, Right Frontal, Right, Occipital, Central, of a signal [25] logarithm of the average activation in 7 groups of 8. Variance of Local Variance of time intervals of 1 electrodes as depicted in Figure 2 s and of 15 s duration (2 separate features) 4. 2DDFT. Pattern without a “smooth” activity dis- 9. Mean Local Variance of time intervals of 1 s tribution do not originate from an easily traceable duration, and of 15 s duration (2 separate features) psychological source and are thus artifacts or mixed 10. Mean Local Skewness, the mean absolute local components. The spatial frequency of a pattern can skewness of time intervals of 1 s and 15 s duration be described by means of a two-dimensional discrete (2 separate features) Fourier transformation. As a first step, the pattern is linearly interpolated to a quadratic 64x64 pattern The above 13 features were all logarithmized in a last matrix. The feature 2DDFT is the average logarith- step. With exception of the Variance feature all were mic band power of higher frequencies of the 1st and calculated after standardization of the time series to var- 4th quadrant (see Figure 3) of the 2D-Fourier spec- iance 1. These features describe outliers in terms of trum of the pattern matrix. unusual high amplitude values, as they are typically pre- 5. Laplace-Filter. Laplace-filtering leads a second sent in blinks and muscle artifacts. Furthermore, they way of finding spatially high frequent patterns, as are sensitive to non-stationarities and non-normal these have more defined edges. Similar to the higher order moments in the time series signal, as they 2DDFT-Feature, the pattern is linearly interpolated can be expected by muscle activity which typically is not to a quadratic 64 × 64 pattern matrix. Then, a 3 × 3 present equally strong over the full duration of 10 min. Laplace filter is applied. The feature is defined as the Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 6 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 Fpz Fpz Fp1 Fp2 Fp1 Fp2 AF7 AF8 AF7 AF8 AFp1 AFp2 AFp1 AFp2 F9 F10 F9 F10 AF3 AF4 AF3 AF4 F7 AF5F AF6F F8 F7 AF5F AF6F F8 AF1F AF2F AF1F AF2F F5 F6 F5 F6 F3 F4 F3 F4 FT9 FFC7 F1 Fz F2 FFC8 FT10 FT9 FFC7 F1 F2 FFC8 FT10 Fz FFC5 FFC6 FFC5 FFC6 FT7 FFC3 FFC1 FFC2 FFC4 FT8 FT7 FFC3 FFC4 FT8 FFC1 FFC2 FC5 FC6 FC5 FC6 FC3 FC1 FC2 FC4 FC3 FC4 FCz FC1 FCz FC2 CFC7 CFC8 CFC7 CFC8 CFC5 CFC3 CFC4 CFC6 CFC5 CFC6 CFC1 CFC2 CFC3 CFC1 CFC2 CFC4 T7 C5 C3 C1 Cz C2 C4 T8 C6 T7 C5 C3 C1 Cz C2 C4 T8 C6 CCP3 CCP1 CCP2 CCP4 CCP3 CCP1 CCP2 CCP4 CCP5 CCP6 CCP5 CCP6 CCP7 CCP8 CCP7 CCP8 CP1 CPz CP2 CP1 CPz CP2 CP3 CP4 CP3 CP4 CP5 CP6 CP5 CP6 TP7 PCP1 PCP2 TP8 PCP1 PCP2 PCP3 PCP4 TP7 PCP3 PCP4 TP8 PCP5 PCP6 PCP5 PCP6 TP9 PCP7 Pz PCP8 TP10 Pz P1 P2 TP9 PCP7 P1 P2 PCP8 TP10 P3 P4 P3 P4 P5 P6 PCP9 PPO1 PPO2 PCP10 PCP9 P5 P6 PCP10 PPO1 PPO2 PPO5 PPO6 P7 P8 P7 PPO5 PPO6 P8 POz PO1 PO2 PO1 POz PO2 PO3 PO4 PO3 PO4 P9 P10 P9 P10 PPO7 OPO1 OPO2 PPO8 OPO1 OPO2 PO7 PO8 PPO7 PPO8 PO7 PO8 O1 Oz O2 O1 O2 Oz OI1 OI2 OI1 OI2 I1 I2 I1 I2 Iz Iz Figure 4 The feature Border Activation. Electrode groups (left) and electrodes (right) used to determine the feature Border Activation. locations of the sources S. However, ICA patterns can be interpreted as EEG potentials for which a 3m physical model is given by a = Fz. Here, z Î ℝ are current moment vectors of unknown sources at m k×3m locations in the brain and F Î ℝ describes the mapping from sources to k sensors, which is deter- Figure 2 Scalp electrode sets. Meanactivationinthe7colored mined by the shape of the head and the conductiv- electrode groups are used as features. ities of brain, skull and skin tissues. We consider m = 2142 sources which are arranged in a 1 cm grid. logarithm of the Frobenius norm of the resulting Source estimation can only be done under additional matrix. constraints since k ≪ m. Commonly, the source dis- 6. Border Activation. This binary feature captures tribution with minimal l -norm (i.e., the “simplest” the spatial distribution at the borders of a pattern. It solution) is sought [26]. This leads to estimates is defined as 1 if either the global maximum of the 2 2 T T −1 T pattern is located at one of the outmost electrodes min ||Fz − a|| + λ||z|| =(F F + λ ) F a := J a λ (3) of the setup in Figure 4 (right), or if the local maxi- where Γ approximately equalizes the cost of dipoles mum of an electrode group in Figure 4 (left) is at different depths [27] and l defines a trade-off located at the outmost electrode of the group and if between the simplicity of the sources and the fidelity that local maximum deviates at least 2 standard of the model. deviations from the group average. Otherwise the feature is defined as -1. The idea behind this feature Since Eq. 3 models only cerebral sources, it is nat- is that a pattern with maximal activation at its bor- ural that noisy patterns and patterns originating out- der is unlikely to be generated by a source inside the side the brain can only be described by rather brain - it thus indicates an artifact. complicated sources, which are characterized by a 7. Current Density Norm of estimated source distri- large l -norm. For an example, see Figure 5. We pro- bution and strongest source’s position x, y, z.ICA pose to use f := log||z|| =log||J a|| as a feature itself does not provide information about the for discriminating physiological from noisy or arti- factual patterns. Here a ˜ := a/||a|| are normalized ICA patterns and l =100 waschosenfrom {0,1,10,100,1000} by cross-validation. To allow for a meaningful comparison of different f values over set- tings of varying numbers of electrodes, we pre-calcu- lated ΓJ on 115 electrodes and used only those rows that corresponded to the recorded electrodes. Note Pattern Pattern Matrix 2D S p ectrum that while this approach is simple, it may not be the Figure 3 The feature 2DDFT. The feature 2DDFT is the average optimal choice when the set of electrodes varies. logarithmic band power of higher frequencies of the 1st and 4th Assuming a pattern is generated by only one source, quadrant of the 2D-Fourier spectrum of the pattern matrix. we can estimate its 3D-coordinates x, y, z as the Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 7 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 task (defined by |E(w /||w ||)| >0.1) while the cross-vali- i i dation error deviates less than one standard error from the minimal cross-validation error. Having obtained a ranking of the features, the addi- tional information needed is how many of the best- ranked features are optimal for classification. With the goal in mind to find a good trade-off between feature size and error we proceed as follows: For every rank position, we compute the cross-validation error obtained by a classification based on the best-ranked features. Then the number of best ranked features is selected to be the minimum number of features yielding a cross- validation error which deviates less than one standard error from the minimal cross-validation error. Obviously, the number of features depends on the classification method. We compare a LPM, a non-linear Support Vector Machine (SVM) with Gaussian kernel [29] and a regularized Linear Discriminant Analysis (RLDA) [24], where we use a recently developed method to analytically calculate the optimal shrinkage parameter for regularization of LDA [30,31]. Since a nested cross-validation is computational expensive, the Figure 5 The feature Current Density Norm.(a)Twoexample hyperparameters of SVM and LPM are set by an outer patterns with high Current Density Norm f = log ||Γz||. (b) Two example patterns with low Current Density Norm. cross-validation, i.e. they are estimated on the whole training set which leads to a slight overfitting on the training data. location of maximal current density. Note that this is As a last step, the final classifier was trained on the only a very simple source localization method. full training data (690 examples) on the selected fea- tures, and tested on unseen test data (1080 examples). Feature Selection and Classification We conduct an embedded feature selection by using the Validation in an auditory ERP study weight vector of a Linear Programming Machine (LPM) To evaluate the artifact detection performance beyond [28]. Like all binary linear classifiers it finds a separating the training domain, data from 18 healthy subjects were d T hyperplane H : ℝ ∋ x ↦ sign(w · x + b) Î {-1, 1} char- used to test the proposed automated component classifi- acterized by a weight vector w and a bias term b.If the cation method in a completely different setup of an features are zero-mean and have same variance, their auditory ERP study. importance for the classification task can be ranked by Experimental Setup their respective absolute weights |w |. The LPM is A group of 18 subjects of 20 to 57 years of age (mean = known to produce a sparse weight vector w by solving 34.1,SD=11.4)underwent an EEGrecording of the following minimization problem: approx. 30 min duration using 64 Ag/AgCl electrodes of n approx. equidistant sensors. EEG was band-pass filtered min ||w|| + C ξ between 0.1-40 Hz. Note that this setup differs from the 1 i i=1 RT experiment, where EEG was recorded from 121 elec- (4) trodes and high-pass filtered at 2 Hz. s.t. y (w · x + b) ≥ 1 − ξ (i =1, ... , n) i i The subjects were situated in the center of a ring of ξ ≥ 0 (i =1, ... , n) six speakers (at ear height). During several short trials We thus apply a LPM to the training data in a 5 × 10 they listened to a rapid sequence (Stimulus Onset Asyn- cross-validation procedure with the goal to obtain a chrony = 175 ms) of six auditory stimuli of 40 ms dura- ranking of the features according to |E(w /||w||)|. tion. The six stimuli varied in pitch and noise. Each Beforehand, the LPM-hyperparameter C was set to C = stimulus type was presented from one speaker only, and 0.1 by a 5 × 10 cross-validation heuristic, such that each speaker emitted one stimulus type only such that LPM yielded good classification results while using a direction was a discriminant cue in addition to the sparse feature vector, i.e. we selected C with the mini- pitch/noise characteristics. Subjects had to count the mal number of features essential for the classification number of appearances of a rare target tone, that was Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 8 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 presented in a pseudo-random sequence together with 5 as a surrogate for the probability of being an artifact. frequent non-target tones (ratio 1:5). Retaining a smaller or larger number of sources corre- Unmixing and Classification sponds to an either very strict or soft policy for the A PCA reduced the dimensionality of the EEG channels removal of potential artifactual sources. We retained 6 to 30 PCA components. Then, the TDSEP algorithm to 30 source components of the most probable true was used to transform the 30 PCA components into 30 neural sources, and removed the others. Further analysis independent source components. The resulting 540 was performed on the remaining sources, i.e. CSP filters source components were hand labeled by two experts were determined on the remaining independent source into artifactual and non-artifactual source components. components and the log-variance of the spatially filtered One of the experts had participated in the rating of the signals were used to train an LDA. RT-study. Both experts rated all independent compo- Note that ICA artifact reduction methods usually nents. On average, the experts identified 28% neuronal reconstruct the EEG from the remaining neural sources. components and 72% artifactual components (expert 1: However, CSP solves an eigenvalue problem and 25% neuronal components, expert 2: 31% neuronal com- requires the covariance matrix of the data to have full ponents). The labeled data was used to test how the rank. Thus, CSP cannot be applied to the reconstructed artifact classifier generalizes to new data acquired in a EEG. different experimental setup by training the classifier The application to the feedback measurement in a solely on the training data from the RT experiment and manner that allows for real-time BCI applications is applying it to this unseen data set. straightforward: After un-mixing the original data according to the ICA filters determined on the calibra- Application to Motor Imagery BCI tion measurement, the previously determined 6 to 30 To investigate the possibility of removing relevant sources were selected for band-pass and CSP filtering neural activity, we incorporated our automatic ICA-clas- and log-variance determination in order to form the test sification step in a motor-imagery BCI system. In this data features. To estimate the influence of the artifact offine analysis we investigate how an ICA-artifact reduc- reduction step on BCI performance, we compared the tion step affects the classification performance of a classification performance with artifact reduction motor imagery BCI system based on the Common Spa- (depending on the number of selected sources) with the tial Patterns (CSP) method. For a detailed discussion of standard CSP procedure using no artifact reduction. CSP the reader is referred to [32]. Experimental Setup Results Eighty healthy BCI-novices performed first motor ima- In the following subsections, the results of the classifier gery with the left hand, right hand and both feet in a model selection and its additional validation on new calibration (i.e. without feedback) measurement. Every 8 data sets is presented. s one of three different visual cues (arrows pointing left, right, down) indicated to the subject which type of Model Selection: RT study motor imagery to perform. Three runs with 25 trials of The ranking of the features obtained by applying a LPM each motor condition were recorded. A classifier was to the training data set of the RT study is shown in trainedusing thepairofclassesthat providedbestdis- Table 1. Figure 7 shows the cross-validation errors for crimination: CSP filters were calculated on the band- SVM, RLDA and LPM plotted against the size of the pass filtered signals and the log-variance of the spatially feature sub-set used for classification. The shape of the filtered signals were used to train a LDA. In a feedback three curves reveals that at first, classification perfor- measurement subjects could control a 1D cursor appli- mance improves when adding features to the feature set. cation in three runs of 100 trials [33]. These features contain necessary but not redundant Motor Imagery BCI preceded by ICA-based artifact information. However, adding more than a certain num- reduction ber of features does not improve classification perfor- The steps conducted to incorporate the artifact reduc- mance-thesefeaturesonlycontain redundant tion are illustrated in Figure 6. The first step consists of information. Classification error slightly increases when a dimensionality reduction from about 90 EEG channels more features are used for the classification task, which in the sensor space into k = 30 PCA components. As in indicates that the classifier overfits on noisy and irrele- the previous experiments, TDSEP was used to transform vant features. the 30 PCA components into 30 independent source The fact that LPM performance is in the range of the components. Then, the component classifier trained on RLDA classifier indicates that the feature ranking was the RT experiment was applied. The components were suitable for our analysis (and not just for the LPM clas- ranked based on the classifiers output, which was used sifier). Given the ranking, the minimum number of $ Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 9 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 Figure 6 Artifact reduction step included in the standard CSP-procedure. The linear artifactreduction transformation of the original EEG into 6 - 30 signal components is calculated in the calibration phase. This transformation is applied to the feedback data. features yielding a cross-validation error which deviates that a high current density norm of a component indi- less than one standard error from the minimal cross- cates an artifactual component. Recall from the defini- validation erroris9forthe SVM and only 6for the tion of the Current Density Norm feature that these RLDA. The SVM classifier slightly outperforms the components are in fact difficult to explain by a promi- RLDA classifier on the training data, but since our goal nent source within the brain. Furthermore, components is to construct a simple linear classifier, we decided to with a high range within the pattern (i.e. outliers in the use the RLDA classifier with the 6 best-ranked features. pattern), a high local skewness (i.e. outliers in the time Notice that while SVM outperforms RLDA on the train- series), high l (i.e. a steep spectrum typical for muscle ing data, this effect might be due to overfitting and dis- artifacts) and low spectral power in the 8-13 Hz range appears on the test data, as is shown in the next section. (i.e. no prominent alpha peak) are rated as artifacts by The 6 best-ranked features are Current Density Norm, the classifier. Interestingly, a low FitError, i.e. a low Range Within Pattern, Mean Local Skewness 15 s, l, 8- error when approximating the spectrum by a 1/f curve, 13 Hz and FitError. They incorporate information from indicates an artifact for the classifier. This is due to the the temporal, spatial and frequency domain. fact that components which have no alpha peak in the spectrum are most probably artifacts. Notice that the Validation 1: RT study FitError feature in itself is not very informative, because Testing the trained classifier on unseen data from the a high FitError cannot distinguish between components RT study (1080 examples from experimental block 3) with a large alpha peak (which contain most probably leads to an mean-squared error (MSE) of 8.9% only, neural activity) and components with an unusual high which corresponds to a high agreement with the expert’s spectrum in higher frequency (which indicates muscle labeling. Interestingly, testing a trained SVM classifier activity). However, in combination with the other five (based on 9 selected features) leads to an error of 9.5%. features, the FitError feature carries additional informa- Thus, after feature selection, the RLDA classifier per- tion which improve classification performance. forms as good as a SVM classifier on unseen test data. It is interesting to take a closer look at the perfor- Let’s take a moment to interpret the obtained classi- mance of single features, which is also given in Table 2. fier: The weight vector w is giveninTable2. It shows The best one, Current Density Norm, leads to a MSE of Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 10 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 Table 1 Ranking of features obtained by LPM. 0.18 SVM Feature Weight RLDA Current Density Norm LPM 0.16 Range Within Pattern Mean Local Skewness 15 s 0.14 6 best-ranked features 8-13 Hz 0.12 FitError Border Activation 0.1 2DDFT Spatial Mean Activation Central 0.08 Max First Derivative Variance k 0.06 0 5 10 15 20 25 30 35 40 Spatial Mean Activation Left number of best-ranked features Spatial Mean Activation Left Frontal Figure 7 Cross-validation error for SVM, RLDA and LPM against Laplace-Filter the number of best-ranked features. A 10-fold cross-validation was repeated 5 times and standard errors are plotted. The SVM and Mean Local Variance 15 s LPM hyperparameters were selected by an outer cross validation. 14-30 Hz The number of 6 best-ranked features was determined for building 4-7 Hz the final classifier, as the estimated error of the RLDA starts to Mean Local Variance 1 s increase significantly for higher numbers of features. Spatial Distance of Extrema Spatial Mean Activation Occipital visual analysis of these cases reveals, that most of them were mixed components that contained both, artifacts Maximum Amplitude and brain activity. Out of the 21 components which were misclassified as neural activity only two were eye Spatial Mean Activation Right movements and none were blinks. In some rare cases, Kurtosis examples which had been mislabeled by the expert could be identified. Figure 8 shows an example of a mis- 0-3 Hz classified mixed component. Deterministic Entropy To quantify the classification performance on muscle Spatial Mean Activation Frontal artifacts, we asked one expert to review the 574 artifac- tual components of the test set for muscle activity. The Variance of Local Variance 1 s expert identified 388 components which contained mus- Range cle activity (which corresponds to 67.5% of the artifac- Spatial Mean Activation Right Frontal tual components and 17.2% of all the components). Out Variance of Local Variance 15 s of the 21 artifactual components which were Mean Local Skewness 1 s 31-45 Hz Shannon Entropy Table 2 Feature weight vector and test errors. The diameter of the black circles visualizes the absolute LPM weights |E(w /|| Feature Feature Test Error Test Error w||)| per feature after learning on the training data set. (The LPM- weight RT ERP hyperparameter C had been set to C = 0.1 based on the cross-validation performance.) Current Density Norm 0.342 0.141 0.488 Range Within Pattern 0.574 0.151 0.186 Mean Local Skewness 0.317 0.309 0.442 14.1% on the test data of the RT study. The combina- 15 s tion of the six features from all three domains improves l 0.569 0.177 0.144 the error substantially compared to even the best single 8-13 Hz -0.219 0.166 0.138 feature. This shows that features which are far from FitError -0.286 0.424 0.640 optimal in single classification have a positive contribu- Combined 0.089 0.147 tion in combination with other features. 6 T Feature weights w for each feature x of the classifier H : ℝ ∋ x a sign(w ·x Looking at the complete test set of 1080 components, i i + b) Î {-1, 1} with 1 ≙ Artifact and -1 ≙ Neuronal activity. Test error (MSE) for 75 of them were misclassified as artifacts and 21 compo- the 6 single features and for the combined classification for the RT nents were misclassified as neural sources. A detailed experiment and the ERP experiment. crossvalidation error Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 11 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 Validation 2: Auditory ERP study The classifier trained on RT data and applied to 540 components of the auditory ERP study leads to an aver- age MSE of 14.7% only for the classification of artifacts (expert 1: 15.7%, expert 2: 13.7%). On average over both experts, 18 of the 540 components were misclassified as artifacts and 61.5 components were misclassifed as neural sources (expert 1: 12 - 73, expert 2: 24 - 50). Table 2 also shows the classification results for every single feature and for the combined classification for the auditory ERP data. The classification performance of the three features Range Within Pattern, l and 8-13 Hz is comparable to those in the RT experiment. They general- ize very well over different experimental setups. However, the single feature classification performance for the remaining three features, Current Density Norm, Mean Local Skewness 15 s,and FitError, was close to chance Figure 8 Example for a misclassified component.Mixed level. This does not imply, however, that these features component that combines central alpha activity and slow Mayer Waves [35]. The human expert considered the mixed component a are unimportant for the classification tasks in the com- neural source, but the classifier labeld it as an artifact. bined feature set. To asses the relevance of each feature in thecombinedfeatureset,wetrained aRLDAonthe ERP data using the labels of expert 1 and report the fea- misclassified as neuronal components, only 12 contained ture weights of the weight vector - Current Density Norm muscle activity (57.1%). This indicates that muscle arti- 0.139; Range Within Pattern 0.355; Mean Local Skewness facts were handled equally well by the classifier as other 15 s 0.255; l 0.531; 8-13 Hz -0.710; FitError 0.059. We types of artifacts. found that while the feature weight of Current Density The performance of a system on the classification task Norm and Mean Local Skewness 15 s slightly decreased has to be judged in the light of the fact that inter-expert compared to the feature vector trained on the RT data, disagreements on EEG signals are often above 10% [34]. these were still far away from zero and thus carry infor- For our data, we asked one expert to re-label the 690 mation for the classification task. components of the training set, two years after the origi- nal labeling. The MSE between the new and the former Validation 3: Application to Motor Imagery BCI rating was 13.2%. Thus, the prediction performance of Figure 9 (left) plots the BCI classification error (1-AUC) our proposed classification method was comparable to against the number of remaining independent the ranking of an human expert. Figure 9 Influence of ICA-based artifact reduction in a motor imagery BCI tested with 80 subjects. Left: box plot of classifition errors (1- AUC) against the number of remaining independent sources compared with no artifact reduction. Right: scatter plot of classification errors (1- AUC) of each subject when removing 10 source components vs. using all 30 source components. Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 12 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 components, including one entry for the standard proce- We could show that the generalization over different dure without artifact reduction. Reducing the dimen- EEG studies is possible, which is in line with the find- sionality of the data to 30 dimensions by PCA does not ings of Mognon et al. [20] who demonstrated the gener- affect BCI performance. Moreover, consecutively remov- alization of an artifact classifier to a different laboratory ing components does not impair BCI performance at and to a different paradigm. Although their method is first, as these are artifactual components according to simple and efficient, it so far does not recognize muscle the classifier. Performance breaks down only when a artifacts. Compared to the classification results of Halder et al. strict removing policy is applied and less than about 12 [22], who reported 8% of error for muscle artifacts and sources (out of ~90 original channels) are retained, which have been ranked as neural sources by the classi- 1% error for eye artifacts, the classification error of our fier. The ranking of the classifier was confirmed by a solution is slightly higher. A major difference between visual analysis of the source components. Following the the two approaches is the way the training data was ranking of very probable artifacts to less probable arti- generated. Halder et al. reported, that subjects had spe- facts, the inspection resulted in clear artifactual compo- cifically been instructed to produce a number of artifacts nents to components that contained mixtures of neural under controlled conditions for the classifier training. It and artifactual activity. can bespeculatedthatsuchatraining set contains Figure 9 (right) shows a scatter plot of classification stronger artifacts and less erroneous labels. Nevertheless, errors (1-AUC) for each subject when removing 10 the results of Halder et al. were generated based on EEG source components vs. using all 30 source components. recordings of only 16 electrodes. Without adjustments, For this soft policy for removing artifactual components it can only be applied to EEG recordings with 16 elec- the variance between subjects is very small, especially trodes. In contrast, our method is applicable to different for subjects with good classification rates. EEG setups. However, we only tested the generalization ability over different EEG studies on electrode sets that Discussion covered the whole scalp with approx. equidistant sen- To summarize, we have constructed a subject-indepen- sors. Whether the classifier is applicable to deal with dent, fast, efficient linear component classification EEG data recorded with further reduced electrode sets method that automates the process of tedious hand- remains an open question that could be analyzed in the selection of e.g. artifactual independent components. future. Theproposedmethodisapplicableonlineand gener- To assess the danger of false positives introduced by alizes to new subjects without re-calibration. It delivers our artifact detection method,weevaluated theinflu- ence of a strong artifact reduction on the classification physiologically interpretable results, generalizes well over different experimental setups and is not limited performance of a standard motor imagery BCI task. An to a specific type of artifact. In particular, muscle arti- offline analysis of data acquired from 80 healthy subjects facts and eye artifacts (besides other types) are demonstrated that removing up to 60% of the sources recognized. (that were ranked according to their artifact classifier The proposed artifact classifier is based on six care- rating) did not impair the overall BCI classification per- fully constructed features that incorporate information formance. Note that we discarded the same number of from the spatial, the temporal and the spectral domain components per subject in order to analyse the effect of of the components and have been selected out of 38 fea- false positives. In a practical BCI-system, it would prob- tures by a thorough feature selection procedure. After ably be beneficial to apply a threshold on the propability its construction on data from a reaction time experi- of being an artifact per component instead. ment, the classifier’s performance was validated on two While the suitability of our approach to remove large different data sets: (1) on unseen data of a second con- artifactualsubspaces ofthedataisawelcomeresult, an dition of the original reaction time study - here the clas- open question remains that addresses the potential per- sifier achieved a classification error of 8.9%, while formance increase by careful artifact removal. Why disagreement between two ratings of experts was 13.2%. didn’t the removal of few artifactual sources improve (2)onunseendataofanauditory oddball ERPstudy - the average motor imagery BCI performance? It is here the classifier showed a classification error of 14.7% known that CSP is rather prone to outliers if the train- in comparison to 10.6% of disagreement between ing data set is small [36]. Strategies to overcome this experts. The classification error is remarkable low given problem include the use of regularization methods for CSP (such as invariant CSP [37] or robust CSP [38]), that the second study has been recorded with half the the explicit removal of outlier trials or channels and reg- number of electrodes, under a completely different para- ularization in the following classification step. As the digm, and contained a significantly higher proportion of artifactual components. standard evaluation procedures for motor imagery data Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 13 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 contained counter measures already (channel rejection classification method seems to be sufficient when the and trial rejection based on variance), and the number feature set is carefully constructed. of training data was considerably large, the overall influ- The classification difficulties of expert raters and of ence of artifacts on the motor imagery data set probably proposed automatic classification methods reflect the was small. Furthermore, we observed, that for subjects fundamental fact that any ICA-based artifact reduction with very good motor imagery classification rates, arti- method depends crucially on the quality of the source facts did not play any role at all. We conjecture that in separation into clear artifactual and neuronal source the other subjects, artifacts either obstructed the rele- components. A good source separation method avoids mixed components that contain both, neural and artifac- vant neural activity (cases where a slight improvement by artifact removal was obtained) or artifacts played tual activity as well as arbitrary splits of a single source some role in the control of the BCI system (cases where into several components. In the following, both type of artifact removal slightly reduced the performance). errors are briefly discussed. In addition to the construction of an efficient, sparse Blind source separation is a difficult problem by itself, and interpretable classifier, our feature-selection metho- and various approaches have been proposed to solve it dology leads to valuable insights into the question of (see, e.g. [39] for a review). In the context of EEG sig- which features are best suited for the discrimination of nals, the goal is to find a source separation that mini- artifactual and neuronal source components. However, mizes the amount of mixed components. The choice of it needs to be kept in mind that while the six identified TDSEP for the pre-processing of the EEG data was features were arguably an exceptionally suitable feature motivated by the ability of the algorithm to utilize tem- set, these features were probably not the overall optimal poral structure in the data. Although this is not a choice. Furthermore, the question remains if the unique feature of TDSEP, this approach seemed to be selected features generalize to other EEG data. Single suitable for the processing of EEG data, which is com- feature classification performance drops on the ERP posed of multidimensional time series signals with tem- data for three of the six features (Current Density Norm, poral dependencies. Moreover, research indicates that Mean Local Skewness 15 s and Fit Error). However, both methods based on second-order statistics might outper- Current Density Norm and Mean Local Skewness 15 s form methods based on higher-order statistics in the carry important information in the combined classifica- removal of ocular artifacts [10,22]. Although, as Fitzgib- tion (when used together with the other four features). bon stated, “the quality of the separation is highly Still, the non-redundant information carried by the Cur- dependent on the type of contamination, the degree of rent Density Norm feature drops substantially – apro- contamination, and the choice of BSS algorithm” [9], a blem that maybe causedbythe useofafixedmatrix thorough test of various ICA methods is out of the ΓJ which had been determined on the RT setup of 115 scope of this paper. electrodes. We found no obvious explanation for the The second kind of error, the arbitrary split of sources importance change for the Fit Error feature, however. into several components, can partially be compensated In any case, several insights can be gained concerning by combining the ICA with a preceding PCA step for the construction of a suitable feature set for the classifi- the dimensionality reduction. This procedure has the cation of artifactual components in general. First, the additional advantage of removing noise in the data. We spatial, the temporal and the spectral domain of the chose to project the original data into a 30 dimensional components contain non-redundant information. Sec- space. The value of 30 was based on rough experience ond, features that quantify aspects of the pattern’s activ- and on a quick visual inspection of the data, and was ity distribution, not its single values, are discriminative. probably not the optimal choice. An improvement of Features that were ranked high in our feature selection the quality of separation might be possible by optimiz- procedure were the range within the pattern, a feature ing the dimensionality reduction, but the effort was not based on the simplicity of a source separation, features undertaken here. Future work is needed to analyze the that analyzed the spatial frequency and a binary feature influence of dimensionality reduction on source which indicates if the maximal activation is on the bor- separation. der of the pattern. Third, features that model the shape To conclude, we hope that the source component of the power spectrum as a 1/f-curve as well as the classification method presented in this study delivers a absolute spectrum in the a range are discriminative. substantial contribution for the BCI community and the Fourth, features that quantify outliers in the time series EEG community in general, as a reliable and practical such as kurtosis, entropy, and mean local skewness, tool for the removal of artifacts. To support the com- munity, to encourage the reproduction of our results, or seem to be important but redundant. We analyzed 12 allow for re-labeling of data we provide the readily such features and only one obtained a high ranking in trained classifier, an implementation of the feature the feature selection. Last but not least, a linear Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 14 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 computer interfacing to mental state monitoring. Journal of neuroscience extraction routines together with example scripts, the methods 2008, 167:82-90. extracted features of the RT data, a visualization of 1770 3. Iwasaki M, Kellinghaus C, Alexopoulos AV, Burgess RC, Kumar AN, Han YH, components together with the expert labels used for the Lüders HO, Leigh RJ: Effects of eyelid closure, blinks, and eye movements on the electroencephalogram. Clinical Neurophysiology 2005, classifier training, and a visualization of components 116(4):878-885. misclassified by our method (see Additional File 1 - 4. Goncharova II, McFarland DJ, Vaughan TM, Wolpaw JR: EMG contamination MatlabCode; Additional File 2 - TrainComponents; of EEG: spectral and topographical characteristics. Clinical Neurophysiology 2003, 114:1580-1593. Additional File 3 - TestComponents; Additional File 4 - 5. Fatourechi M, Bashashati A, Ward RK, EBirch G: EMG and EOG artifacts in Misclassifications). brain computer interface systems: A survey. Clinical Neurophysiology 2007, 118:480-494. 6. Croft RJ, Barry RJ: Removal of ocular artifact from the EEG: a review. Additional material Clinical Neurophysiology 2000, 30:5-19. 7. Makeig S, Bell AJ, Jung TP, Sejnowski TJ: Independent Component Additional file 1: MatlabCode. A Matlab implementation of the feature Analysis of Electroencephalographic Data. Advances in neural information extraction routines together with example scripts, the readily trained processing systems 1996, 8:145-151. classifier and the extracted features for all the components of the RT 8. Jung TP, Makeig S, Humphries C, Lee TW, Mckeown MJ, Iragui V, data set. Sejnowski TJ: Removing electroencephalographic artifacts by blind source separation. Psychophysiology 2000, 37:163-178. Additional file 2: TrainComponents. Visualization of the 690 9. Fitzgibbon SP, Powers DMW, Pope KJ, Clark CR: Removal of EEG Noise and independent components in the training RT data, together with the Artifact Using Blind Source Separation. Clinical Neurophysiology 2007, expert’s labels. 24(3):232-243. Additional file 3: TestComponents. Visualization of the 1080 10. Romero S, Mañanas MA, Barbanoj MJ: A comparative study of automatic independent components in the RT test data, together with the expert’s techniques for ocular artifact reduction in spontaneous EEG signals labels. based on clinical target variables: A simulation case. Computers in Biology Additional file 4: Misclassifications. Visualization of the 75 + 21 and Medicine 2008, 38:348-360. misclassified components of the RT test data 11. Crespo-Garcia M, Atienza M, LCantero J: Muscle artifact remval from human sleep EEG by using independent component analysis. Ann Biomed Eng 2008, 36:467-475. 12. McMenamin BW, Shackman AJ, Maxwell JS, Bachhuber DRW, Koppenhaver AM, Greischar LL, Davidson RJ: Validation of ICA-based Acknowledgements myogenic artifact correction for scalp and source-localized EEG. The authors thank Martijn Schreuder for his help with recording and NeuroImage 2010, 49:2416-2432. preparing the auditory ERP data set, Claudia Sannelli for her help with 13. Barbati G, Porcaro C, Zappasodi F, Rossini PM, Tecchio F: Optimization of recording and preparing the RT EEG data set, the authors of [33] for an independent component analysis approach for artifact identifaction providing the motor imagery data set. and removal in magnetoencephalographic signals. Clinical Funding by the European Community under the PASCAL Network of Neurophysiology 2004, 115:1220-1232. Excellence (IST-2002-506778) and under the FP7 Programme (TOBI ICT-2007- 14. Delorme A, Makeig S, Sejnowski T: Automatic artifact rejection for EEG 224631 and ICT-216886), by the Bundesministerium für Bildung und Forschung data using high-order statistics and independent component analysis. (BMBF) (FKZ 01IBE01A, FKZ 16SV2231, 01GQ0850 and 01IB001A) and by the Proceedings of third international independent component analysis and blind Deutsche Forschungsgemeinschaft (DFG) (VitalBCI MU 987/3-1) is gratefully source decomposition conference, San Diego, CA 2001, 457-462. acknowledged. (This publication only reflects the authors’ views. Funding 15. Romero S, Mañanas M, Riba J, Giménez S, Clos S, Barbanoj M: Evaluation of agencies are not liable for any use that may be made of the information an automatic ocular filtering method for awake spontaneos EEG signals contained herein.) Last but not least, we would like to thank our reviewers based on Independt Component Analysis. 26th Annual International for their valuable comments. Conference of the Engineering in Medicine and Biology Society (EMBS) 2004, 925-928. Authors’ contributions 16. Shoker L, Sanei S, Chambers J: Artifact Removal from The EEG studies were performed in cooperation with colleagues mentioned Electroencephalograms using a hybrid BSS-SVM algorithm. IEEE Signal in the acknowledgments section. Authors MT and SH designed and carried Processing Letters 2005, 12(10):721-724. out the RT EEG study. MT designed and carried out the ERP study. IW and 17. James CJ, Gibson OJ: Temporally constrained ICA: an application to MT designed the feature extraction and feature selection algorithms and the artifact rejection in electromagnetic brain signal analysis. IEEE Trans artifact classification method. SH provided the current density norm feature Biomed Eng 2003, 50:1108-1116. method and implementation. All other implementations were carried out by 18. Joyce CA, Gorodnitsky IF, Kutas M: Automatic removal of eye movement IW. IW and MT analyzed and evaluated the overall methodology and wrote and blink artifacts from EEG data using blind component separation. the manuscript. All authors proof-read and approved the final manuscript. Psychophysiology 2004, 41:331-325. 19. Viola FC, Thorne J, Edmonds B, Schneider T, Eichele T, Debener S: Semi- Competing interests automatic identification of independent components representing EEG The authors declare that they have no competing interests. artifact. Clinical Neurophysiology 2009, 120:868-877. 20. Mognon A, Jovicich J, Bruzzone L, Buiatti M: ADJUST: An automatic EEG Received: 29 April 2011 Accepted: 2 August 2011 artifact detector based on the joint use of spatial and temporal features. Published: 2 August 2011 Psychophysiology 2010, 1-12. 21. LeVan P, Urrestarazu E, Gotman J: A system for automatic artifact removal References in ictal scalp EEG based on independent component analysis and 1. del R, Millán J, Rupp R, Mueller-Putz G, Murray-Smith R, Giugliemma C, Bayesian classification. Clinical Neurophysiology 2006, 117:912-927. Tangermann M, Vidaurre C, Cincotti F, Kübler A, Leeb R, Neuper C, 22. Halder S, Bensch M, Mellinger J, Bogdan M, Kübler A, Birbaumer N, Müller KR, Mattia D: Combining Brain-Computer Interfaces and Assistive Rosenstiel W: Online Artifact Removal for Brain-Computer Interfaces Technologies: State-of-the-Art and Challenges. Frontiers in Neuroprosthetics Using Support Vector Machines and Blind Source Separation. 2010, 4. Computational Intelligence and Neuroscience 2007, 7(3):1-10. 2. Müller KR, Tangermann M, Dornhege G, Krauledat M, Curio G, Blankertz B: 23. Ziehe A, Laskov P, Nolte G, Müller KR: A Fast Algorithm for Joint Machine Learning for real-time single-trial EEG-analysis: From brain- Diagonalization with Non-orthogonal Transformations and its Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 15 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 Application to Blind Source Separation. Journal of Machine Learning Research 2004, 5:801-818. 24. Blankertz B, Lemm S, Treder MS, Haufe S, Müller KR: Single-trial analysis and classification of ERP components - a tutorial. NeuroImage 2011. 25. Titchener MR: T-entropy of EEG/EOG Sensitive to Sleep State. International Symposium on Nonlinear Theory and Applications (NOLTA) 2006, 859-862. 26. Hämäläinen MS, Ilmoniemi R: Interpreting magnetic fields of the brain: minimum-norm estimates. Med Biol Eng Comput 1994, 32:35-42. 27. Haufe S, Nikulin VV, Ziehe A, Müller KR, Nolte G: Combining sparsity and rotational invariance in EEG/MEG source reconstruction. NeuroImage 2008, 42:726-738. 28. Bennett KP, Mangasarian OL: Robust linear programming discrimination of two linearly inseparable sets. Optimizaion Methods and Software 1994. 29. Chang CC, Lin CJ: LIBSVM. a library for support vector machines 2001. 30. Ledoit O, Wolf M: A well-conditioned estimator for large-dimensional covariance matrices. Journal of Multivariate Analysis 2004, 88(2):365-411. 31. Schäfer J, Strimmer K: A Shrinkage Approach to Large-Scale Covariance Matrix Estimation and Implications for Functional Genomics. Statistical Applications in Genetics and Molecular Biology 2005, 4(32). 32. Blankertz B, Tomioka R, Lemm S, Kawanabe M, Müller KR: Optimizing spatial filters for robust EEG single-trial analysis. IEEE Signal Proc Magazine 2008, 25:41-56. 33. Blankertz B, Sannelli C, Halder S, Hammer EM, Kübler A, Müller KR, Curio G, Dickhaus T: Neurophysiological predictor of SMR-based BCI performance. NeuroImage 2010, 51(4):1303-1309. 34. Klekowicz H, Malinow U, Niemcewi S, Wakarow A, Wolynczyk-Gmaj D, Piotrowski T, Durka P: Automatic analysis of sleep EEG. Frontiers in Neuroinformatics. Conference Abstract. Neuroinformatics 2008. 35. Lugaresi E, Coccagna G, Mantovani M, Lebrun R: Some periodic phenomena arising during drowsiness and sleeping in humans. Electroenceph clin Neurophysiol 1972, 32:701-705. 36. Krauledat M, Dornhege G, Blankertz B, Müller KR: Robustifying EEG data analysis by removing outliers. Chaos and Complexity Letters 2007, 2(3):259-274. 37. Blankertz B, Kawanabe M, Tomioka R, Hohlefeld F, Nikulin V, Müller KR: Invariant Common Spatial Patterns: Alleviating Nonstationarities in Brain-Computer Interfacing. In Advances in Neural Information Processing Systems 20. Edited by: Platt J, Koller D, Singer Y, Roweis S. Cambridge, MA: MIT Press; 2008:113-120. 38. Kawanabe M, Vidaurre C, Scholler S, Blankertz B, Müller KR: Robust Common Spatial Filters with a Maxmin Approach. EMBS-Conference 2009, 2470-2473. 39. Choi S, Cichocki A, Park HM, Lee SY: Blind Source Separation and Independent Component Analysis: A Review. Neural Information Processing - Letters and Reviews 2005, 6:1-57. doi:10.1186/1744-9081-7-30 Cite this article as: Winkler et al.: Automatic Classification of Artifactual ICA-Components for Artifact Removal in EEG Signals. Behavioral and Brain Functions 2011 7:30. Submit your next manuscript to BioMed Central and take full advantage of: • Convenient online submission • Thorough peer review • No space constraints or color figure charges • Immediate publication on acceptance • Inclusion in PubMed, CAS, Scopus and Google Scholar • Research which is freely available for redistribution Submit your manuscript at www.biomedcentral.com/submit http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Behavioral and Brain Functions Springer Journals

Automatic Classification of Artifactual ICA-Components for Artifact Removal in EEG Signals

Download PDF
15 pages

Loading next page...
 
/lp/springer-journals/automatic-classification-of-artifactual-ica-components-for-artifact-JPKtEYvlk9

References (95)

Publisher
Springer Journals
Copyright
Copyright © 2011 by Winkler et al; licensee BioMed Central Ltd.
Subject
Biomedicine; Neurosciences; Neurology; Behavioral Therapy; Psychiatry
eISSN
1744-9081
DOI
10.1186/1744-9081-7-30
pmid
21810266
Publisher site
See Article on Publisher Site

Abstract

Background: Artifacts contained in EEG recordings hamper both, the visual interpretation by experts as well as the algorithmic processing and analysis (e.g. for Brain-Computer Interfaces (BCI) or for Mental State Monitoring). While hand-optimized selection of source components derived from Independent Component Analysis (ICA) to clean EEG data is widespread, the field could greatly profit from automated solutions based on Machine Learning methods. Existing ICA-based removal strategies depend on explicit recordings of an individual’s artifacts or have not been shown to reliably identify muscle artifacts. Methods: We propose an automatic method for the classification of general artifactual source components. They are estimated by TDSEP, an ICA method that takes temporal correlations into account. The linear classifier is based on an optimized feature subset determined by a Linear Programming Machine (LPM). The subset is composed of features from the frequency-, the spatial- and temporal domain. A subject independent classifier was trained on 640 TDSEP components (reaction time (RT) study, n = 12) that were hand labeled by experts as artifactual or brain sources and tested on 1080 new components of RT data of the same study. Generalization was tested on new data from two studies (auditory Event Related Potential (ERP) paradigm, n = 18; motor imagery BCI paradigm, n = 80) that used data with different channel setups and from new subjects. Results: Based on six features only, the optimized linear classifier performed on level with the inter-expert disagreement (<10% Mean Squared Error (MSE)) on the RT data. On data of the auditory ERP study, the same pre- calculated classifier generalized well and achieved 15% MSE. On data of the motor imagery paradigm, we demonstrate that the discriminant information used for BCI is preserved when removing up to 60% of the most artifactual source components. Conclusions: We propose a universal and efficient classifier of ICA components for the subject independent removal of artifacts from EEG data. Based on linear methods, it is applicable for different electrode placements and supports the introspection of results. Trained on expert ratings of large data sets, it is not restricted to the detection of eye- and muscle artifacts. Its performance and generalization ability is demonstrated on data of different EEG studies. Background artifactual sources. Typical artifacts of the EEG are Signals of the electroencephalogram (EEG) can reflect caused either by the non-neural physiological activities the electrical background activity of the brain as well as of the subject or by external technical sources. Eye the activity which is specific for a cognitive task during blinks, eye movements, muscle activity in the vicinity of the head (e.g. face muscles, jaws, tongue, neck), heart an experiment. As the electrical field generated by neural activity is very small, it can only be recognized beat, pulse and Mayer waves are examples for physiolo- by EEG if large assemblies of neurons show a similar gical artifact sources, while swaying cables in the mag- behavior. Resulting neural EEG signals are in the range netic field of the earth, line humming, power supplies or of micro volts only and can easily be masked by transformers can be the cause of technical artifacts. Brain-Computer Interfaces (BCI) are based on the sin- * Correspondence: irene.winkler@tu-berlin.de gle trial classification of the ongoing EEG signal and can Machine Learning Laboratory, Berlin Institute of Technology, Franklinstr. 28/ improve the life quality of disabled individuals especially 29, 10587 Berlin, Germany © 2011 Winkler et al; licensee BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 2 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 in combination with other assitive technology [1]. The independent components, stationarity of the sources and exclusion of artifacts is of special interest for BCI appli- the mixture, and prior knowledge about the number of cations, as the intended or unconscious use of artifacts components), their application usually leads to a good for BCI control are usually not desirable when the BCI separation, with only a small number of hybrid compo- system is tested on healthy subjects. Furthermore, as nents that contain both, artifacts and neural signals averaging methods have to be avoided, these real-time [9-12]. Existing methods for artifact rejection can be sepa- systems BCIs rely on relatively clean EEG signals. The rated into hand-optimized, semi-automatic and fully same holds true for other Mental State Monitoring automatic approaches. Semi-automatic approaches applications, that monitor a subject’s mental state con- tinuouslyand onafinegranulartimeresolutionto require user interaction for ambiguous or outlier com- detect changes e.g. of wakefulness, responsiveness or ponents [13,14]. While fully automated methods were mental workload as early as possible [2]. proposed for the classification of eye artifacts [15,16], The two physiological artifacts most problematic for these methods do not easily generalize to non-eye arti- BCI applications are ocular (EOG) and muscle (EMG) facts or even require the additional recording of the artifacts. EOG activity is either caused by rolling of the EOG [17,18]. Viola et al. and Mognon et al. [19,20] both eyes or by eye blinks which occur approx. 20 times per developed an EEGLAB plug-in which finds artifactual minute [3]. Both result in a low-frequency activity most independent components. Both plug-ins have a fully prominent over the anterior head regions, with maximal automatic mode that has been shown to recognize and frequencies below 4 Hz. In contrast, EMG activity reject major artifacts like eye blinks, eye movements and (caused by chewing, swallowing, head or tongue move- heart beats, while the detection of muscular or more ments) is usually a high-frequency activity (>20 Hz) subtle artifacts has not been reported. The plug-in which ranges from rather small to very large amplitudes developed by Viola et al. relies on a user-defined tem- [4]. plate, while Mognon’s approach does not require user For an extensive review of artifact reduction techni- interaction. ques in the context of BCI-systems, the reader can refer Existing more flexible approaches for the general clas- to Fatourechi et al. [5]. Since the rejection of artifactual sification of different artifact types were reported for trials amounts to a considerable loss of data, a method EEG data of epileptic patients [21], where the authors that removes the artifacts while preserving the underly- report a Mean Squared Error (MSE) of approx. 20% for their system based on a Bayesian classifier. Halder et al. ing neural activity is needed. For example, linear filter- [22] report a classification error below 10% for their ing is a simple and effective method if artifactual and neural activity are located in non-overlapping frequency Support Vector Machine (SVM) based system for a bands. Unfortunately, artifacts and the brain signal of fixed number of electrodes if dedicated artifact record- interest do usually overlap. Nevertheless, ocular activity ings are available for the classifier training. But even if can be partially removed by regression-based methods, such optimized conditions are present, difficulties of which subtract a part of the activity measured at addi- separating muscle artifact components from neural com- tional electrooculogram (EOG) channels from the EEG ponents are common [22]. (see [6] for a review). Regression-based methods require The review of the existing literature did not reveal a the reliable recording of additional EOG channels and systematic screening of potentially discriminant features are limited by the fact that the EOG is contaminated by for the general task of artifact detection/removal. More- brain activity which is removed as well. Furthermore, over, most approaches restrict themselves to part of the they cannot eliminate non-eye activity. available information, e.g. rely on spatial patterns only If artifactual signal components and neural activity of [19], or spatial patterns and spectral features [22], or interest are not systematically co-activated due to a dis- spatial pattern and temporal features [20]. advantageous experimental design, methods of Blind Our proposed solution for a general artifact detection Source Separation (BSS) like Independent Component method is motivated by the needs of EEG practitioners. Analysis (ICA) are promising approaches for their First, it is desirable that a method efficiently and reliably separation [7,8]. A common approach is the transforma- detects all classes of artifacts,e.g.isnotrestrictedto tion of the EEG signals into a space of independent eye-, heart beat-, or muscle artifacts. Second, a practical source components, the hand-selection of non-artifac- method must be applicable post-hoc, i.e. without the need of dedicated artifact recordings at the time of the tual neural sources and the reconstruction of the EEG experiment. Third, it is difficult to convince EEG practi- without the artifactual components (for an example of tioners to use a method of artifact rejection if it is a independent source components, see Figure 1). While assumptions for the application of ICA methods are black box and refuses introspection. As the goal must only approximately met in practice (linear mixture of be to develop a method, that delivers interpretable and Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 3 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 (a) (b) (c) sec Hz Filter P atte rn Figure 1 Three example independent source components. Time series (first column), spectrum (second column), filter (third column) and pattern (fourth column) of three components. The first row (a) shows an alpha generator in the occipital lobe. The second row (b) shows a rare muscle artifact component with an increased spectrum in higher frequencies. The third line (c) shows an eye artifact component that appears regularly, has an increases spectrum in lower frequencies and a typical front-back distribution in the pattern. easy to understand results, we decided for a linear clas- signal pre-processing methods (including a temporal sification method. Luckily, linear methods have proven a variant of ICA) are introduced, we describe 38 features high performance for a number of classification tasks in that are candidates for the artifact discrimination task. the field of EEG-based BCI systems. However, to be able BasedonlabelsprovidedbyEEGexperts, athorough to estimate the performance loss compared to a poten- feature selection procedure is described, that is used to tially better, but difficult to interpret, non-linear classifi- condense the 38 features to a small subset. Furthermore, cation method, the results of a Gaussian SVM are classification methods are introduced. The methods sec- reported in parallel. tion ends with a description of two other EEG para- We decided to use a sparse approach (sparsity in the digms (auditory Event Related Potential (ERP) and features) although it is a mixed blessing. It leads to a motor imagery for BCI), that will be used to validate the generalization approach of the proposed artifact classi- trade-off between efficiency and interpretability, as redundant but slightly less discriminative features are fier. In the results section, the outcome of the feature removed with high probability from the overall set of selection procedure is given, together with the artifact features. This has to be kept in mind during the analysis classification performance on unseen data of the RT of results. To reach the goal of a sparse method that paradigm, data of a unseen auditory ERP paradigm. delivers physiologically interpretable results, we decided Finally the method is applied in the context of a motor to incorporate a thorough feature selection procedure in imagery BCI setup, before the paper closes with a combination with a linear classification method that is discussion. based on features of all three available information domains of EEG data: the spatial domain (e.g. patterns Methods of independent components), the frequency domain and In the following subsections, we will describe how the the temporal domain. proposed new artifact classification method is set up. The paper is organized as follows: In the methods sec- Then we will introduce two further studies that are uti- tion, a reaction time (RT) paradigm is introduced, as lized to test the classifier’s generalizability. data from this study forms the basis for the construction Participants of the studies described below provided of the proposed artifact detection method. After the verbal and written informedconsent andwere freeto dB Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 4 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 stop their participation at any time. All collected data of a given multivariate time series X =(x ,..., x ) into 1 k was anonymized before any subsequent analysis or pre- unknown, assumed mutually independent source com- sentation took place. ponents S =(s ,..., s ) . Note that both the demixing W 1 k and the source components S are unknown, and that Classifier Construction using a RT study BSS algorithms differ in the definition of independence The artifact classifier is set up based on labeled indepen- between components. While ICA algorithms exploit dent components gained from a reaction time (RT) higher order statistics, TDSEP relies on second-order study. statistics by taking the temporal structure of the time Experimental Setup series into account. TDSEP amounts to finding a demix- Data from 12 healthy right-handed male subjects were ing W which leads to minimal cross-covariances over used to train and to test the proposed automated com- several time-lags between all pairs of components of S. ponent classification method. Every subject participated For a mathematical discussion, let in one EEG recording session of approx. 5 hours dura- be the cross-covariance (τ):= E(X (t)X (t − τ )) tion. EEG was recorded from 121 approx. equidistant matrix of the whitened data X at time-lag τ,where the sensors and high pass filtered at 2 Hz. During this ses- whitening transformation linearly decorrellates and sion, 4 repeated blocks of 3 different conditions (C0, C1, scales the data such that Σ(0) = I. Consider now that (1) C2) were performed. Each block lasted approx. 45 min- Whitening reduces the BSS problem to finding an utes. During all three conditions, subjects performed a ˜ ˜ orthogonal demixing matrix ;(2) W(τ )W equals forced-choice left or right key press reaction time task the cross-covariance matrix of the source components S upon two auditory stimuli in an oddball paradigm. The at time-lag τ; and (3) The independence assumption key press actions were performed with micro switches yields that the cross-covariance matrix of the source attached to the index fingers. During condition C0 sub- components S at time-lag τ is a diagonal matrix. TDSEP jects had to gaze at a fixation cross without any further thus computes ˜ as the matrix that jointly diagonalizes visual task. Condition C1 introduced an additional dis- a set of whitened cross-covariances Σ(τ). Here we use τ traction, as a video of a driving scene had to be watched = 1,..., 99. passively on a screen. Condition C2 introduced an addi- In the context of EEG signals, TDSEP finds k indepen- tional second task: subjects infrequently had to follow dent components contributing to the scalp EEG. They simple lane change instructions and control a steering are now characterized by their time course, a spatial wheel. By design, EEG recordings under condition C2 pattern given by the respective column of the mixing -1 were inevitably more prone to muscle and eye artifacts, matrix A := W , and a spatial filter given by the respec- whileC1 possiblystimulated eye movement artifacts, tive row of the demixing matrix W.The pattern con- but not muscle artifacts. However, all subjects had been tains the projection strengths of the respective instructed during all conditions to avoid producing component onto the scalp electrodes, whereas the filter artifacts. gives the projection strength of the scalp sensors onto Unmixing and data split the source component (see, e.g [24]). All resulting To avoid the artificial split of signal components due to source components were hand labeled into artifactual the high dimensionality of the data, the separation of and non-artifactual components by two experts who the EEG signals by an ICA method was preceded by a each labeled one half of the ICA components based on dimensionality reduction by Principal Component Ana- four plots per component, namely the time series, the lysis (PCA) from 121 EEG channels in the sensor space frequency spectrum and one scalp plot of the compo- into k = 30 PCA components. This choice of k was nent’s filter and one of its pattern. Not all components based on previous experience, but was probably not the were unambigious but instead contained a mixture of optimal choice. The TDSEP algorithm (Temporal Dec- neural and artifactual activity. Discarding all those com- orrelation source SEParation) [23] was used to trans- ponents which contain traces of artifacts would remove form the 30 PCA components into 30 independent too much of the relevant neural activity. Therefore, only source components. PCA and TDSEP were applied in a those mixed components were labeled as artifacts, that subject specific way, i.e. PCA and TDSEP matrices were revealed a relatively small amount of neural activity calculated seperately for each subject. compared to the strength of the artifact contained. TDSEP is a BSS algorithm to estimate a linear demix- For the training of the proposed automated classifica- ing tion method, 23 EEG recordings of 10 minutes duration were taken from the first experimental block only, lead- (1) WX = S ing to 690 labeled source components. Neural Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 5 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 components and artifact components were approx. Features derived from a component’s spectrum equally distributed (46% vs. 54%). Figure 1 shows typical 1. k , l, k and Fit Error describe the deviation of a 1 2 examples of two artifacts and one neural component. component’s spectrum from a prototypical 1/fre- quency curve and its shape. The parameters k , l, k The trained classifier was tested on 36 unseen EEG 1 2 >0 of the curve recordings from the third experimental blocks. Among these 1080 source components were 47% neuronal com- ponents and 53% artifact components. f → − k (2) Feature Extraction In order to provide substantial information to an auto- mated classification method, we construct an initial fea- ture set that contains 13 features from a component’s are determined by three points of the log spectrum: time series, 9 features from its spectrum and 16 from its (1) value at 2 Hz, (2) local minimum in the band 5- pattern. Based on this collection of 38 features a subset 13 Hz, (3) local minimum in the band 33-39 Hz. of the most discriminative features is determined in a The logarithm of k , l, k and of the mean squared 1 2 feature selection procedure. error of the approximation to the real spectrum are Features derived from a component’s time series used as features. 1. Variance of a component’s time series. It is not The spectrum of muscle artifacts, characterized by possible to determine the variances of the indepen- unusual high values in the 20-50 Hz range, are thus -1 dent components, as both S and A := W are approximated by a comparatively steep curve with unknown, and the solution is thus undetermined up high l and low k . to scaling. We estimate the impact one independent 2. 0-3 Hz, 4-7 Hz, 8-13 Hz, 14-30 Hz, 31-45 Hz, component s has on the original EEG by calculating the average log band power of the δ (0-3 Hz), θ (4-7 Var(std(A )· s )where A denotes the respective pat- i i i Hz), a (8-13 Hz), b (14-30 Hz) and g (31-45 Hz) tern. The idea here is to calculate the standard band. deviation of one independent component when its corresponding pattern has unit variance. Features derived from a component’s pattern 2. Maximum Amplitude 1. Range Within Pattern, logarithm of the differ- 3. Range of the signal amplitude ence between the minimal and the maximal activa- 4. Max First Derivative, approximated for the dis- tion in a pattern s(t )−s(t ) i+1 i crete signal s(t)in t by s (t ) ≈ i 2. Spatial Distance of Extrema, logarithm of the 5. Kurtosis Euclidean norm of the 2D-coordinates of the mini- 6. Shannon Entropy mal and maximal activation in a pattern 7. Deterministic Entropy, a computationally tract- 3. Spatial Mean Activation Left, Left Frontal, able measure related to the Kolmogorov complexity Frontal, Right Frontal, Right, Occipital, Central, of a signal [25] logarithm of the average activation in 7 groups of 8. Variance of Local Variance of time intervals of 1 electrodes as depicted in Figure 2 s and of 15 s duration (2 separate features) 4. 2DDFT. Pattern without a “smooth” activity dis- 9. Mean Local Variance of time intervals of 1 s tribution do not originate from an easily traceable duration, and of 15 s duration (2 separate features) psychological source and are thus artifacts or mixed 10. Mean Local Skewness, the mean absolute local components. The spatial frequency of a pattern can skewness of time intervals of 1 s and 15 s duration be described by means of a two-dimensional discrete (2 separate features) Fourier transformation. As a first step, the pattern is linearly interpolated to a quadratic 64x64 pattern The above 13 features were all logarithmized in a last matrix. The feature 2DDFT is the average logarith- step. With exception of the Variance feature all were mic band power of higher frequencies of the 1st and calculated after standardization of the time series to var- 4th quadrant (see Figure 3) of the 2D-Fourier spec- iance 1. These features describe outliers in terms of trum of the pattern matrix. unusual high amplitude values, as they are typically pre- 5. Laplace-Filter. Laplace-filtering leads a second sent in blinks and muscle artifacts. Furthermore, they way of finding spatially high frequent patterns, as are sensitive to non-stationarities and non-normal these have more defined edges. Similar to the higher order moments in the time series signal, as they 2DDFT-Feature, the pattern is linearly interpolated can be expected by muscle activity which typically is not to a quadratic 64 × 64 pattern matrix. Then, a 3 × 3 present equally strong over the full duration of 10 min. Laplace filter is applied. The feature is defined as the Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 6 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 Fpz Fpz Fp1 Fp2 Fp1 Fp2 AF7 AF8 AF7 AF8 AFp1 AFp2 AFp1 AFp2 F9 F10 F9 F10 AF3 AF4 AF3 AF4 F7 AF5F AF6F F8 F7 AF5F AF6F F8 AF1F AF2F AF1F AF2F F5 F6 F5 F6 F3 F4 F3 F4 FT9 FFC7 F1 Fz F2 FFC8 FT10 FT9 FFC7 F1 F2 FFC8 FT10 Fz FFC5 FFC6 FFC5 FFC6 FT7 FFC3 FFC1 FFC2 FFC4 FT8 FT7 FFC3 FFC4 FT8 FFC1 FFC2 FC5 FC6 FC5 FC6 FC3 FC1 FC2 FC4 FC3 FC4 FCz FC1 FCz FC2 CFC7 CFC8 CFC7 CFC8 CFC5 CFC3 CFC4 CFC6 CFC5 CFC6 CFC1 CFC2 CFC3 CFC1 CFC2 CFC4 T7 C5 C3 C1 Cz C2 C4 T8 C6 T7 C5 C3 C1 Cz C2 C4 T8 C6 CCP3 CCP1 CCP2 CCP4 CCP3 CCP1 CCP2 CCP4 CCP5 CCP6 CCP5 CCP6 CCP7 CCP8 CCP7 CCP8 CP1 CPz CP2 CP1 CPz CP2 CP3 CP4 CP3 CP4 CP5 CP6 CP5 CP6 TP7 PCP1 PCP2 TP8 PCP1 PCP2 PCP3 PCP4 TP7 PCP3 PCP4 TP8 PCP5 PCP6 PCP5 PCP6 TP9 PCP7 Pz PCP8 TP10 Pz P1 P2 TP9 PCP7 P1 P2 PCP8 TP10 P3 P4 P3 P4 P5 P6 PCP9 PPO1 PPO2 PCP10 PCP9 P5 P6 PCP10 PPO1 PPO2 PPO5 PPO6 P7 P8 P7 PPO5 PPO6 P8 POz PO1 PO2 PO1 POz PO2 PO3 PO4 PO3 PO4 P9 P10 P9 P10 PPO7 OPO1 OPO2 PPO8 OPO1 OPO2 PO7 PO8 PPO7 PPO8 PO7 PO8 O1 Oz O2 O1 O2 Oz OI1 OI2 OI1 OI2 I1 I2 I1 I2 Iz Iz Figure 4 The feature Border Activation. Electrode groups (left) and electrodes (right) used to determine the feature Border Activation. locations of the sources S. However, ICA patterns can be interpreted as EEG potentials for which a 3m physical model is given by a = Fz. Here, z Î ℝ are current moment vectors of unknown sources at m k×3m locations in the brain and F Î ℝ describes the mapping from sources to k sensors, which is deter- Figure 2 Scalp electrode sets. Meanactivationinthe7colored mined by the shape of the head and the conductiv- electrode groups are used as features. ities of brain, skull and skin tissues. We consider m = 2142 sources which are arranged in a 1 cm grid. logarithm of the Frobenius norm of the resulting Source estimation can only be done under additional matrix. constraints since k ≪ m. Commonly, the source dis- 6. Border Activation. This binary feature captures tribution with minimal l -norm (i.e., the “simplest” the spatial distribution at the borders of a pattern. It solution) is sought [26]. This leads to estimates is defined as 1 if either the global maximum of the 2 2 T T −1 T pattern is located at one of the outmost electrodes min ||Fz − a|| + λ||z|| =(F F + λ ) F a := J a λ (3) of the setup in Figure 4 (right), or if the local maxi- where Γ approximately equalizes the cost of dipoles mum of an electrode group in Figure 4 (left) is at different depths [27] and l defines a trade-off located at the outmost electrode of the group and if between the simplicity of the sources and the fidelity that local maximum deviates at least 2 standard of the model. deviations from the group average. Otherwise the feature is defined as -1. The idea behind this feature Since Eq. 3 models only cerebral sources, it is nat- is that a pattern with maximal activation at its bor- ural that noisy patterns and patterns originating out- der is unlikely to be generated by a source inside the side the brain can only be described by rather brain - it thus indicates an artifact. complicated sources, which are characterized by a 7. Current Density Norm of estimated source distri- large l -norm. For an example, see Figure 5. We pro- bution and strongest source’s position x, y, z.ICA pose to use f := log||z|| =log||J a|| as a feature itself does not provide information about the for discriminating physiological from noisy or arti- factual patterns. Here a ˜ := a/||a|| are normalized ICA patterns and l =100 waschosenfrom {0,1,10,100,1000} by cross-validation. To allow for a meaningful comparison of different f values over set- tings of varying numbers of electrodes, we pre-calcu- lated ΓJ on 115 electrodes and used only those rows that corresponded to the recorded electrodes. Note Pattern Pattern Matrix 2D S p ectrum that while this approach is simple, it may not be the Figure 3 The feature 2DDFT. The feature 2DDFT is the average optimal choice when the set of electrodes varies. logarithmic band power of higher frequencies of the 1st and 4th Assuming a pattern is generated by only one source, quadrant of the 2D-Fourier spectrum of the pattern matrix. we can estimate its 3D-coordinates x, y, z as the Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 7 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 task (defined by |E(w /||w ||)| >0.1) while the cross-vali- i i dation error deviates less than one standard error from the minimal cross-validation error. Having obtained a ranking of the features, the addi- tional information needed is how many of the best- ranked features are optimal for classification. With the goal in mind to find a good trade-off between feature size and error we proceed as follows: For every rank position, we compute the cross-validation error obtained by a classification based on the best-ranked features. Then the number of best ranked features is selected to be the minimum number of features yielding a cross- validation error which deviates less than one standard error from the minimal cross-validation error. Obviously, the number of features depends on the classification method. We compare a LPM, a non-linear Support Vector Machine (SVM) with Gaussian kernel [29] and a regularized Linear Discriminant Analysis (RLDA) [24], where we use a recently developed method to analytically calculate the optimal shrinkage parameter for regularization of LDA [30,31]. Since a nested cross-validation is computational expensive, the Figure 5 The feature Current Density Norm.(a)Twoexample hyperparameters of SVM and LPM are set by an outer patterns with high Current Density Norm f = log ||Γz||. (b) Two example patterns with low Current Density Norm. cross-validation, i.e. they are estimated on the whole training set which leads to a slight overfitting on the training data. location of maximal current density. Note that this is As a last step, the final classifier was trained on the only a very simple source localization method. full training data (690 examples) on the selected fea- tures, and tested on unseen test data (1080 examples). Feature Selection and Classification We conduct an embedded feature selection by using the Validation in an auditory ERP study weight vector of a Linear Programming Machine (LPM) To evaluate the artifact detection performance beyond [28]. Like all binary linear classifiers it finds a separating the training domain, data from 18 healthy subjects were d T hyperplane H : ℝ ∋ x ↦ sign(w · x + b) Î {-1, 1} char- used to test the proposed automated component classifi- acterized by a weight vector w and a bias term b.If the cation method in a completely different setup of an features are zero-mean and have same variance, their auditory ERP study. importance for the classification task can be ranked by Experimental Setup their respective absolute weights |w |. The LPM is A group of 18 subjects of 20 to 57 years of age (mean = known to produce a sparse weight vector w by solving 34.1,SD=11.4)underwent an EEGrecording of the following minimization problem: approx. 30 min duration using 64 Ag/AgCl electrodes of n approx. equidistant sensors. EEG was band-pass filtered min ||w|| + C ξ between 0.1-40 Hz. Note that this setup differs from the 1 i i=1 RT experiment, where EEG was recorded from 121 elec- (4) trodes and high-pass filtered at 2 Hz. s.t. y (w · x + b) ≥ 1 − ξ (i =1, ... , n) i i The subjects were situated in the center of a ring of ξ ≥ 0 (i =1, ... , n) six speakers (at ear height). During several short trials We thus apply a LPM to the training data in a 5 × 10 they listened to a rapid sequence (Stimulus Onset Asyn- cross-validation procedure with the goal to obtain a chrony = 175 ms) of six auditory stimuli of 40 ms dura- ranking of the features according to |E(w /||w||)|. tion. The six stimuli varied in pitch and noise. Each Beforehand, the LPM-hyperparameter C was set to C = stimulus type was presented from one speaker only, and 0.1 by a 5 × 10 cross-validation heuristic, such that each speaker emitted one stimulus type only such that LPM yielded good classification results while using a direction was a discriminant cue in addition to the sparse feature vector, i.e. we selected C with the mini- pitch/noise characteristics. Subjects had to count the mal number of features essential for the classification number of appearances of a rare target tone, that was Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 8 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 presented in a pseudo-random sequence together with 5 as a surrogate for the probability of being an artifact. frequent non-target tones (ratio 1:5). Retaining a smaller or larger number of sources corre- Unmixing and Classification sponds to an either very strict or soft policy for the A PCA reduced the dimensionality of the EEG channels removal of potential artifactual sources. We retained 6 to 30 PCA components. Then, the TDSEP algorithm to 30 source components of the most probable true was used to transform the 30 PCA components into 30 neural sources, and removed the others. Further analysis independent source components. The resulting 540 was performed on the remaining sources, i.e. CSP filters source components were hand labeled by two experts were determined on the remaining independent source into artifactual and non-artifactual source components. components and the log-variance of the spatially filtered One of the experts had participated in the rating of the signals were used to train an LDA. RT-study. Both experts rated all independent compo- Note that ICA artifact reduction methods usually nents. On average, the experts identified 28% neuronal reconstruct the EEG from the remaining neural sources. components and 72% artifactual components (expert 1: However, CSP solves an eigenvalue problem and 25% neuronal components, expert 2: 31% neuronal com- requires the covariance matrix of the data to have full ponents). The labeled data was used to test how the rank. Thus, CSP cannot be applied to the reconstructed artifact classifier generalizes to new data acquired in a EEG. different experimental setup by training the classifier The application to the feedback measurement in a solely on the training data from the RT experiment and manner that allows for real-time BCI applications is applying it to this unseen data set. straightforward: After un-mixing the original data according to the ICA filters determined on the calibra- Application to Motor Imagery BCI tion measurement, the previously determined 6 to 30 To investigate the possibility of removing relevant sources were selected for band-pass and CSP filtering neural activity, we incorporated our automatic ICA-clas- and log-variance determination in order to form the test sification step in a motor-imagery BCI system. In this data features. To estimate the influence of the artifact offine analysis we investigate how an ICA-artifact reduc- reduction step on BCI performance, we compared the tion step affects the classification performance of a classification performance with artifact reduction motor imagery BCI system based on the Common Spa- (depending on the number of selected sources) with the tial Patterns (CSP) method. For a detailed discussion of standard CSP procedure using no artifact reduction. CSP the reader is referred to [32]. Experimental Setup Results Eighty healthy BCI-novices performed first motor ima- In the following subsections, the results of the classifier gery with the left hand, right hand and both feet in a model selection and its additional validation on new calibration (i.e. without feedback) measurement. Every 8 data sets is presented. s one of three different visual cues (arrows pointing left, right, down) indicated to the subject which type of Model Selection: RT study motor imagery to perform. Three runs with 25 trials of The ranking of the features obtained by applying a LPM each motor condition were recorded. A classifier was to the training data set of the RT study is shown in trainedusing thepairofclassesthat providedbestdis- Table 1. Figure 7 shows the cross-validation errors for crimination: CSP filters were calculated on the band- SVM, RLDA and LPM plotted against the size of the pass filtered signals and the log-variance of the spatially feature sub-set used for classification. The shape of the filtered signals were used to train a LDA. In a feedback three curves reveals that at first, classification perfor- measurement subjects could control a 1D cursor appli- mance improves when adding features to the feature set. cation in three runs of 100 trials [33]. These features contain necessary but not redundant Motor Imagery BCI preceded by ICA-based artifact information. However, adding more than a certain num- reduction ber of features does not improve classification perfor- The steps conducted to incorporate the artifact reduc- mance-thesefeaturesonlycontain redundant tion are illustrated in Figure 6. The first step consists of information. Classification error slightly increases when a dimensionality reduction from about 90 EEG channels more features are used for the classification task, which in the sensor space into k = 30 PCA components. As in indicates that the classifier overfits on noisy and irrele- the previous experiments, TDSEP was used to transform vant features. the 30 PCA components into 30 independent source The fact that LPM performance is in the range of the components. Then, the component classifier trained on RLDA classifier indicates that the feature ranking was the RT experiment was applied. The components were suitable for our analysis (and not just for the LPM clas- ranked based on the classifiers output, which was used sifier). Given the ranking, the minimum number of $ Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 9 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 Figure 6 Artifact reduction step included in the standard CSP-procedure. The linear artifactreduction transformation of the original EEG into 6 - 30 signal components is calculated in the calibration phase. This transformation is applied to the feedback data. features yielding a cross-validation error which deviates that a high current density norm of a component indi- less than one standard error from the minimal cross- cates an artifactual component. Recall from the defini- validation erroris9forthe SVM and only 6for the tion of the Current Density Norm feature that these RLDA. The SVM classifier slightly outperforms the components are in fact difficult to explain by a promi- RLDA classifier on the training data, but since our goal nent source within the brain. Furthermore, components is to construct a simple linear classifier, we decided to with a high range within the pattern (i.e. outliers in the use the RLDA classifier with the 6 best-ranked features. pattern), a high local skewness (i.e. outliers in the time Notice that while SVM outperforms RLDA on the train- series), high l (i.e. a steep spectrum typical for muscle ing data, this effect might be due to overfitting and dis- artifacts) and low spectral power in the 8-13 Hz range appears on the test data, as is shown in the next section. (i.e. no prominent alpha peak) are rated as artifacts by The 6 best-ranked features are Current Density Norm, the classifier. Interestingly, a low FitError, i.e. a low Range Within Pattern, Mean Local Skewness 15 s, l, 8- error when approximating the spectrum by a 1/f curve, 13 Hz and FitError. They incorporate information from indicates an artifact for the classifier. This is due to the the temporal, spatial and frequency domain. fact that components which have no alpha peak in the spectrum are most probably artifacts. Notice that the Validation 1: RT study FitError feature in itself is not very informative, because Testing the trained classifier on unseen data from the a high FitError cannot distinguish between components RT study (1080 examples from experimental block 3) with a large alpha peak (which contain most probably leads to an mean-squared error (MSE) of 8.9% only, neural activity) and components with an unusual high which corresponds to a high agreement with the expert’s spectrum in higher frequency (which indicates muscle labeling. Interestingly, testing a trained SVM classifier activity). However, in combination with the other five (based on 9 selected features) leads to an error of 9.5%. features, the FitError feature carries additional informa- Thus, after feature selection, the RLDA classifier per- tion which improve classification performance. forms as good as a SVM classifier on unseen test data. It is interesting to take a closer look at the perfor- Let’s take a moment to interpret the obtained classi- mance of single features, which is also given in Table 2. fier: The weight vector w is giveninTable2. It shows The best one, Current Density Norm, leads to a MSE of Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 10 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 Table 1 Ranking of features obtained by LPM. 0.18 SVM Feature Weight RLDA Current Density Norm LPM 0.16 Range Within Pattern Mean Local Skewness 15 s 0.14 6 best-ranked features 8-13 Hz 0.12 FitError Border Activation 0.1 2DDFT Spatial Mean Activation Central 0.08 Max First Derivative Variance k 0.06 0 5 10 15 20 25 30 35 40 Spatial Mean Activation Left number of best-ranked features Spatial Mean Activation Left Frontal Figure 7 Cross-validation error for SVM, RLDA and LPM against Laplace-Filter the number of best-ranked features. A 10-fold cross-validation was repeated 5 times and standard errors are plotted. The SVM and Mean Local Variance 15 s LPM hyperparameters were selected by an outer cross validation. 14-30 Hz The number of 6 best-ranked features was determined for building 4-7 Hz the final classifier, as the estimated error of the RLDA starts to Mean Local Variance 1 s increase significantly for higher numbers of features. Spatial Distance of Extrema Spatial Mean Activation Occipital visual analysis of these cases reveals, that most of them were mixed components that contained both, artifacts Maximum Amplitude and brain activity. Out of the 21 components which were misclassified as neural activity only two were eye Spatial Mean Activation Right movements and none were blinks. In some rare cases, Kurtosis examples which had been mislabeled by the expert could be identified. Figure 8 shows an example of a mis- 0-3 Hz classified mixed component. Deterministic Entropy To quantify the classification performance on muscle Spatial Mean Activation Frontal artifacts, we asked one expert to review the 574 artifac- tual components of the test set for muscle activity. The Variance of Local Variance 1 s expert identified 388 components which contained mus- Range cle activity (which corresponds to 67.5% of the artifac- Spatial Mean Activation Right Frontal tual components and 17.2% of all the components). Out Variance of Local Variance 15 s of the 21 artifactual components which were Mean Local Skewness 1 s 31-45 Hz Shannon Entropy Table 2 Feature weight vector and test errors. The diameter of the black circles visualizes the absolute LPM weights |E(w /|| Feature Feature Test Error Test Error w||)| per feature after learning on the training data set. (The LPM- weight RT ERP hyperparameter C had been set to C = 0.1 based on the cross-validation performance.) Current Density Norm 0.342 0.141 0.488 Range Within Pattern 0.574 0.151 0.186 Mean Local Skewness 0.317 0.309 0.442 14.1% on the test data of the RT study. The combina- 15 s tion of the six features from all three domains improves l 0.569 0.177 0.144 the error substantially compared to even the best single 8-13 Hz -0.219 0.166 0.138 feature. This shows that features which are far from FitError -0.286 0.424 0.640 optimal in single classification have a positive contribu- Combined 0.089 0.147 tion in combination with other features. 6 T Feature weights w for each feature x of the classifier H : ℝ ∋ x a sign(w ·x Looking at the complete test set of 1080 components, i i + b) Î {-1, 1} with 1 ≙ Artifact and -1 ≙ Neuronal activity. Test error (MSE) for 75 of them were misclassified as artifacts and 21 compo- the 6 single features and for the combined classification for the RT nents were misclassified as neural sources. A detailed experiment and the ERP experiment. crossvalidation error Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 11 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 Validation 2: Auditory ERP study The classifier trained on RT data and applied to 540 components of the auditory ERP study leads to an aver- age MSE of 14.7% only for the classification of artifacts (expert 1: 15.7%, expert 2: 13.7%). On average over both experts, 18 of the 540 components were misclassified as artifacts and 61.5 components were misclassifed as neural sources (expert 1: 12 - 73, expert 2: 24 - 50). Table 2 also shows the classification results for every single feature and for the combined classification for the auditory ERP data. The classification performance of the three features Range Within Pattern, l and 8-13 Hz is comparable to those in the RT experiment. They general- ize very well over different experimental setups. However, the single feature classification performance for the remaining three features, Current Density Norm, Mean Local Skewness 15 s,and FitError, was close to chance Figure 8 Example for a misclassified component.Mixed level. This does not imply, however, that these features component that combines central alpha activity and slow Mayer Waves [35]. The human expert considered the mixed component a are unimportant for the classification tasks in the com- neural source, but the classifier labeld it as an artifact. bined feature set. To asses the relevance of each feature in thecombinedfeatureset,wetrained aRLDAonthe ERP data using the labels of expert 1 and report the fea- misclassified as neuronal components, only 12 contained ture weights of the weight vector - Current Density Norm muscle activity (57.1%). This indicates that muscle arti- 0.139; Range Within Pattern 0.355; Mean Local Skewness facts were handled equally well by the classifier as other 15 s 0.255; l 0.531; 8-13 Hz -0.710; FitError 0.059. We types of artifacts. found that while the feature weight of Current Density The performance of a system on the classification task Norm and Mean Local Skewness 15 s slightly decreased has to be judged in the light of the fact that inter-expert compared to the feature vector trained on the RT data, disagreements on EEG signals are often above 10% [34]. these were still far away from zero and thus carry infor- For our data, we asked one expert to re-label the 690 mation for the classification task. components of the training set, two years after the origi- nal labeling. The MSE between the new and the former Validation 3: Application to Motor Imagery BCI rating was 13.2%. Thus, the prediction performance of Figure 9 (left) plots the BCI classification error (1-AUC) our proposed classification method was comparable to against the number of remaining independent the ranking of an human expert. Figure 9 Influence of ICA-based artifact reduction in a motor imagery BCI tested with 80 subjects. Left: box plot of classifition errors (1- AUC) against the number of remaining independent sources compared with no artifact reduction. Right: scatter plot of classification errors (1- AUC) of each subject when removing 10 source components vs. using all 30 source components. Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 12 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 components, including one entry for the standard proce- We could show that the generalization over different dure without artifact reduction. Reducing the dimen- EEG studies is possible, which is in line with the find- sionality of the data to 30 dimensions by PCA does not ings of Mognon et al. [20] who demonstrated the gener- affect BCI performance. Moreover, consecutively remov- alization of an artifact classifier to a different laboratory ing components does not impair BCI performance at and to a different paradigm. Although their method is first, as these are artifactual components according to simple and efficient, it so far does not recognize muscle the classifier. Performance breaks down only when a artifacts. Compared to the classification results of Halder et al. strict removing policy is applied and less than about 12 [22], who reported 8% of error for muscle artifacts and sources (out of ~90 original channels) are retained, which have been ranked as neural sources by the classi- 1% error for eye artifacts, the classification error of our fier. The ranking of the classifier was confirmed by a solution is slightly higher. A major difference between visual analysis of the source components. Following the the two approaches is the way the training data was ranking of very probable artifacts to less probable arti- generated. Halder et al. reported, that subjects had spe- facts, the inspection resulted in clear artifactual compo- cifically been instructed to produce a number of artifacts nents to components that contained mixtures of neural under controlled conditions for the classifier training. It and artifactual activity. can bespeculatedthatsuchatraining set contains Figure 9 (right) shows a scatter plot of classification stronger artifacts and less erroneous labels. Nevertheless, errors (1-AUC) for each subject when removing 10 the results of Halder et al. were generated based on EEG source components vs. using all 30 source components. recordings of only 16 electrodes. Without adjustments, For this soft policy for removing artifactual components it can only be applied to EEG recordings with 16 elec- the variance between subjects is very small, especially trodes. In contrast, our method is applicable to different for subjects with good classification rates. EEG setups. However, we only tested the generalization ability over different EEG studies on electrode sets that Discussion covered the whole scalp with approx. equidistant sen- To summarize, we have constructed a subject-indepen- sors. Whether the classifier is applicable to deal with dent, fast, efficient linear component classification EEG data recorded with further reduced electrode sets method that automates the process of tedious hand- remains an open question that could be analyzed in the selection of e.g. artifactual independent components. future. Theproposedmethodisapplicableonlineand gener- To assess the danger of false positives introduced by alizes to new subjects without re-calibration. It delivers our artifact detection method,weevaluated theinflu- ence of a strong artifact reduction on the classification physiologically interpretable results, generalizes well over different experimental setups and is not limited performance of a standard motor imagery BCI task. An to a specific type of artifact. In particular, muscle arti- offline analysis of data acquired from 80 healthy subjects facts and eye artifacts (besides other types) are demonstrated that removing up to 60% of the sources recognized. (that were ranked according to their artifact classifier The proposed artifact classifier is based on six care- rating) did not impair the overall BCI classification per- fully constructed features that incorporate information formance. Note that we discarded the same number of from the spatial, the temporal and the spectral domain components per subject in order to analyse the effect of of the components and have been selected out of 38 fea- false positives. In a practical BCI-system, it would prob- tures by a thorough feature selection procedure. After ably be beneficial to apply a threshold on the propability its construction on data from a reaction time experi- of being an artifact per component instead. ment, the classifier’s performance was validated on two While the suitability of our approach to remove large different data sets: (1) on unseen data of a second con- artifactualsubspaces ofthedataisawelcomeresult, an dition of the original reaction time study - here the clas- open question remains that addresses the potential per- sifier achieved a classification error of 8.9%, while formance increase by careful artifact removal. Why disagreement between two ratings of experts was 13.2%. didn’t the removal of few artifactual sources improve (2)onunseendataofanauditory oddball ERPstudy - the average motor imagery BCI performance? It is here the classifier showed a classification error of 14.7% known that CSP is rather prone to outliers if the train- in comparison to 10.6% of disagreement between ing data set is small [36]. Strategies to overcome this experts. The classification error is remarkable low given problem include the use of regularization methods for CSP (such as invariant CSP [37] or robust CSP [38]), that the second study has been recorded with half the the explicit removal of outlier trials or channels and reg- number of electrodes, under a completely different para- ularization in the following classification step. As the digm, and contained a significantly higher proportion of artifactual components. standard evaluation procedures for motor imagery data Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 13 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 contained counter measures already (channel rejection classification method seems to be sufficient when the and trial rejection based on variance), and the number feature set is carefully constructed. of training data was considerably large, the overall influ- The classification difficulties of expert raters and of ence of artifacts on the motor imagery data set probably proposed automatic classification methods reflect the was small. Furthermore, we observed, that for subjects fundamental fact that any ICA-based artifact reduction with very good motor imagery classification rates, arti- method depends crucially on the quality of the source facts did not play any role at all. We conjecture that in separation into clear artifactual and neuronal source the other subjects, artifacts either obstructed the rele- components. A good source separation method avoids mixed components that contain both, neural and artifac- vant neural activity (cases where a slight improvement by artifact removal was obtained) or artifacts played tual activity as well as arbitrary splits of a single source some role in the control of the BCI system (cases where into several components. In the following, both type of artifact removal slightly reduced the performance). errors are briefly discussed. In addition to the construction of an efficient, sparse Blind source separation is a difficult problem by itself, and interpretable classifier, our feature-selection metho- and various approaches have been proposed to solve it dology leads to valuable insights into the question of (see, e.g. [39] for a review). In the context of EEG sig- which features are best suited for the discrimination of nals, the goal is to find a source separation that mini- artifactual and neuronal source components. However, mizes the amount of mixed components. The choice of it needs to be kept in mind that while the six identified TDSEP for the pre-processing of the EEG data was features were arguably an exceptionally suitable feature motivated by the ability of the algorithm to utilize tem- set, these features were probably not the overall optimal poral structure in the data. Although this is not a choice. Furthermore, the question remains if the unique feature of TDSEP, this approach seemed to be selected features generalize to other EEG data. Single suitable for the processing of EEG data, which is com- feature classification performance drops on the ERP posed of multidimensional time series signals with tem- data for three of the six features (Current Density Norm, poral dependencies. Moreover, research indicates that Mean Local Skewness 15 s and Fit Error). However, both methods based on second-order statistics might outper- Current Density Norm and Mean Local Skewness 15 s form methods based on higher-order statistics in the carry important information in the combined classifica- removal of ocular artifacts [10,22]. Although, as Fitzgib- tion (when used together with the other four features). bon stated, “the quality of the separation is highly Still, the non-redundant information carried by the Cur- dependent on the type of contamination, the degree of rent Density Norm feature drops substantially – apro- contamination, and the choice of BSS algorithm” [9], a blem that maybe causedbythe useofafixedmatrix thorough test of various ICA methods is out of the ΓJ which had been determined on the RT setup of 115 scope of this paper. electrodes. We found no obvious explanation for the The second kind of error, the arbitrary split of sources importance change for the Fit Error feature, however. into several components, can partially be compensated In any case, several insights can be gained concerning by combining the ICA with a preceding PCA step for the construction of a suitable feature set for the classifi- the dimensionality reduction. This procedure has the cation of artifactual components in general. First, the additional advantage of removing noise in the data. We spatial, the temporal and the spectral domain of the chose to project the original data into a 30 dimensional components contain non-redundant information. Sec- space. The value of 30 was based on rough experience ond, features that quantify aspects of the pattern’s activ- and on a quick visual inspection of the data, and was ity distribution, not its single values, are discriminative. probably not the optimal choice. An improvement of Features that were ranked high in our feature selection the quality of separation might be possible by optimiz- procedure were the range within the pattern, a feature ing the dimensionality reduction, but the effort was not based on the simplicity of a source separation, features undertaken here. Future work is needed to analyze the that analyzed the spatial frequency and a binary feature influence of dimensionality reduction on source which indicates if the maximal activation is on the bor- separation. der of the pattern. Third, features that model the shape To conclude, we hope that the source component of the power spectrum as a 1/f-curve as well as the classification method presented in this study delivers a absolute spectrum in the a range are discriminative. substantial contribution for the BCI community and the Fourth, features that quantify outliers in the time series EEG community in general, as a reliable and practical such as kurtosis, entropy, and mean local skewness, tool for the removal of artifacts. To support the com- munity, to encourage the reproduction of our results, or seem to be important but redundant. We analyzed 12 allow for re-labeling of data we provide the readily such features and only one obtained a high ranking in trained classifier, an implementation of the feature the feature selection. Last but not least, a linear Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 14 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 computer interfacing to mental state monitoring. Journal of neuroscience extraction routines together with example scripts, the methods 2008, 167:82-90. extracted features of the RT data, a visualization of 1770 3. Iwasaki M, Kellinghaus C, Alexopoulos AV, Burgess RC, Kumar AN, Han YH, components together with the expert labels used for the Lüders HO, Leigh RJ: Effects of eyelid closure, blinks, and eye movements on the electroencephalogram. Clinical Neurophysiology 2005, classifier training, and a visualization of components 116(4):878-885. misclassified by our method (see Additional File 1 - 4. Goncharova II, McFarland DJ, Vaughan TM, Wolpaw JR: EMG contamination MatlabCode; Additional File 2 - TrainComponents; of EEG: spectral and topographical characteristics. Clinical Neurophysiology 2003, 114:1580-1593. Additional File 3 - TestComponents; Additional File 4 - 5. Fatourechi M, Bashashati A, Ward RK, EBirch G: EMG and EOG artifacts in Misclassifications). brain computer interface systems: A survey. Clinical Neurophysiology 2007, 118:480-494. 6. Croft RJ, Barry RJ: Removal of ocular artifact from the EEG: a review. Additional material Clinical Neurophysiology 2000, 30:5-19. 7. Makeig S, Bell AJ, Jung TP, Sejnowski TJ: Independent Component Additional file 1: MatlabCode. A Matlab implementation of the feature Analysis of Electroencephalographic Data. Advances in neural information extraction routines together with example scripts, the readily trained processing systems 1996, 8:145-151. classifier and the extracted features for all the components of the RT 8. Jung TP, Makeig S, Humphries C, Lee TW, Mckeown MJ, Iragui V, data set. Sejnowski TJ: Removing electroencephalographic artifacts by blind source separation. Psychophysiology 2000, 37:163-178. Additional file 2: TrainComponents. Visualization of the 690 9. Fitzgibbon SP, Powers DMW, Pope KJ, Clark CR: Removal of EEG Noise and independent components in the training RT data, together with the Artifact Using Blind Source Separation. Clinical Neurophysiology 2007, expert’s labels. 24(3):232-243. Additional file 3: TestComponents. Visualization of the 1080 10. Romero S, Mañanas MA, Barbanoj MJ: A comparative study of automatic independent components in the RT test data, together with the expert’s techniques for ocular artifact reduction in spontaneous EEG signals labels. based on clinical target variables: A simulation case. Computers in Biology Additional file 4: Misclassifications. Visualization of the 75 + 21 and Medicine 2008, 38:348-360. misclassified components of the RT test data 11. Crespo-Garcia M, Atienza M, LCantero J: Muscle artifact remval from human sleep EEG by using independent component analysis. Ann Biomed Eng 2008, 36:467-475. 12. McMenamin BW, Shackman AJ, Maxwell JS, Bachhuber DRW, Koppenhaver AM, Greischar LL, Davidson RJ: Validation of ICA-based Acknowledgements myogenic artifact correction for scalp and source-localized EEG. The authors thank Martijn Schreuder for his help with recording and NeuroImage 2010, 49:2416-2432. preparing the auditory ERP data set, Claudia Sannelli for her help with 13. Barbati G, Porcaro C, Zappasodi F, Rossini PM, Tecchio F: Optimization of recording and preparing the RT EEG data set, the authors of [33] for an independent component analysis approach for artifact identifaction providing the motor imagery data set. and removal in magnetoencephalographic signals. Clinical Funding by the European Community under the PASCAL Network of Neurophysiology 2004, 115:1220-1232. Excellence (IST-2002-506778) and under the FP7 Programme (TOBI ICT-2007- 14. Delorme A, Makeig S, Sejnowski T: Automatic artifact rejection for EEG 224631 and ICT-216886), by the Bundesministerium für Bildung und Forschung data using high-order statistics and independent component analysis. (BMBF) (FKZ 01IBE01A, FKZ 16SV2231, 01GQ0850 and 01IB001A) and by the Proceedings of third international independent component analysis and blind Deutsche Forschungsgemeinschaft (DFG) (VitalBCI MU 987/3-1) is gratefully source decomposition conference, San Diego, CA 2001, 457-462. acknowledged. (This publication only reflects the authors’ views. Funding 15. Romero S, Mañanas M, Riba J, Giménez S, Clos S, Barbanoj M: Evaluation of agencies are not liable for any use that may be made of the information an automatic ocular filtering method for awake spontaneos EEG signals contained herein.) Last but not least, we would like to thank our reviewers based on Independt Component Analysis. 26th Annual International for their valuable comments. Conference of the Engineering in Medicine and Biology Society (EMBS) 2004, 925-928. Authors’ contributions 16. Shoker L, Sanei S, Chambers J: Artifact Removal from The EEG studies were performed in cooperation with colleagues mentioned Electroencephalograms using a hybrid BSS-SVM algorithm. IEEE Signal in the acknowledgments section. Authors MT and SH designed and carried Processing Letters 2005, 12(10):721-724. out the RT EEG study. MT designed and carried out the ERP study. IW and 17. James CJ, Gibson OJ: Temporally constrained ICA: an application to MT designed the feature extraction and feature selection algorithms and the artifact rejection in electromagnetic brain signal analysis. IEEE Trans artifact classification method. SH provided the current density norm feature Biomed Eng 2003, 50:1108-1116. method and implementation. All other implementations were carried out by 18. Joyce CA, Gorodnitsky IF, Kutas M: Automatic removal of eye movement IW. IW and MT analyzed and evaluated the overall methodology and wrote and blink artifacts from EEG data using blind component separation. the manuscript. All authors proof-read and approved the final manuscript. Psychophysiology 2004, 41:331-325. 19. Viola FC, Thorne J, Edmonds B, Schneider T, Eichele T, Debener S: Semi- Competing interests automatic identification of independent components representing EEG The authors declare that they have no competing interests. artifact. Clinical Neurophysiology 2009, 120:868-877. 20. Mognon A, Jovicich J, Bruzzone L, Buiatti M: ADJUST: An automatic EEG Received: 29 April 2011 Accepted: 2 August 2011 artifact detector based on the joint use of spatial and temporal features. Published: 2 August 2011 Psychophysiology 2010, 1-12. 21. LeVan P, Urrestarazu E, Gotman J: A system for automatic artifact removal References in ictal scalp EEG based on independent component analysis and 1. del R, Millán J, Rupp R, Mueller-Putz G, Murray-Smith R, Giugliemma C, Bayesian classification. Clinical Neurophysiology 2006, 117:912-927. Tangermann M, Vidaurre C, Cincotti F, Kübler A, Leeb R, Neuper C, 22. Halder S, Bensch M, Mellinger J, Bogdan M, Kübler A, Birbaumer N, Müller KR, Mattia D: Combining Brain-Computer Interfaces and Assistive Rosenstiel W: Online Artifact Removal for Brain-Computer Interfaces Technologies: State-of-the-Art and Challenges. Frontiers in Neuroprosthetics Using Support Vector Machines and Blind Source Separation. 2010, 4. Computational Intelligence and Neuroscience 2007, 7(3):1-10. 2. Müller KR, Tangermann M, Dornhege G, Krauledat M, Curio G, Blankertz B: 23. Ziehe A, Laskov P, Nolte G, Müller KR: A Fast Algorithm for Joint Machine Learning for real-time single-trial EEG-analysis: From brain- Diagonalization with Non-orthogonal Transformations and its Winkler et al. Behavioral and Brain Functions 2011, 7:30 Page 15 of 15 http://www.behavioralandbrainfunctions.com/content/7/1/30 Application to Blind Source Separation. Journal of Machine Learning Research 2004, 5:801-818. 24. Blankertz B, Lemm S, Treder MS, Haufe S, Müller KR: Single-trial analysis and classification of ERP components - a tutorial. NeuroImage 2011. 25. Titchener MR: T-entropy of EEG/EOG Sensitive to Sleep State. International Symposium on Nonlinear Theory and Applications (NOLTA) 2006, 859-862. 26. Hämäläinen MS, Ilmoniemi R: Interpreting magnetic fields of the brain: minimum-norm estimates. Med Biol Eng Comput 1994, 32:35-42. 27. Haufe S, Nikulin VV, Ziehe A, Müller KR, Nolte G: Combining sparsity and rotational invariance in EEG/MEG source reconstruction. NeuroImage 2008, 42:726-738. 28. Bennett KP, Mangasarian OL: Robust linear programming discrimination of two linearly inseparable sets. Optimizaion Methods and Software 1994. 29. Chang CC, Lin CJ: LIBSVM. a library for support vector machines 2001. 30. Ledoit O, Wolf M: A well-conditioned estimator for large-dimensional covariance matrices. Journal of Multivariate Analysis 2004, 88(2):365-411. 31. Schäfer J, Strimmer K: A Shrinkage Approach to Large-Scale Covariance Matrix Estimation and Implications for Functional Genomics. Statistical Applications in Genetics and Molecular Biology 2005, 4(32). 32. Blankertz B, Tomioka R, Lemm S, Kawanabe M, Müller KR: Optimizing spatial filters for robust EEG single-trial analysis. IEEE Signal Proc Magazine 2008, 25:41-56. 33. Blankertz B, Sannelli C, Halder S, Hammer EM, Kübler A, Müller KR, Curio G, Dickhaus T: Neurophysiological predictor of SMR-based BCI performance. NeuroImage 2010, 51(4):1303-1309. 34. Klekowicz H, Malinow U, Niemcewi S, Wakarow A, Wolynczyk-Gmaj D, Piotrowski T, Durka P: Automatic analysis of sleep EEG. Frontiers in Neuroinformatics. Conference Abstract. Neuroinformatics 2008. 35. Lugaresi E, Coccagna G, Mantovani M, Lebrun R: Some periodic phenomena arising during drowsiness and sleeping in humans. Electroenceph clin Neurophysiol 1972, 32:701-705. 36. Krauledat M, Dornhege G, Blankertz B, Müller KR: Robustifying EEG data analysis by removing outliers. Chaos and Complexity Letters 2007, 2(3):259-274. 37. Blankertz B, Kawanabe M, Tomioka R, Hohlefeld F, Nikulin V, Müller KR: Invariant Common Spatial Patterns: Alleviating Nonstationarities in Brain-Computer Interfacing. In Advances in Neural Information Processing Systems 20. Edited by: Platt J, Koller D, Singer Y, Roweis S. Cambridge, MA: MIT Press; 2008:113-120. 38. Kawanabe M, Vidaurre C, Scholler S, Blankertz B, Müller KR: Robust Common Spatial Filters with a Maxmin Approach. EMBS-Conference 2009, 2470-2473. 39. Choi S, Cichocki A, Park HM, Lee SY: Blind Source Separation and Independent Component Analysis: A Review. Neural Information Processing - Letters and Reviews 2005, 6:1-57. doi:10.1186/1744-9081-7-30 Cite this article as: Winkler et al.: Automatic Classification of Artifactual ICA-Components for Artifact Removal in EEG Signals. Behavioral and Brain Functions 2011 7:30. Submit your next manuscript to BioMed Central and take full advantage of: • Convenient online submission • Thorough peer review • No space constraints or color figure charges • Immediate publication on acceptance • Inclusion in PubMed, CAS, Scopus and Google Scholar • Research which is freely available for redistribution Submit your manuscript at www.biomedcentral.com/submit

Journal

Behavioral and Brain Functions – Springer Journals

Published: Aug 2, 2011

There are no references for this article.