Original Research ARTICLE
A Riemannian Modification of Artifact Subspace Reconstruction for EEG Artifact Handling
- Neuropsychology Lab, Department of Psychology, European Medical School, Carl von Ossietzky University of Oldenburg, Oldenburg, Germany
Artifact Subspace Reconstruction (ASR) is an adaptive method for the online or offline correction of artifacts comprising multichannel electroencephalography (EEG) recordings. It repeatedly computes a principal component analysis (PCA) on covariance matrices to detect artifacts based on their statistical properties in the component subspace. We adapted the existing ASR implementation by using Riemannian geometry for covariance matrix processing. EEG data that were recorded on smartphone in both outdoors and indoors conditions were used for evaluation (N = 27). A direct comparison between the original ASR and Riemannian ASR (rASR) was conducted for three performance measures: reduction of eye-blinks (sensitivity), improvement of visual-evoked potentials (VEPs) (specificity), and computation time (efficiency). Compared to ASR, our rASR algorithm performed favorably on all three measures. We conclude that rASR is suitable for the offline and online correction of multichannel EEG data acquired in laboratory and in field conditions.
Brain-Computer Interfaces (BCIs) translate brain signals into control signals and this translation is typically based on detecting different brain states in near real-time (Zander and Kothe, 2011; Yger et al., 2015). Non-invasive BCIs often rely on electroencephalography (EEG) for brain signal acquisition. EEG has the advantage that it can be recorded with small, unobtrusive and wireless hardware, incorporating off-the-shelf smartphones. Accordingly, portable EEG technology can bring BCI applications closer to real-life applications (Gwin et al., 2010; De Vos et al., 2014; Gramann et al., 2014; Debener et al., 2015; Jungnickel and Gramann, 2016; Bleichner and Debener, 2017; Blum et al., 2017). Moreover, EEG can be acquired while walking freely outdoors (Piñeyro Salvidegoitia et al., 2019) or indoors on treadmills (Solis-Escalante et al., 2012; Wagner et al., 2012; Seeber et al., 2015). However, while such recordings allow for the investigation of human behavior during motion and in natural environments, the resulting EEG signals are contaminated with even more artifacts than traditionally recorded, stationary EEG. In our previous studies, for instance, we observed that eye movements are more prominent in outdoors mobile compared to indoors stationary conditions. In the present paper, we considered stationary and mobile EEG data, recorded indoors and outdoors. Specifically, we evaluated two versions of an EEG artifact attenuation approach which can be applied online, that is, it can clean-up EEG artifacts in near real-time, a necessity for advancing BCI applications.
A widely used approach for artifact attenuation is independent component analysis (ICA) (e.g., Makeig et al., 1996; Delorme et al., 2007). Some studies show that ICA can be used online (Akhtar et al., 2012; Hsu et al., 2014; Pion-Tonachini et al., 2018), but it is computationally demanding and designed for offline use. ICA decomposes the EEG data into a set of maximally temporally independent components each representing different brain and non-brain sources. Components representing artifacts can either be removed manually or automatically (Campos Viola et al., 2009; Bigdely-Shamlo et al., 2013; Pion-Tonachini et al., 2018). ICA may not deal with all types of artifacts equally well, this holds in particular for non-biological artifacts that are difficult to model, such as cap movement or cable movement (Campos Viola et al., 2009; Chang et al., 2018). However, these artifact sources may compromise mobile EEG signal quality in particular.
Mullen et al. (2015) proposed an EEG artifact attenuation algorithm aimed at correcting anomalies in artifact-contaminated EEG data online. The ASR approach learns a statistical model on clean calibration data and attenuates (in theory any type of) artifact by decomposing short segments of EEG data and contrasting them to calibration data in the component subspace. ASR, like other statistical anomaly detection methods (Chandola et al., 2009; Islam et al., 2016), assumes that non-brain signals induce a large amount of variance in the EEG data and can therefore be detected based on their deviant statistical properties in the principal component analysis (PCA) subspace. Similar to ICA, data containing artifacts are reconstructed using an unmixing matrix to recover EEG data. ASR has been designed for online and offline use. In both cases, the data are processed in chunks of 500 ms. As a result, ASR has a very short processing delay and relatively low computational complexity. It is well-suited for online applications and it has been shown to perform very well on EEG data in different recording conditions (Arad et al., 2018; Kim and Kim, 2018; Plechawska-Wojcik et al., 2019).
In its original implementation, ASR decomposes covariance matrices of the EEG data for the detection of artifacts using traditional Euclidean geometry. Due to their mathematical properties though, covariance matrices should be processed using methods applicable for curved, high-dimensional data spaces (Moakher and Batchelor, 2006; Yger et al., 2015; Lotte et al., 2018). For such computations, Riemannian geometry has been shown to be more precise than Euclidean methods (Barachant et al., 2012, 2013; Kalunga et al., 2015; Yger et al., 2015; Horev et al., 2017; Lotte et al., 2018) and has been used in online artifact detection and correction previously (Barthélemy et al., 2017, 2019). Therefore, we speculated that a Riemannian-adapted ASR could outperform the original ASR algorithm. In this report, we describe the implementation and evaluation of a Riemannian version of ASR, hereafter referred to as rASR. Specifically, we compared rASR to ASR and to uncorrected EEG data. As measures of artifact reduction sensitivity, we evaluated the performance by which ASR and rASR reduce eye-blink amplitudes and topographies. As measures of artifact reduction specificity, we calculated and compared visual-evoked potential (VEP) signal-to-noise ratios (SNR) and amplitudes. Moreover, computation times for each procedure were measured to assess computational efficiency.
We used 24-channel EEG data recorded from N = 27 participants performing a memory experiment similar to the one described in Piñeyro Salvidegoitia et al. (2019). EEG acquisition and experimental control were implemented on a smartphone, and artifact correction performance was evaluated for two conditions, standing and walking outdoors on the university campus, and standing and walking indoors on a treadmill.
Materials and Methods
The Original ASR Algorithm and the EEGLAB Plugin clean_rawdata
The ASR algorithm is explained in detail in Chang et al. (2018), Mullen et al. (2015), and Pion-Tonachini et al. (2018) and is available as part of the open source EEGLAB plugin clean_rawdata1. Briefly, ASR learns statistical properties of clean calibration data and compares these statistics during the processing with statistics of new data, potentially containing artifacts. It is recommended that the calibration data have a length of at least 1 min and be recorded from the respective participant during rest under comparable recording conditions as the subsequent processing. The ASR algorithm consists of two parts, the calibration and the processing. During the calibration, data are filtered, then a robust covariance matrix U, ∈ ℝc×c, with c being the number of channels, is computed using the L1-median of subsequent sample covariance matrices of the incoming data segment. The geometric median, or L1-median is defined as U = , where argmin minimizes the Euclidean distance of X to U and n denotes the amount of covariance matrices computed so far during the calibration. The single sample covariance matrices are computed as Cov (s) = X (s) X (s)T, Cov (s) ∈ ℝc×c, with c being the channel number, X (s) ∈ ℝt×c the data segment, s consisting of c channels and t samples. U is then used to compute the mixing matrix M ∈ ℝc×c such that MMT = U. While learning the statistical model, this mixing matrix M is decomposed using PCA to obtain V ∈ ℝ1×c and D ∈ ℝ1×c, the eigenvectors and eigenvalues of the mixing matrix. At the end of the calibration phase, data is projected into component space where statistical properties of each component are computed. These statistics are used to determine a threshold operator T = μ + k × σ, T ∈ ℝc×c, with μ, σ being mean and standard deviation of the components and k a tuning parameter, evaluated in Chang et al. (2018). The threshold matrix T defines the limits of normal data during the processing. During the processing, the current covariance matrix Cov(s) is averaged with the recording from the previous data segment s that was being processed. Those matrices are temporally smoothed, U(s) = f(cov(s), U(s – 1)) where U(s) ∈ ℝc×c depends on the chosen step size parameter for the data segment s and f is a weighted Euclidean average function. The amount of covariance matrices in U that are decomposed in succession depends on the window length parameter and the amount of samples currently in the method: range = sampling rate ∗ windowlength (500 ms by default). range defines the size of U : C × C × range, the update interval used for decomposition is defined by which translates to several thousand iterations of decomposition and reconstruction for every data segment in the ASR method. During this main processing loop, U is indexed in chunks of size stepsize, those U (s) ∈ ℝc×c are used to detect and reconstruct artifacts. In component space, data are reconstructed using the mixing matrix M from the calibration and a portion of the eigenvectors Vclean of U(s), namely those of the clean components determined by the threshold operator : . As a result, the artifact is removed from the data. For a complete description of implementation details, please refer to the asr_process function in the clean_rawdata plugin.
Overview of Riemannian Geometry for EEG Data
The central step in both ASR and rASR is the computation and decomposition of covariance matrices using an eigendecomposition. Data is thereby represented on a new set of axes which fulfill a statistical criterion (Clifford, 2008). In the case of eigendecomposition of covariance matrices, this criterion is variance. Covariance matrices are symmetric positive definite (SPD) matrices which lie in the space M(n) of all n x n real matrices (Lang, 1999), which is a differentiable Riemannian manifold (e.g., Barachant et al., 2012). In this space, they belong to a subset which forms a convex cone (Dattorro, 2017). Several studies have shown the beneficial effects of using Riemannian distance and average measures, as well as geometry-aware PCA methods, such as principal geodesic analysis (PGA) (Fletcher et al., 2004; Horev et al., 2017). Since the space of SPD matrices is a (negatively) curved space, the use of traditional Euclidean geometry, which implies that distances are computed along straight lines in the data space, turns out to be disadvantageous. Riemannian geometry has shown to be more efficient and yields more precise results in the analysis of EEG data (Shyu et al., 2006; Kalunga et al., 2015; Yger et al., 2015; Congedo et al., 2017; Horev et al., 2017; Barthélemy et al., 2019).
The rASR Algorithm: Geometry Aware Methods
The rASR algorithm differs from original ASR in all computations or decomposition of covariance matrices during the processing. In rASR, a sample covariance matrix is computed as a robust, unbiased estimator of the covariance matrix of the current data segment (Kalunga et al., 2015). It is defined as: where t is the number of samples in the current data segment and X ∈ ℝt×c is the current channel matrix consisting of t samples and c channels (Barachant et al., 2012; Barthélemy et al., 2019). By computing an estimator of the covariance matrix of the current data segment held in the method, the rASR method omits the necessity of computing individual covariance matrices for every small chunk of the data defined by stepsize, as described in the ASR paragraph. As a result, instead of decomposing several thousand covariance matrices, this expensive operation is only done once on the estimator covariance matrix. The reconstruction of the data segment is then done successively in small chunks, the same as in the original ASR method. Here, the stepsize parameter is set to a smaller value in rASR to ensure a good coverage of reconstruction and counterbalance the less-sensitive covariance matrix computation. To ensure a robust covariance matrix for successive data segments, the currently computed covariance matrix is, both in rASR as well as in ASR, averaged with the covariance matrix computed with the previous data segment. The running mean averaging method used in ASR potentially induces imprecisions in the result because of the so-called swelling effect, which can occur when using Euclidean geometry in non-linear spaces (Horev et al., 2017). In rASR, the method of averaging the current and previous covariance matrices was therefore exchanged by a geometry-aware averaging method, namely the Riemannian center of mass method, also known as Karcher mean (Horev et al., 2017): , where d is the geodesic (Riemannian) distance which is minimized, the analogous concept to Euclidean distance in the arithmetic mean. In rASR, the computation of the sample covariance matrices is therefore done differently and in addition, it is done less often because the sample covariance estimator is used instead of covariance matrices of many samples of the current data segment. Further, the eigendecomposition is computed using a PGA instead of a PCA, described in Fletcher et al. (2004) and Horev et al. (2017) and implemented in the Manopt toolbox (Boumal et al., 2013). PGA extends PCAs dimensionality reduction to SPD matrices. The goal in PGA is to project data into a lower dimensional subspace which best preserves the variance, analogous to PCA. In PGA, the data are projected to geodesic submanifolds which are computed such that they maximize the projected variance of the data. Distances between submanifolds are given by geodesics, the analogous of a straight line in curved space. Variance is defined as the Riemannian distance from the mean, linear subspaces from PCA are now extended to geodesic manifolds (an extensive explanation of PGA can be found in Fletcher et al. (2004)). These distances between components are especially large when data contain artifacts and component matrices differ considerably from each other, as assumed in the case of an artifactual data segment, leading to large distances on the manifold. In this case, PGA preserves more variance than PCA (Horev et al., 2017).
Experimental Design and Procedure
An unpublished data set recorded for a different purpose was used in the present report (Jacobsen, unpublished). This study followed a previous work of our group investigating the neural correlates of episodic memory formation outdoors (Piñeyro Salvidegoitia et al., 2019). The study was approved by the local ethics committee, participants provided their written informed consent prior to participation. We recorded EEG data from 27 subjects (7 male, 20 female, mean age 23 years, ± 2.5 years SD) recorded in two conditions, one indoors and one outdoors. A smartphone (“Sony Xperia Z1,” model: C6903; OS: Android 5.1.1) running the recording application Smarting (version 1.6.0, mBrainTrain, 2016 Fully Mobile EEG Devices) and a 24-channel mobile EEG cap (Easycap, Herrsching, Germany) were used to collect and record EEG data. The Presentation mobile app (Version 1.2.1, Neurobehavioral Systems Inc., Albany, CA, United States, RRID:SCR_002521) controlled experimental events. In the outdoors condition, participants were instructed to walk naturally next to the experimenter who guided them to different locations on the university campus. Indoors, participants walked on a mechanic treadmill at their preferred speed. In both conditions, words that subjects had to remember were presented on the smartphone display. The experimental design involved self-initiated presentation of those words, subjects first stopped walking and then pressed a button on the screen to start a trial. Outdoors, participants stopped walking when indicated by the experimenter, indoors, participants stopped walking when indicated by a cue on the smartphone display. Note that the data used for the analyses of the VEP elicited by the word on the screen were therefore recorded when subjects stood still and looked at the smartphone screen. Eye-blink data on the other hand were extracted from the entire indoor and outdoor dataset, containing data recorded while participants were walking and while they were standing. In the beginning of each block, 1 min of resting EEG was recorded which served as calibration data for ASR and rASR. Matlab [Version 9.0 (2016a), The Mathworks Inc., Natick, MA, United States, RRID:SCR_001622] and EEGLAB (Version 13.6.5b, Delorme and Scott, 2004, RRID:SCR_007292) were used for all analyses. The rASR toolbox was developed by the first author of this report (SB) and is freely available alongside all processing scripts used for our analyses2.
All data were recorded with a sampling rate of 250 Hz, digitized data recordings had a resolution of 24 bit. Passive Ag/AgCl electrodes were used and arranged according to the international 10/20 setup (reference: FCz, DRL: AFz). Impedances were kept below 10 kΩ. EEG data, experimental events and additional sensor readings from the smartphone and IMUs in the amplifier were recorded into an xdf file on the smartphone3. For subsequent analyses, the EEG data were low-pass filtered at 40 Hz (FIR, filter order 166) and high-pass filtered at 0.25 Hz (FIR, filter order 3300). Data were then submitted to a subset of functions from the clean_rawdata toolbox which contains the ASR algorithm (version 0.34). rASR is implemented in its own plugin which is identical to clean_rawdata in all but the two core methods asr_calibrate and asr_process. Figure 1 shows the toolbox and its wrapper and preprocessing functions for the ASR algorithm. The ASR algorithm is implemented in the core functions asr_calibrate and asr_process which are therefore the only functions that were changed to develop a Riemannian-adapted ASR algorithm. In Figure 1, adapted parts in rASR indicated by dashed outlines. During the processing, segments of raw EEG data are corrected in several steps which have all been adapted in the rASR algorithm. The covariance matrix computation and decomposition now use geometry-aware methods, the choice of clean components and the reconstruction are called in different order but have not been methodologically adapted. This rearrangement and replacement of parts of the ASR algorithm requires different parameters in the call of clean_asr for ASR and rASR, whereby rASR needs substantially more aggressive values in particular because of the scarcer update of the covariance matrix. For both ASR and rASR, we cleaned the calibration data recorded at the beginning of the experimental blocks using ASR or rASR before submitting them to processing. Data from each block were then processed using clean_flatlines, clean_drifts, and finally clean_asr using the dedicated cleaned calibration data. Since we are aiming for an online artifact cleaning solution, we refrained from using clean_channels, a powerful, yet computationally expensive and time-consuming method which cannot be used in online processing of EEG data. ASR was first used with the pre-defined default parameters, but after cleaning of the calibration data before submitting it to the processing, these parameters yielded results of clear overcorrection, so they were set to less aggressive parameters. Parameters for both methods were evaluated before starting the analysis by means of maximizing SNR and preserving morphology and are specified hereafter with name (value for rASR, value for ASR): flatline (1, 5), hp [(0.25, 0.95), (0.25, 0.95)], channel (0.9, -1), noisy (3, -1), burst (2, -1), window (0.3, 0.5), cutoff (1, 5), stepsize (16, 32), maxdims (1, 0.66). All remaining parameters were used with their default values.
Figure 1. The EEGLAB toolbox methods are listed on the left (in gray) and include the two core functions of ASR and rASR that are called by the function clean_asr. Parts of the algorithm that were adapted using geometry-aware methods are marked by dashed sections. The toolbox consists of several preprocessing and correction functions that together comprise artifact correction for multi-channel EEG data. Some of the functions marked in square brackets are computationally expensive and profit from many data, but they are typically not used for online applications, we therefore refrained from using them herein. When using the toolbox, all functions are called once, clean_asr is then calling the calibration and subsequently the processing repeatedly with short segments of the uncorrected data. The procedure is the same for the online and offline use of ASR and rASR.
Apart from the artifact attenuation step, the analysis pipeline was identical for both methods. Data were epoched according to events in the data, either stemming from experimental events (VEP) or from detected artifact events (eye-blinks). For comparison reasons, we also included data which were temporally filtered the same way as the artifact-corrected data sets but otherwise not further processed. These data are further referred to as uncorrected data.
Specificity Analysis: VEP
An artifact correction method that reduces the overall data amplitude will inevitably reduce artifact amplitudes, as well as the amplitudes of any signal of interest. In order to examine whether a signal of interest is retained in the corrected EEG data, we analyzed the amplitude and the SNR of VEPs. We compared ASR, rASR and uncorrected data for both indoors and outdoors data. We focused our analyses on VEPs elicited by words appearing on the smartphone screen. These onset/offset VEPs evoked by a static stimulus have a distinct morphology in time and space (Odom et al., 2016). We analyzed data recorded by occipital channels O1 and O2, referenced to the analog reference (located at FCz). Furthermore, every subject’s average VEP was used to determine individual indices for the latencies of the early positive deflection and the following negative deflection. The SNR was computed by dividing the ERP peak amplitude by a pre-baseline noise estimator of the standard deviation in the pre-stimulus interval (200 ms before stimulus trigger). SNR is reported in decibel (dB) throughout this paper: . To compare the specificity in retaining the signal of interest for the artifact attenuation procedures, VEP SNR values and amplitudes were statistically analyzed using a 2 × 3 repeated measures analysis of variance (ANOVA) with the factors condition (indoors, outdoors) and method (uncorrected, ASR, rASR) and followed up by paired t-tests where appropriate.
Sensitivity Analysis: Eye-Blink Reduction
In order to objectively identify eye-blinks, we used the freely available Blinker toolbox (Kleifges et al., 2017) to detect eye-related artifacts on data from the combined frontal channels Fp1, Fp2. Time indices of detected eye-blinks were added to the data sets as eye-blink events and were then used to create eye-blink epochs. The average number of eye-blinks per subject was 1165 during the outdoors condition (range: 642 to 2013) and 962 indoors (range: 361 to 1813). Blink amplitudes were statistically analyzed using a 2 × 2 repeated measures ANOVA with factors condition (indoors, outdoors) and method (ASR, rASR) to assess whether one method outperformed the other, and whether this was the case for both indoors and outdoors data. Significant ANOVA effects were followed up by paired t-tests.
In addition to the eye-blink amplitude analysis, we also evaluated eye-blink topography residual variance. Assuming spatial orthogonality between ongoing EEG and eye-blinks, a complete reduction of eye-blinks should reduce the correlation of eye-blink topographies before and after the correction to a minimum, whereas an incomplete (or over complete) correction would be indicated by residual absolute correlation between topographies of uncorrected and corrected data. Accordingly, we calculated average eye-blinks for each condition and derived the topography at eye-blink peak latency in uncorrected and corrected data sets. Pearson correlation values were Fisher z-transformed and submitted to a 2 × 2 ANOVA with the factors condition (indoors, outdoors) and method (ASR, rASR).
Specificity: Visual-Evoked Potential
Figure 2 shows the averaged single subject VEPs together with the grand averaged VEP P1 and N1 topographies for both conditions, uncorrected data and ASR and rASR corrected data. We observed typical bilateral occipital VEP topographies in corrected and uncorrected data. As can be seen, ASR and rASR did not affect VEPs very much. The similarity of VEPs between indoor and outdoor conditions complements previous reports on outdoor ERP acquisitions (Debener et al., 2012; De Vos et al., 2014) and suggests that VEPs can be obtained from visual stimuli presented on smartphone, even in otherwise uncontrolled outdoor conditions. ASR and rASR preserved the morphology of the VEP and reduced the standard deviation of the signal. The VEP maps revealed the expected early positive and late negative occipital activation in both conditions, indoors and outdoors. The amplitudes of the VEP N1 peaks were lower in ASR corrected data sets than in uncorrected and rASR corrected data sets (uncorrected: -8.5 μV indoors, -10.35 μV outdoors, ASR: -8.1 μV indoors, -9.56 μV outdoors, rASR: -10.44 μV indoors, -10.61 μV outdoors). A 2 × 3 ANOVA of the N1 peak amplitude with the two within subjects factors, condition (indoors, outdoors) and method (uncorrected, ASR, rASR) yielded a significant main effect for condition [F(1,26) = 5.33, p = 0.03, ηp2 = 0.01]. VEP amplitudes were smaller indoors (-9.11 μV) than outdoors (-10.4 μV). We found no evidence for a difference between methods [F(2,52) = 1.49, n.s.]. A comparable 2 × 3 ANOVA of the N1 SNR with two within subjects factors, condition (indoors, outdoors), and method (uncorrected, ASR, rASR) yielded no significant results, neither for the main effect condition [F(1,26) = 1.19, n.s.], nor for the main effect method [F(2,52) = 1.25, n.s.]. The interaction effect missed significance as well [F(2,52) = 2.91, p = 0.06, n.s.] and was therefore not further examined.
Figure 2. Single subject (gray lines) and grand average (black lines) visual-evoked potentials, shown for uncorrected data (left) and following correction with ASR and rASR (center column and right). VEPs are shown for indoors and outdoors data and show the arithmetic mean of occipital channels O1 and O2, referenced to a central reference near Cz. The voltage maps show the averaged single subject topographies at subject-specific latencies for the early P1 and the later N1 response, all topographies are scaled identically. Individual latencies are indicated by red (P1) and blue (N1) dots. In the uncorrected data, one channel (right temporal T8) was interpolated for the voltage map.
Sensitivity: Blink Amplitudes
Figure 3 shows the ASR and rASR comparison to uncorrected data for the analysis of the blink amplitudes. As expected, the blink amplitude of uncorrected data were orders of magnitudes larger than typical EEG signals of interest with a mean of 254.92 μV indoors and 289.02 μV outdoors. A 2 × 2 repeated measures ANOVA with two within subject factors condition (indoors, outdoors) and method (ASR, rASR) yielded a significant main effect for condition. The residual blink amplitude was significantly smaller indoors (2.33 μV) than outdoors [5.37 μV , F(1,26) = 8.28, p = 0.008, ηp2 = 0.07] but no significant difference between methods was found [F(1,26) = 2.56, p = 0.12].
Figure 3. Blink artifact topographies and time courses for all methods and conditions. Note the different y-axis scale for uncorrected and corrected data. The bold black lines show the grand average blink signal, the thin gray lines show single subject averages. Eye-blink potentials were measured as the mean between channels Fp1 and Fp2, marked as white points in the voltage maps. Blink potentials and blink voltage maps illustrate eye-blink morphologies and topographies before (left) and after (center, right) artifact correction. Smaller maps show the voltage distribution 132 ms post eye-blink maximum. In uncorrected data, the map reflects residual eye-blink activity, in corrected datasets, it is compatible with the interpretation of eye-blink generated event-related brain potentials.
In addition to the blink amplitude on the frontopolar channels alone, we also evaluated blink topographies, similar to Campos Viola et al. (2009) and Bigdely-Shamlo et al. (2013), following the notion of a multidimensional description of artifacts (Tamburro et al., 2018). Figure 4 shows similarity values (R2), reflecting the shared variance between blink topographies of uncorrected and corrected data sets. This was calculated separately for indoors and outdoors conditions and both artifact correction methods. Statistical evaluation was performed with Fisher z correlation values of corrected and uncorrected topographies. The corresponding 2 × 2 ANOVA with two within subject factors (indoors, outdoors) and method (ASR, rASR) yielded a main effect for method [F(1,26) = 10.42, p = 0.003, ηp2 = 0.08]: the average similarity value was significantly higher for the ASR method (mean R2 = 0.283) than for the rASR method (meanR2 = 0.15). The main effect of condition was not significant [F(1,26) = n.s.], however, the interaction between method and condition was significant [F(1,26) = 0.83, p = 0.004, ηp2 = 0.06], indicating a significant difference of methods outdoors [t(26) = -4.51, p = 0.00012, d = 1.03] but not indoors [t(26) = n.s.].
Figure 4. Similarity values (R2) of eye-blink topographies obtained by correlating uncorrected and corrected topographies. A low R2 value reflects very little residual topographical influence of blinks, thereby indicating better artifact suppression. Comparison of correlation values between uncorrected and ASR-corrected topographies (x-axis) and between uncorrected and rASR-corrected topographies (y-axis). The diagonal line marks the border between favorable values for rASR (lower section) and favorable values for ASR (upper section). rASR showed lower similarity values than ASR.
Following eye-blink correction we observed small post eye-blink deflections that could either indicate residual artifact or eye-blink related cortical potentials. As indicated by topographical inspection, a negative deflection approximately 132 ms post eye-blink peak latency was evident. Similar effects have been described before (Berg and Davies, 1988; Heuser-Link et al., 1992) as indicating a cortical response to the lid closure. Since this effect was unpredicted, we did not follow it up statistically.
Matlab was used to measure the processing time to correct the data sets using the different methods, once including all toolbox calls and once regarding only the core functions (i.e., the ones that were adapted). As illustrated in Figure 5, when measuring the computation time including all wrapper and preprocessing functions contained in the toolbox, rASR took 44.5% less time to process 158 data sets (27 participants with 4 to 6 blocks of recording) than ASR (mean rASR = 5.6 s, SD = 0.7 s, mean ASR = 10 s, SD = 1.5 s per data set). When only measuring the time for the core algorithm, rASR was even more efficient. This is because this measurement only includes those functions that contained the Riemannian adaptation. Here, rASR took 82% less time to process the data sets than ASR (mean time rASR = 1.4 s, SD = 0.1 s, mean time ASR = 7.7 s, SD = 1.1 s per data set). A paired t-test confirmed significant differences between the methods, both for the complete computation time [t(157) = -55.15, p < 0.01, d = 4.39] and the minimal computation time [t(157) = -77.51, p < 0.01, d = 6.17].
Figure 5. Artifact correction computation times per method and extend of used functions for the two artifact-correction procedures, measured on a laptop running Matlab 2016a, Intel®i7 (quadcore) processor and 16GB RAM. The median is shown by the line which divides the boxes, outliers are indicated by crosses outside the whiskers. The core functions contained the Riemannian adaptation and therefore show the biggest difference in computation time. When computation time for all preprocessing functions was measured, the difference was less pronounced since these functions include unchanged computations.
Fast and powerful artifact correction methods are needed for EEG applications, especially for mobile recordings. This requires a computationally efficient procedure and an approach that can deal with the complexity of EEG recorded in mobile conditions, such as walking outdoors. We have extended the original ASR algorithm with Riemannian geometry. The use of Riemannian methods in the computations involving covariance matrices proved to be beneficial, rASR was faster than ASR and corrected blink artifacts favorably to ASR. At the same time, rASR preserved the signal of interest, here the visually-evoked potential VEP. However, the following issues require consideration.
Regarding artifacts, we focused on eye-blinks and disregarded other physiological or non-physiological artifacts. The reason for this choice was that eye-blinks have a strong impact on event-related EEG analysis quality as they have a large amplitude, dominating multi-channel recordings at time of occurrence. Moreover, eye-blinks have been well characterized physiologically and tools exist that can automatically detect eye-blinks with reasonably good accuracy. We used the Blinker toolbox (Kleifges et al., 2017) to identify eye-blinks and created eye-blink artifact related epochs, in order to evaluate eye-blink correction quality. By focusing on automatically detected, we may have missed non-stereotypical blinks that were not detected by the Blinker toolbox because a coherent blink template could not be learned for them. We did not evaluate how well ASR and rASR perform on these non-stereotypical eye-related artifacts. While future work may benefit from improved blink detection procedures (Barthélemy et al., 2017; Jas et al., 2017) and the inclusion of other EEG artifacts to provide a more complete picture, the current analyses suggest that rASR implements eye-artifact correction very successfully for common eye-blinks, regardless of whether they were recorded indoors or outdoors.
Interestingly, eye-blink correction did not result in isoelectric activity, even after averaging over a larger number of eye-blink epochs. This suggests that either the correction procedure applied under- or overcorrected the artifact, or that an event-related brain potential in response to eye-blinks was found. The topographical analysis as illustrated in Figure 3 supports the latter interpretation. Blink-related event-related potentials have been reported before (Berg and Davies, 1988; Heuser-Link et al., 1992) but did not receive much attention lately. However, eye-blinks provide information about cognitive processing (Wascher et al., 2015) and, in combination with a thorough analysis of blink-generated event-related potentials, may provide a rich source of information in applied fields, such as neuroergonomics.
It is worth noting that ASR and rASR benefit from the recording of suitable calibration data. One minute of resting EEG containing few artifacts is recommended (Mullen et al., 2015). In our experience the quality of the calibration data plays an important role in the quality of ASR-based artifact reduction, and this probably also holds for the modified version of rASR. However, it was beyond the scope of this paper to compare the role of the calibration data on ASR and rASR performance. We recommend using proper calibration data, that is, calibration data taken from the same recording session and based on resting conditions. ASR performance on attenuating particular artifacts, such as free walking-related artifacts may critically depend on the same artifact being absent in the calibration data. In a future project, it would be very interesting to specifically investigate the artifact reduction capabilities of rASR on movement-related artifacts in general and walking-related artifacts in particular (Gwin et al., 2010; Oliveira et al., 2016). The mobile EEG system used in the present study features the concurrent recording of 24-channel EEG, 3D gyroscope signals from the head and 3D accelerometer signals from the smartphone, and these motion sensor signals, in particular when located at the feet, may be instrumental for gait-artifact detection evaluation (Hwang et al., 2018). The current gold standard in EEG artifact reduction seems to be ICA, which was shown in a number of studies to perform favorably when compared to other artifact reduction procedures (Hoffmann and Falkenstein, 2008; McMenamin et al., 2010) and with suitable parameters (Winkler et al., 2015). However, only few attempts were made to apply ICA online (Akhtar et al., 2012; Hsu et al., 2014; Pion-Tonachini et al., 2018). These particular approaches are computationally expensive and require careful validation of block size and sensitivity measures. To our knowledge, previously presented online EEG artifact suppression procedures appear to be more complex and computationally more demanding than ASR and rASR. Accordingly, it is not clear whether these procedures could be operated on handheld devices such as smartphones. The rASR algorithm presented here is computationally less expensive than ASR and introduces only a small processing lag. We are therefore convinced that future mobile EEG applications will benefit from efficient online artifact correction. Of course, online artifact correction may not perform well under any circumstances, therefore it should be good practice to apply state-of-the-art artifact correction offline and compare the results. This requires software solutions to save both artifact-corrected as well as uncorrected raw data. To support a widespread use of online artifact correction, the original ASR artifact correction method was implemented in Java and will be included in SCALA (Blum et al., 2017), a fully mobile framework for BCI applications on smartphone. SCALA is freely available on GitHub for download and is considered an open beta. It uses the Lab Streaming Layer (Swartz Center for Computational Neuroscience and Kothe, 2015) to receive raw data and event markers. SCALA classifies incoming time series such as EEG online and returns a classification result via LSL. The rASR algorithm presented here will be included into SCALA in the near future.
SB implemented rASR and wrote the analysis scripts. NJ collected the data in a previous study in our group. SD supervised the project during which the data were originally collected. SD and MB supervised this project. All authors contributed to the analysis of the results and to the writing of the manuscript.
This work was supported by the Cluster of Excellence “Hearing4All” EXC 1077, University of Oldenburg, Germany and by the BMBF project NeuroCommTrainer (16SV7790), Oldenburg, Germany.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The handling Editor declared a past co-authorship with one of the authors SD.
The authors thank F. Yger for discussions on PGA, F. Dehais for piloting rASR on his data, and D. Grieser and C. Neurohr for their fruitful discussions and helpful explanations regarding Riemannian manifolds.
Akhtar, M. U. T., Jung, T. P., Makeig, S., and Cauwenberghs, G. (2012). “Recursive independent component analysis for online blind source separation,” in Proceedings of the IEEE International Symposium on Circuits and Systems, (Seoul: IEEE), 2813–2816.
Arad, E., Bartsch, R. P., Kantelhardt, J. W., and Plotnik, M. (2018). Performance-based approach for movement artifact removal from electroencephalographic data recorded during locomotion. PLoS One 13:0197153. doi: 10.1371/journal.pone.0197153
Barachant, A., Andreev, A., and Congedo, M. (2013). “The Riemannian Potato: an automatic and adaptive artifact detection method for online experiments using Riemannian geometry,” in Proceedings of the TOBI Workshop, (Sion), 19–20.
Barachant, A., Bonnet, S., Congedo, M., and Jutten, C. (2012). Multiclass brain-computer interface classification by Riemannian geometry. IEEE Trans. Biomed. Eng. 59, 920–928. doi: 10.1109/TBME.2011.2172210
Barthélemy, Q., Mayaud, L., Ojeda, D., and Congedo, M. (2019). The Riemannian potato field: a tool for online signal quality index of EEG. IEEE Trans. Neural Syst. Rehabil. Eng. 27, 244–255. doi: 10.1109/TNSRE.2019.2893113
Barthélemy, Q., Mayaud, L., Renard, Y., Kim, D., Kang, S.-W., Gunkelman, J., et al. (2017). Online denoising of eye-blinks in electroencephalography. Neurophysiol. Clin. 47, 371–391. doi: 10.1016/j.neucli.2017.10.059
Bigdely-Shamlo, N., Kreutz-Delgado, K., Kothe, C., Makeig, S., Ken, K.-D., Christian, K., et al. (2013). “EyeCatch: data-mining over half a million EEG independent components to construct a fully-automated eye-component detector,” in Proceedings of the 35th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), (Osaka: IEEE), 5845–5848. doi: 10.1109/EMBC.2013.6610881
Blum, S., Debener, S., Emkes, R., Volkening, N., Fudickar, S., and Bleichner, M. G. (2017). EEG recording and online signal processing on android: a multiapp framework for brain-computer interfaces on smartphone. BioMed Res. Int. 2017:3072870. doi: 10.1155/2017/3072870
Campos Viola, F., Thorne, J., Edmonds, B., Schneider, T., Eichele, T., and Debener, S. (2009). Semi-automatic identification of independent components representing EEG artifact. Clin. Neurophysiol. 120, 868–877. doi: 10.1016/j.clinph.2009.01.015
Chang, C.-Y., Hsu, S.-H., Pion-Tonachini, L., and Jung, T.-P. (2018). “Evaluation of artifact subspace reconstruction for automatic EEG artifact removal,” in Proceedings of the 40th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), (Honolulu, HI: IEEE).
Congedo, M., Barachant, A., and Bhatia, R. (2017). Riemannian geometry for EEG-based brain-computer interfaces; a primer and a review. Brain Comput. Interfaces 4, 155–174. doi: 10.1080/2326263X.2017.1297192
De Vos, M., Gandras, K., and Debener, S. (2014). Towards a truly mobile auditory brain–computer interface: exploring the P300 to take away. Int. J. Psychophysiol. 91, 46–53. doi: 10.1016/j.ijpsycho.2013.08.010
Debener, S., Emkes, R., De Vos, M., Bleichner, M., Rodriguez-Villegas, E., and Reiner, E. (2015). Unobtrusive ambulatory EEG using a smartphone and flexible printed electrodes around the ear. Sci. Rep. 5:16743. doi: 10.1038/srep16743
Debener, S., Minow, F., Emkes, R., Gandras, K., and de Vos, M. (2012). How about taking a low-cost, small, and wireless EEG for a walk? Psychophysiology 49, 1617–1621. doi: 10.1111/j.1469-8986.2012.01471.x
Delorme, A., and Scott, M. (2004). EEGLAB: an open source toolbox for analysis of single-trial EEG dynamics including independent component analysis. J. Neurosci. Methods 134, 9–21. doi: 10.1016/j.jneumeth.2003.10.009
Delorme, A., Sejnowski, T., and Makeig, S. (2007). Enhanced detection of artifacts in EEG data using higher-order statistics and independent component analysis. NeuroImage 34, 1443–1449. doi: 10.1016/j.neuroimage.2006.11.004
Fletcher, P. T., Lu, C., Pizer, S. M., and Joshi, S. (2004). Principal geodesic analysis for the study of nonlinear statistics of shape. IEEE Trans. Med. Imaging 23, 995–1005. doi: 10.1109/TMI.2004.831793
Gwin, J. T., Gramann, K., Makeig, S., and Ferris, D. P. (2010). Removal of movement artifact from high-density EEG recorded during walking and running. J. Neurophysiol. 103, 3526–3534. doi: 10.1152/jn.00105.2010
Hwang, T.-H., Reh, J., Effenberg, A. O., and Blume, H. (2018). Real-time gait analysis using a single head-worn inertial measurement unit. IEEE Trans. Consum. Electron. 64, 240–248. doi: 10.1109/TCE.2018.2843289
Islam, M. K., Rastegarnia, A., and Yang, Z. (2016). Methods for artifact detection and removal from scalp EEG: a review. Neurophysiol. Clin. -Clin. Neurophysiol. 46, 287–305. doi: 10.1016/j.neucli.2016.07.002
Jas, M., Engemann, D. A., Bekhti, Y., Raimondo, F., and Gramfort, A. (2017). Autoreject: automated artifact rejection for MEG and EEG data. NeuroImage 159, 417–429. doi: 10.1016/j.neuroimage.2017.06.030
Kim, M., and Kim, S.-P. (2018). “A comparsion of artifact rejection methods for a BCI using event related potentials,” in Proceedings of the 6th International Conference on Brain-Computer Interface (BCI), (GangWon: IEEE), 1–4. doi: 10.1109/IWW-BCI.2018.8311530
Kleifges, K., Bigdely-Shamlo, N., Kerick, S. E., and Robbins, K. A. (2017). BLINKER: automated extraction of ocular indices from EEG enabling large-scale analysis. Front. Neurosci. 11:12. doi: 10.3389/fnins.2017.00012
Lotte, F., Bougrain, L., Cichocki, A., Clerc, M., Congedo, M., Rakotomamonjy, A., et al. (2018). A review of classification algorithms for EEG-based brain – computer interfaces: a 10 year update. J. Neural Eng. 4, 1–24. doi: 10.1088/1741-2552/aab2f2
Makeig, S., Jung, T.-P., Bell, A. J., and Sejnowski, T. J. (1996). “Independent component analysis of electroencephalographic data,” in Proceedings of the 8th International Conference on Neural Information Processing Systems, (Cambridge, MA: MIT Press).
mBrainTrain. (2016). Fully Mobile EEG Devices. Available at: https://mbraintrain.com/. (accessed April 27, 2017).
McMenamin, B. W., Shackman, A. J., Maxwell, J. S., Bachhuber, D. R. W., Koppenhaver, A. M., Greischar, L. L., et al. (2010). Validation of ICA-based myogenic artifact correction for scalp and source-localized EEG. NeuroImage 49, 2416–2432. doi: 10.1016/j.neuroimage.2009.10.010
Moakher, M., and Batchelor, P. G. (2006). “Symmetric positive-definite matrices: from geometry to applications and visualization,” in Visualization and Processing of Tensor Fields Mathematics and Visualization, eds J. Weickert and H. Hagen (Berlin: Springer), 285–298. doi: 10.1007/3-540-31272-2_17
Mullen, T. R., Kothe, C. A. E., Chi, Y. M., Ojeda, A., Kerth, T., Makeig, S., et al. (2015). Real-time neuroimaging and cognitive monitoring using wearable dry EEG. IEEE Trans. Bio-Med. Eng. 62, 2553–2567. doi: 10.1109/TBME.2015.2481482
Odom, J. V., Bach, M., Brigell, M., Holder, G. E., McCulloch, D. L., Mizota, A., et al. (2016). ISCEV standard for clinical visual evoked potentials: (2016 update). Doc. Ophthalmol. 133, 1–9. doi: 10.1007/s10633-016-9553-y
Oliveira, A. S., Schlink, B. R., Hairston, W. D., König, P., and Ferris, D. P. (2016). Induction and separation of motion artifacts in EEG data using a mobile phantom head device. J. Neural Eng. 13:036014. doi: 10.1088/1741-2560/13/3/036014
Piñeyro Salvidegoitia, M., Jacobsen, N., Bauer, A.-K. R., Griffiths, B., Hanslmayr, S., and Debener, S. (2019). Out and about: subsequent memory effect captured in a natural outdoor environment with smartphone EEG. Psychophysiology 56:e13331. doi: 10.1111/psyp.13331
Pion-Tonachini, L., Hsu, S.-H., Chang, C.-Y., and Jung, T.-P. (2018). “Online automatic artifact rejection using the real-time EEG source-mapping toolbox (REST),” in Proceedings of the IEEE Engineering in Medicine and Biology Society (EMBC), (Honolulu, HI: IEEE).
Plechawska-Wojcik, M., Kaczorowska, M., and Zapala, D. (2019). “The artifact subspace reconstruction (ASR) for EEG signal correction. a comparative study,” in Proceedings of 39th International Conference on Information Systems Architecture and Technology – ISAT 2018, eds J. Świątek, L. Borzemski, and Z. Wilimowska (Cham: Springer), 125–135. doi: 10.1007/978-3-319-99996-8_12
Seeber, M., Scherer, R., Wagner, J., Solis-Escalante, T., and Müller-Putz, G. R. (2015). High and low gamma EEG oscillations in central sensorimotor areas are conversely modulated during the human gait cycle. NeuroImage 112, 318–326. doi: 10.1016/j.neuroimage.2015.03.045
Hsu, S.-H., Mullen, T., Tzyy-Ping Jung, T.-P., and Cauwenberghs, G. (2014). “Online recursive independent component analysis for real-time source separation of high-density EEG,” in Proceedings of the 36th Annual International Conference of the IEEE Engineering in Medicine and Biology Society, (Chicago, IL: IEEE), 3845–3848. doi: 10.1109/EMBC.2014.6944462
Shyu, M.-L., Chen, S.-C., Sarinnapakorn, K., and Chang, L. (2006). “Principal component-based anomaly detection scheme,” in Foundations and Novel Approaches in Data Mining, eds T. Young Lin, S. Ohsuga, C. J. Liau, and X. Hu (Berlin: Springer-Verlag), 311–329. doi: 10.1007/11539827_18
Solis-Escalante, T. S., Wagner, J., Scherer, R., Grieshofer, P., and Müller-Putz, G. (2012). “Assessing participation during robotic assisted gait training based on EEG: feasibility study,” in Proceedings of the 3rd TOBI workshop Bringing BCIs to End-Users: Facing the Challenge-Evaluation, User Perspective, User Needs and Ethical Questions, (Wuerzburg).
Swartz Center for Computational Neuroscience, and Kothe, C. (2015). Lab Streaming Layer (LSL). Available at: https://github.com/sccn/labstreaminglayer. (Accessed April 27, 2017).
Tamburro, G., Fiedler, P., Stone, D., Haueisen, J., and Comani, S. (2018). A new ICA-based fingerprint method for the automatic removal of physiological artifacts from EEG recordings. PeerJ 6:e4380. doi: 10.7717/peerj.4380
Wagner, J., Solis-Escalante, T., Grieshofer, P., Neuper, C., Müller-Putz, G., and Scherer, R. (2012). Level of participation in robotic-assisted treadmill walking modulates midline sensorimotor EEG rhythms in able-bodied subjects. NeuroImage 63, 1203–1211. doi: 10.1016/j.neuroimage.2012.08.019
Wascher, E., Heppner, H., Möckel, T., Kobald, S. O., and Getzmann, S. (2015). Eye-blinks in choice response tasks uncover hidden aspects of information processing. EXCLI J. 14, 1207–1218. doi: 10.17179/excli2015-696
Winkler, I., Debener, S., Muller, K.-R., and Tangermann, M. (2015). “On the influence of high-pass filtering on ICA-based artifact reduction in EEG-ERP,” in Proceedings of the 37th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC), (Milan: IEEE), 4101–4105. doi: 10.1109/EMBC.2015.7319296
Yger, F., Lotte, F., and Sugiyama, M. (2015). “Averaging covariance matrices For EEG signal classification based on the CSP: an empirical study,” in Proceeding of the 23rd European Signal Processing Conference (EUSIPCO), (Nice: IEEE), 2721–2725.
Keywords: Riemann, mobile EEG, ASR, BCI, neuroergonomics, PCA, PGA
Citation: Blum S, Jacobsen NSJ, Bleichner MG and Debener S (2019) A Riemannian Modification of Artifact Subspace Reconstruction for EEG Artifact Handling. Front. Hum. Neurosci. 13:141. doi: 10.3389/fnhum.2019.00141
Received: 29 August 2018; Accepted: 11 April 2019;
Published: 26 April 2019.
Edited by:Fabien Lotte, Institut National de Recherche en Informatique et en Automatique (INRIA), France
Reviewed by:Johanna Wagner, University of California, San Diego, United States
Sylvain Chevallier, Université de Versailles Saint-Quentin-en-Yvelines, France
Copyright © 2019 Blum, Jacobsen, Bleichner and Debener. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Sarah Blum, email@example.com