Improved EOG Artifact Removal Using Wavelet Enhanced Independent Component Analysis

Electroencephalography (EEG) signals are frequently contaminated with unwanted electrooculographic (EOG) artifacts. Blinks and eye movements generate large amplitude peaks that corrupt EEG measurements. Independent component analysis (ICA) has been used extensively in manual and automatic methods to remove artifacts. By decomposing the signals into neural and artifactual components and artifact components can be eliminated before signal reconstruction. Unfortunately, removing entire components may result in losing important neural information present in the component and eventually may distort the spectral characteristics of the reconstructed signals. An alternative approach is to correct artifacts within the independent components instead of rejecting the entire component, for which wavelet transform based decomposition methods have been used with good results. An improved, fully automatic wavelet-based component correction method is presented for EOG artifact removal that corrects EOG components selectively, i.e., within EOG activity regions only, leaving other parts of the component untouched. In addition, the method does not rely on reference EOG channels. The results show that the proposed method outperforms other component rejection and wavelet-based EOG removal methods in its accuracy both in the time and the spectral domain. The proposed new method represents an important step towards the development of accurate, reliable and automatic EOG artifact removal methods.


Introduction
Electroencephalography (EEG) is a non-invasive method for measuring brain activity. Due to its low cost and high temporal resolution, it is routinely used in clinical diagnostics, epilepsy surgery and cognitive psychology research. A major concern when processing EEG measurement data is the presence of various artifacts that are generated by extra-cerebral sources, such as eye blinks and eye movements (electrooculographic/EOG artifacts), muscle movement (neck, jaw and face muscles; electromyogram/EMG artifact) or heart-related EEG disturbances (electrocardiography/ECG artifact and pulse artifact). Unfortunately, these artifacts distort the measured EEG signals and, in the worst case, can make entire measurement datasets unusable. Artifact removal is therefore an essential step in correct EEG data pre-processing [1].
The simplest artifact removal approach is to discard artifact contaminated data segments from the measurement, and process only the remaining clean segments. This approach, however, requires visual data inspection as well as manual rejection of artifact contaminated data segments or epochs. Besides being a very labor-intensive task that requires a trained expert, this method cannot be automated.

Independent Component Analysis
Independent component analysis (ICA) was originally developed to solve the blind source separation (BSS) problem [28] and normally refers to a class of algorithms that can recover statistically independent signals (components) from a linear mixture, based on higher-order statistics as a measure of independence. ICA is considered a robust method for identifying and removing artifacts normally found in EEG signals.
A brief formal introduction is as follows. Assume N statistically independent sources, s i (t), i = 1, . . . , N. Suppose that the sources cannot be observed directly, only via N sensors that obtain N observation signals, x(t). The observed signals are mixtures of the original sources. Sensors must be spatially separated (e.g., as the electrodes on the scalp), as each sensor must measure a mixture different from the others. The mixing process than can be described as where A is the square mixing matrix (spatial weight matrix, channels × components) and W = A −1 is the "unmixing matrix" that must be obtained in order to calculate an estimateŜ t of the original sources asŜ The following restrictions apply to ICA in order to produce a solution: (i) sources must be statistically independent, (ii) sources cannot have Gaussian distribution and (iii) the mixing matrix must be invertible. The estimation ofŜ t requires pre-processing steps (dimensionality reduction, centering and uncorrelation). Various ICA variants exist due to differences in the statistical measures used in the separation step [29,30] One of the most popular ICA variants in the EEG community is the Infomax ICA algorithm [31]. Infomax ICA uses a contrast function based on neural network theory and maximizes the output entropy of the neural network. Assuming x as the input to the neural network with outputs φ i w T i x , where φ i is some non-linear function, the goal is to maximize the entropy of the output L 2 = H φ 1 w T 1 x , . . . , φ n w T n x using a stochastic gradient ascend algorithm. For a more detailed description of the theoretical foundations of ICA and ICA algorithm variants, their convergence properties, the quality of source separation or their runtime complexity, the interested reader is referred to the literature [2,3,28,[32][33][34][35].

EEG Datasets
Three different sets of EEG measurement data have been selected for the evaluation of our proposed method. These include publicly available datasets as well as data recorded in our laboratory.
Semi-simulated dataset: The publicly available Klados EEG dataset [36] was created for the purpose of EOG artifact removal validations; to serve as a reference dataset that can be used for comparison purposes. Data were recorded from 27 subjects (males and females), using the standard 19 electrode 10-20 layout EEG system, with sampling frequency of 200 Hz, resulting in 54 datasets. Simulated EOG artifacts were then added to the pure, artifact-free data using the following expression: where Pure_EEG i,j is the signal obtained with eyes closed (no EOG artifacts), and the V EOG and H EOG terms are the additive vertical and horizontal EOG activities. PhysioNet EEG datasets: The PhysioNet database contains brain-computer interface datasets [37,38] that were recorded during Brain Computer Interface (BCI) experiments to measure the event-related potential (ERP) of the P300 waves in a spelling experiment. Data were collected using the BioSemi Active Two EEG system, with 64 EEG electrodes and additional vertical and horizontal (VEOG, HEOG) ocular electrodes at 2048 Hz sampling rate. Laboratory resting-state dataset: We have recorded 2-3 min closed and open eye resting state EEG in our laboratory from 22 adult volunteers (males, age from 16 to 21 years). During the experiment, subjects had to sit and relax in a silent room. Data were recorded using a Biosemi ActiveTwo EEG system (f s = 2048 Hz) using 128 active electrodes arranged in the ABC radial electrode layout. The volunteers gave their written consent for participating in the experiments.

EOG Artifact Removal Algorithm
In this section, the key steps of our proposed EOG artifact removal method are shown in algorithmic form and as a flowchart in Figure 1, followed by a detailed description of each step.
Algorithm: EOG removal: Step 1 Each measured dataset is bandpass filtered (1-47 Hz, zero phase 4th order Butterworth), then re-referenced to the average reference.
Step 2 Infomax ICA is applied to the signal to estimate the source independent components.
Step 3 Automatic identification of the EOG component: the EOG component is identified based on the correlation between each component and data of each frontal EEG channel. The component with the highest correlation and above a threshold weight is selected as an EOG component.
Step 4 The identified EOG components are searched for EOG peaks.
Step 5 One-second windows are placed around the detected EOG peaks.

(a)
If the windows cover more than 60 percent of the given component, the entire component is marked for rejection. Continue at Step 7. (b) Otherwise, the EOG windows in the component are set as the target of artifact removal.
Step 6 Wavelet decomposition using Symlet sym4 [17,39,40] wavelets of five levels is applied to decompose signals in each target window to different wavelet components, and only the high frequency components are retained for the signal reconstruction process. These retained components are used in the inverse wavelet transform to reconstruct the cleaned independent component. Step 7 Using the inverse ICA process, the artifact free signals are estimated from the corrected components.

EOG Artifact Removal Algorithm
In this section, the key steps of our proposed EOG artifact removal method are shown in algorithmic form and as a flowchart in Figure 1, followed by a detailed description of each step. Algorithm: EOG removal: Step 1 Each measured dataset is bandpass filtered (1-47 Hz, zero phase 4th order Butterworth), then re-referenced to the average reference.
Step 2 Infomax ICA is applied to the signal to estimate the source independent components.
Step 3 Automatic identification of the EOG component: the EOG component is identified based on the correlation between each component and data of each frontal EEG channel. The component with the highest correlation and above a threshold weight is selected as an EOG component.
Step 4 The identified EOG components are searched for EOG peaks.
Step 5 One-second windows are placed around the detected EOG peaks. a) If the windows cover more than 60 percent of the given component, the entire component is marked for rejection. Continue at Step 7. b) Otherwise, the EOG windows in the component are set as the target of artifact removal.
Step 6 Wavelet decomposition using Symlet sym4 [17,39,40] wavelets of five levels is applied to decompose signals in each target window to different wavelet components, and only the high frequency components are retained for the signal reconstruction process. These retained components are used in the inverse wavelet transform to reconstruct the cleaned independent component. Step 7 Using the inverse ICA process, the artifact free signals are estimated from the corrected components.

Method Details
In this section, the key steps of the proposed method were described in detail. Once the input signal is filtered and the ICA process was executed, the first key task was to identify which component represents EOG artifacts (Algorithm, Step 3

Method Details
In this section, the key steps of the proposed method were described in detail. Once the input signal is filtered and the ICA process was executed, the first key task was to identify which component represents EOG artifacts (Algorithm, Step 3). The Pearson correlation R X,Y between a given independent component Y and each of the frontal channels X (shown in Figure 2) was computed based on the underlying assumption that EOG artifacts appear primarily in the frontal channels and the component describing EOG activity should have high correlation with some of these channels. The Pearson-correlation between the frontal channels and the components was computed as where σ X and σ Y are the standard deviations of channel X and component Y, respectively. Components with the highest R value are identified as candidate EOG components to be examined further. Naturally, a different set of frontal electrodes must be selected for different electrode layouts. The number of frontal channels does not affect the accuracy of the proposed method as long as there are at least two frontal channels on the forehead, one close to the left and one to the right eye. This ensures that high correlation between ocular artifacts and EOG components can be found.
Brain Sci. 2019, 9, 355 6 of 21 where and are the standard deviations of channel X and component Y, respectively. Components with the highest value are identified as candidate EOG components to be examined further. Naturally, a different set of frontal electrodes must be selected for different electrode layouts. The number of frontal channels does not affect the accuracy of the proposed method as long as there are at least two frontal channels on the forehead, one close to the left and one to the right eye. This ensures that high correlation between ocular artifacts and EOG components can be found. The candidate components were further examined for weight value distribution and only those with weights greater than a threshold were kept as EOG components. Elements of the weight vector w are defined as: where is the average weight of component j over the frontal channels, is the weight element of the mixing matrix A, K is the number of the frontal channels and N is the number of components. The distribution of values in the weight vectors was used to calculate a statistical threshold. The distributions are shown for all three datasets in Figures 3-5 as boxplots. Red crosses represent weights for the EOG components. Note that the maximum value of each distribution acts as a reliable threshold for detecting the EOG component (outliers). The candidate components were further examined for weight value distribution and only those with weights greater than a threshold were kept as EOG components. Elements of the weight vector w are defined as: where w j is the average weight of component j over the frontal channels, w ij is the weight element of the mixing matrix A, K is the number of the frontal channels and N is the number of components.     The threshold was computed from the distribution of weights and each weight vector element was tested against it to decide whether the component is in fact an EOG component: where ̅ is the weight of component Yi, w is the averaged weight vector and 3 and are the upper quartile and interquartile range, respectively [41]. The result of this step is illustrated in Figures  6 and 7. Figure 6 shows the independent components of a selected dataset. Components 1 and 2 contained EOG artifacts (blinks and eye movements, respectively). Figure 7 shows the result of the component selection and threshold application that identified the EOG components for the sample datasets 1-4. The next step in the algorithm (Step 4) was the detection of EOG peaks within the components. First a normal peak detection was performed on the component values (finding local maxima [42]), then the peaks were further examined whether they were, in fact, EOG peaks. The decision whether   The threshold was computed from the distribution of weights and each weight vector element was tested against it to decide whether the component is in fact an EOG component: where ̅ is the weight of component Yi, w is the averaged weight vector and 3 and are the upper quartile and interquartile range, respectively [41]. The result of this step is illustrated in Figures  6 and 7. Figure 6 shows the independent components of a selected dataset. Components 1 and 2 contained EOG artifacts (blinks and eye movements, respectively). Figure 7 shows the result of the component selection and threshold application that identified the EOG components for the sample datasets 1-4. The next step in the algorithm (Step 4) was the detection of EOG peaks within the components. First a normal peak detection was performed on the component values (finding local maxima [42]), then the peaks were further examined whether they were, in fact, EOG peaks. The decision whether a local maximum belongs to the set of EOG peaks P was made using the following rule containing amplitude and duration constraints. The threshold was computed from the distribution of weights and each weight vector element was tested against it to decide whether the component is in fact an EOG component: where w i is the weight of component Y i , w is the averaged weight vector and Q 3 and IQR are the upper quartile and interquartile range, respectively [41]. The result of this step is illustrated in Figures 6 and 7. Figure 6 shows the independent components of a selected dataset. Components 1 and 2 contained EOG artifacts (blinks and eye movements, respectively). Figure 7 shows the result of the component selection and threshold application that identified the EOG components for the sample datasets 1-4. Brain Sci. 2019, 9, x FOR PEER REVIEW 8 of 22  After locating the EOG peaks, target windows were placed around the peaks for EOG artifact removal (Algorithm: Step 5). A window size of 1 second duration was used, as this spans the length of the EOG artifact waveforms [43,44]). These windows would equally designate vertical-EOG (VEOG) and horizontal-EOG (HEOG) sections. Figure 8 illustrates the results of this step showing the windows marking blink and eye movement EOGs, respectively.  After locating the EOG peaks, target windows were placed around the peaks for EOG artifact removal (Algorithm: Step 5). A window size of 1 second duration was used, as this spans the length of the EOG artifact waveforms [43,44]). These windows would equally designate vertical-EOG (VEOG) and horizontal-EOG (HEOG) sections. Figure 8 illustrates the results of this step showing the windows marking blink and eye movement EOGs, respectively. The next step in the algorithm (Step 4) was the detection of EOG peaks within the components. First a normal peak detection was performed on the component values (finding local maxima [42]), then the peaks were further examined whether they were, in fact, EOG peaks. The decision whether a local maximum m k belongs to the set of EOG peaks P was made using the following rule containing amplitude and duration constraints.
where m k is the kth peak in component Y i and E{|Y i |}. is the expected value of the component vector Y i and t refers to the timestamp of peak m k . Each two consecutive selected peaks must satisfy the peak amplitude condition and the between-peak time distance of 0.5 s to correctly classify peaks as EOG artifacts. After locating the EOG peaks, target windows were placed around the peaks for EOG artifact removal (Algorithm: Step 5). A window size of 1 s duration was used, as this spans the length of the EOG artifact waveforms [43,44]). These windows would equally designate vertical-EOG (VEOG) and horizontal-EOG (HEOG) sections. Figure 8 illustrates the results of this step showing the windows marking blink and eye movement EOGs, respectively. Artifact removal was performed on the EOG components selectively, only within the target windows, using wavelet decomposition (Algorithm: Step 6). The discrete wavelet transform (DWT) of a signal f(t) is defined as where is the wavelet basis function, j is the scale parameter and k is the shift parameter. The success of EOG detection in a component is dependent on the choice of wavelet basis function [17] and the level of decomposition [45]. Several wavelet basis functions, e.g., Haar, Daubechies, coiflet and Symlet, can be used to detect and correct EOG waveforms [39,46,47]. It has been shown [47] that the Symlet wavelet family (sym2 to sym20) is the most suitable for EOG peaks and has been used successfully in several artifact removal applications. The sym-4 wavelet was selected as final basis function due to its smallest error (RMSE) between the corrected and artifact-free signals [47]. Our tests with the Symlet wavelets confirmed the same results (mean RMSE-Haar: 9.85, db4: 7.42, sym3: 7.37, sym4: 6.29, sym5: 6.54, sym6: 6.96).
The ICA component signal was decomposed into wavelet components by passing through a quadrature mirror filter performing low-pass and high-pass filtering followed by downsampling the input signal at each level of decomposition and generating the output coefficients related to lower and higher frequencies [48]. The details of this process are shown in Figure 9.  Artifact removal was performed on the EOG components selectively, only within the target windows, using wavelet decomposition (Algorithm: Step 6). The discrete wavelet transform (DWT) of a signal f (t) is defined as where ϕ is the wavelet basis function, j is the scale parameter and k is the shift parameter. The success of EOG detection in a component is dependent on the choice of wavelet basis function [17] and the level of decomposition [45]. Several wavelet basis functions, e.g., Haar, Daubechies, coiflet and Symlet, can be used to detect and correct EOG waveforms [39,46,47]. It has been shown [47] that the Symlet wavelet family (sym2 to sym20) is the most suitable for EOG peaks and has been used successfully in several artifact removal applications. The sym-4 wavelet was selected as final basis function due to its smallest error (RMSE) between the corrected and artifact-free signals [47]. Our tests with the Symlet wavelets confirmed the same results (mean RMSE-Haar: 9.85, db4: 7.42, sym3: 7.37, sym4: 6.29, sym5: 6.54, sym6: 6.96). The ICA component signal was decomposed into wavelet components by passing through a quadrature mirror filter performing low-pass and high-pass filtering followed by downsampling the input signal at each level of decomposition and generating the output coefficients related to lower and higher frequencies [48]. The details of this process are shown in Figure 9.
The ICA component signal was decomposed into wavelet components by passing through a quadrature mirror filter performing low-pass and high-pass filtering followed by downsampling the input signal at each level of decomposition and generating the output coefficients related to lower and higher frequencies [48]. The details of this process are shown in Figure 9. Different levels of wavelet decomposition were tested to find the optimal parameters. Five levels of DWT were used to decompose the component into detail (D1:D5) and approximation coefficients (A), as illustrated in Figure 10. Coefficients D1:D3 represent the higher frequency components while coefficients D4:D5 while A represent low frequency components. Since the spectrum of the EOG artifacts is concentrated in the frequencies below 7 Hz [49], the signals were reconstructed only from coefficients D1:D3, which represent the high frequencies related to the EEG signal; the other components were discarded. The reconstructed signals were then projected back to the EOG components, and inverted to obtain the artifact free data. The proposed method was able to automatically detect and correct both the vertical EOG activity (blinks) and horizontal EOG artifacts (eyes movements), which made it suitable for unsupervised artifact removal applications. Different levels of wavelet decomposition were tested to find the optimal parameters. Five levels of DWT were used to decompose the component into detail (D1:D5) and approximation coefficients (A), as illustrated in Figure 10. Coefficients D1:D3 represent the higher frequency components while coefficients D4:D5 while A represent low frequency components. Since the spectrum of the EOG artifacts is concentrated in the frequencies below 7 Hz [49], the signals were reconstructed only from coefficients D1:D3, which represent the high frequencies related to the EEG signal; the other components were discarded. The reconstructed signals were then projected back to the EOG components, and inverted to obtain the artifact free data. The proposed method was able to automatically detect and correct both the vertical EOG activity (blinks) and horizontal EOG artifacts (eyes movements), which made it suitable for unsupervised artifact removal applications. Figure 10. Wavelet decomposition of a target EOG peak signal window within an independent component.

Performance Metrics
The quality of artifact removal methods can be quantified by two basic types of metrics; metrics that describe the amount of artifact removed by a given cleaning method, and metrics that measure the distortion introduced in the signal by the cleaning process [50]. Two metrics of the first type are the artifact removal percentage and the signal-to-noise ratio difference [51]. When the true, uncontaminated EEG and the added artifact signals are known, the artifact removal percentage can be calculated as where is the autocorrelation of the true EEG signal with time lag 1, is correlation between the true EEG and the cleaned signals, while is the correlation between the true EEG and the artifactual signals. When is close to the reference , the negative term tends to 0, Figure 10. Wavelet decomposition of a target EOG peak signal window within an independent component.

Performance Metrics
The quality of artifact removal methods can be quantified by two basic types of metrics; metrics that describe the amount of artifact removed by a given cleaning method, and metrics that measure the distortion introduced in the signal by the cleaning process [50]. Two metrics of the first type are the artifact removal percentage λ and the signal-to-noise ratio difference [51]. When the true, uncontaminated EEG and the added artifact signals are known, the artifact removal percentage can be calculated as where R re f is the autocorrelation of the true EEG signal with time lag 1, R cleaned is correlation between the true EEG and the cleaned signals, while R contam is the correlation between the true EEG and the artifactual signals. When R cleaned is close to the reference R re f , the negative term tends to 0, hence a high lambda value indicates high efficacy in artifact removal. The difference in signal-to-noise ratio ∆SNR [51] is a similar measure characterizing the amount of artifact removed from the signals. It is defined as where σ 2 x is the variance of the true EEG signal, and σ 2 e contam and σ 2 e cleaned are the variances of the error signals e contam (n) = r(n) − x(n) and e cleaned (n) = r (n) − x(n) with x(n), r(n) and r (n) representing the true EEG, contaminated and the artifact cleaned signals, respectively. Distortion in the time domain can be quantified using the root mean square error calculated between the true EEG x(n) and the cleaned signals r (n).
Spectral distortion can be measured by the magnitude squared coherence (MSC) [52] that computes the frequency-domain correlation between the pure and the cleaned EEG signals: where R xy ( f ) is the cross spectral density between the two signals x and y at frequency f , and R xx ( f ) and R yy ( f ) are the autospectral density of x and y, respectively. MSC is a frequently used metric for evaluating frequency-related distortions after artifact removal [6,[53][54][55][56][57].

Results
This section presents the performance evaluation of the proposed EOG removal method. Three datasets mentioned in Section 3 were used; the Klados, the PhysioNet and the laboratory resting-state datasets. For each dataset, the proposed method (PM) is compared to the traditional full component rejection method (ICArej) [58] and the wavelet-enhanced ICA (wICA) [6] component correction methods using the performance metrics specified in Section 3.4. wICA is also compared to rejection ICA to confirm its claimed higher performance.

Semi-Simulated EEG Dataset
The performance of the proposed method was first evaluated on the Klados datasets [36]. These measurements contain semi-simulated signals, containing resting-state measured signals with and without added simulated EOG contamination. Access to the pure EEG signal allows for calculating accurate performance metrics. For illustrative purposes, Figure 11 shows the contaminated and pure EEG signals, as well as the absolute difference between the wICA-cleaned signal and the pure EEG, and the difference of the signal cleaned with the proposed method and the pure EEG signal. Note that the amplitude scales are different in order to make the difference signals visible. The contaminated segment shows three strong blink (ch [1][2][3][4][17][18][19] and two eye movement (ch [11][12] artifacts. Note the difference between the difference signals (wICA-EEG true , PM-EEG true ) obtained after cleaning with the wICA and the proposed method. The high-frequency content in the wICA difference signal indicates the removal of non-EOG signal components. Figure 12 shows a zoomed-in section of dataset12 (channel Fp1) illustrating how the PM cleaning method leaves the EEG signal intact outside the EOG zones, and how it follows the true EEG within the zones. The figures qualitatively indicate the improved removal quality of the proposed method.  In addition to the statistical analysis, for enabling side-by-side comparison with the wICA method, Table 1 lists the RMSE values for the exact same datasets and channels that were reported in [6]. better (19.1%, p = 0.00236) than the wICA and 32.6% better (p = 1.43 × 10 −5 ) than the reject ICA methods. With respect to the Δ metric, the wICA method was significantly better than the reject ICA method (50.05%, p = 2.08 × 10 −5 ). The proposed method, however, resulted in significantly increased SNR compared to wICA (79.5%, p = 7.78 × 10 −15 ) and reject ICA better (169.34%, p = 7.96 × 10 −36 ). The RMSE results were similarly positive; wICA improved upon reject ICA by 39.1% (p = 3.89 × 10 −18 ), while the proposed method showed 36.32% improvement (p = 5.84 × 10 −33 ) over the wICA and 61.22% over the reject ICA (p = 5.80 × 10 −9 ) methods, reducing the average RMSE from 5.579 (wICA) to 3.553 V. Figure 11. Illustration of the cleaning performance on one artifact contaminated section of the Klados dataset9. The two subplots on the right show the difference of the pure electroencephalography (EEG) data and the wavelet-enhanced ICA (wICA) and proposed method (PM) cleaned signals, respectively. Amplitude scales are different to make difference signal visible.  A quantitative statistical comparison was performed on the entire dataset (54 measurements), in which the λ, ∆SNR, RMSE and MSC metrics were computed for each channel in each dataset with the three removal methods (rejection ICA, wICA and the proposed method) under study. After the performing channel, the distributions of the metrics for the dataset population are shown in Figure 13. Each metric value set was checked for normality and equal variance (F-test). A two-sample t-test (α = 0.05) was performed to decide whether there is a significant difference in performance between the PM and the wICA/ICArej methods for any metric. Performance of the wICA with respect to the rejection ICA method was also examined to verify claims that wICA outperforms rejection-based removal.
In addition to the statistical analysis, for enabling side-by-side comparison with the wICA method, Table 1 lists the RMSE values for the exact same datasets and channels that were reported in [6]. In addition to the statistical analysis, for enabling side-by-side comparison with the wICA method, Table 1 lists the RMSE values for the exact same datasets and channels that were reported in [6]. While the RMSE result indicates improved removal quality in the time domain, a key question remained as to how the spectral characteristics of the signal change after cleaning. Figure 14 illustrates the effect of artifact removal on the power spectral density of the EEG signals. The frontal channel Fp1 of dataset12 was used to show the difference among the difference methods. Note how the contaminated signal introduces strong − frequency band distortions. The reject ICA and wICA methods decreased this low frequency distortion but introduced higher and band frequency power increase. The proposed method, on the other hand, removed low frequency artifact-related distortions and followed the power density distribution of the pure EEG signal for higher frequencies with very little error.
Performing the analysis for the entire dataset, the magnitude squared coherence after cleaning with the proposed method was 13.69% better (p = 4.20 × 10 −8 ) than the wICA results and 15.93% better (p = 3.91 × 10 −8 ) than the reject ICA values. No significant difference was found between the wICA and reject ICA results (p = 0.335). Figure 15 shows the overall grand average MSC results for the three methods. The performance advantage of our proposed method over the rejection ICA and wICA methods was clearly demonstrated.  While the RMSE result indicates improved removal quality in the time domain, a key question remained as to how the spectral characteristics of the signal change after cleaning. Figure 14 illustrates the effect of artifact removal on the power spectral density of the EEG signals. The frontal channel Fp1 of dataset12 was used to show the difference among the difference methods. Note how the contaminated signal introduces strong δ − θ frequency band distortions. The reject ICA and wICA methods decreased this low frequency distortion but introduced higher α and β band frequency power increase. The proposed method, on the other hand, removed low frequency artifact-related distortions and followed the power density distribution of the pure EEG signal for higher frequencies with very little error.
Performing the analysis for the entire dataset, the magnitude squared coherence after cleaning with the proposed method was 13.69% better (p = 4.20 × 10 −8 ) than the wICA results and 15.93% better (p = 3.91 × 10 −8 ) than the reject ICA values. No significant difference was found between the wICA and reject ICA results (p = 0.335). Figure 15 shows the overall grand average MSC results for the three methods. The performance advantage of our proposed method over the rejection ICA and wICA methods was clearly demonstrated. Figure 16 shows, for a selected single frontal channel (Fp1, dataset12), the magnitude squared coherence in order to compare the spectral accuracy of the different EOG removal methods in a non-averaged manner. The results indicate that the different EOG artifact cleaning methods produced different spectral distortion in frequencies below 7 Hz. Coherence was the lowest for the uncleaned, EOG contaminated signal. The rejection-based ICA and wICA methods both reduced this distortion, but it was our proposed method that produced coherence values closest to the ideal value of 1. Note that wICA also introduced a slight distortion in the 7-17 Hz range as well, which might be the result of unnecessary removal of higher frequency wavelet components.  Figure 16 shows, for a selected single frontal channel (Fp1, dataset12), the magnitude squared coherence in order to compare the spectral accuracy of the different EOG removal methods in a nonaveraged manner. The results indicate that the different EOG artifact cleaning methods produced different spectral distortion in frequencies below 7 Hz. Coherence was the lowest for the uncleaned, EOG contaminated signal. The rejection-based ICA and wICA methods both reduced this distortion, but it was our proposed method that produced coherence values closest to the ideal value of 1. Note that wICA also introduced a slight distortion in the 7-17 Hz range as well, which might be the result of unnecessary removal of higher frequency wavelet components.      (Fp1, dataset12), the magnitude squared coherence in order to compare the spectral accuracy of the different EOG removal methods in a nonaveraged manner. The results indicate that the different EOG artifact cleaning methods produced different spectral distortion in frequencies below 7 Hz. Coherence was the lowest for the uncleaned, EOG contaminated signal. The rejection-based ICA and wICA methods both reduced this distortion, but it was our proposed method that produced coherence values closest to the ideal value of 1. Note that wICA also introduced a slight distortion in the 7-17 Hz range as well, which might be the result of unnecessary removal of higher frequency wavelet components.     Figure 16 shows, for a selected single frontal channel (Fp1, dataset12), the magnitude squared coherence in order to compare the spectral accuracy of the different EOG removal methods in a nonaveraged manner. The results indicate that the different EOG artifact cleaning methods produced different spectral distortion in frequencies below 7 Hz. Coherence was the lowest for the uncleaned, EOG contaminated signal. The rejection-based ICA and wICA methods both reduced this distortion, but it was our proposed method that produced coherence values closest to the ideal value of 1. Note that wICA also introduced a slight distortion in the 7-17 Hz range as well, which might be the result of unnecessary removal of higher frequency wavelet components.

Resting State EEG Dataset
To evaluate the performance of our method on real EEG data, 2-3 min long 128-channel resting state EEG measurements of 10 subjects (obtained in our laboratory) were used. Since the true, artifact-free EEG signals are unknown in this case, modified performance metrics were used. The true EEG signal was estimated for each subject from a manually selected 5-s long artifact-free segment.
The datasets were then cleaned with the three different methods, and partitioned into 5-s long segments. The performance metrics were subsequently calculated by using the entire signal (all 5-s segments) with respect to the reference segment in the corresponding formulae. The distribution of the results for each method is shown in Figure 17. While the range of values are lower (λ and ∆SNR) or higher (RMSE) than those obtained for the semi-simulated Klados dataset (due to the different estimation of the true EEG), the trend in performance was the same. lower ( and Δ ) or higher ( ) than those obtained for the semi-simulated Klados dataset (due to the different estimation of the true EEG), the trend in performance was the same.
Similar results were obtained for the spectral distortion. As shown in Figure 18, the Magnitude Squared Coherence (MSC) values for the proposed method were significantly higher than for the ICArej and wICA methods.
As a qualitative illustration of the effect of our removal method on real measurement data, Figure 19 shows a 20-s section of the contaminated resting state EEG before and after EOG removal (proposed method). Figure 20 illustrates the same effect on a 2D scalp potential map at the peak of an EOG artifact. The EOG artifact was clearly visible in the frontal area that disappeared after cleaning. Note also the emerging parietal topography in the cleaned version, which was almost completely hidden in the contaminated map. Brain Sci. 2019, 9, x FOR PEER REVIEW 16 of 22 Figure 18. MSC values obtained with different cleaning methods for the resting state laboratory dataset (20 subjects).
As a qualitative illustration of the effect of our removal method on real measurement data, Figure 19 shows a 20-second section of the contaminated resting state EEG before and after EOG removal (proposed method). Figure 20 illustrates the same effect on a 2D scalp potential map at the peak of an EOG artifact. The EOG artifact was clearly visible in the frontal area that disappeared after cleaning. Note also the emerging parietal topography in the cleaned version, which was almost completely hidden in the contaminated map.

Peak Detection Performance
The accuracy of peak detection is crucial in the proposed method. Since the Klados and PhysioNet datasets contain annotations for EOG events, these were used to verify the performance of our EOG peak detection approach. Peak detection performance was characterized by the sensitivity measure, The PhysioNet P300 dataset was originally created to detect and classify P300 peaks in the BCI  As a qualitative illustration of the effect of our removal method on real measurement data, Figure 19 shows a 20-second section of the contaminated resting state EEG before and after EOG removal (proposed method). Figure 20 illustrates the same effect on a 2D scalp potential map at the peak of an EOG artifact. The EOG artifact was clearly visible in the frontal area that disappeared after cleaning. Note also the emerging parietal topography in the cleaned version, which was almost completely hidden in the contaminated map.

Peak Detection Performance
The accuracy of peak detection is crucial in the proposed method. Since the Klados and PhysioNet datasets contain annotations for EOG events, these were used to verify the performance of our EOG peak detection approach. Peak detection performance was characterized by the sensitivity measure,

Artifact Removal Performance
The PhysioNet P300 dataset was originally created to detect and classify P300 peaks in the BCI speller experiment [37,38] and as such, can be used to examine the proposed method for cleaning task-oriented event related potential data. Two tests were conducted to verify whether or not the cleaning methods distort ERP waveforms and peaks. First, a statistical analysis was performed on the RMSE values to verify the presence of significant improvements; second, the distortion effects of the different cleaning methods were examined.
For the statistical analysis, from among the target and non-target epochs, the target epochs were selected that elicit the P300 component. These resulted in 21 stimulus-locked epochs of length 500 ms extracted from the original contaminated 64-channel measurements for each subject (subjects s03, s04, s08 and recording rc02). From these 21 epochs, the artifact free epochs were selected and averaged for estimating the reference, pure P300 ERP signal ERP re f (number of epochs varied from 16 to 19) and averaged to generate a pure reference ERP signal. The contaminated P300 (ERP contam ) was computed by averaging the 21 uncleaned epochs. Then, the original recordings were cleaned with the three removal methods in question (rejection ICA, wICA and PM), and an ERP signal for each method was generated by averaging the 21 segments of the cleaned signals resulting in ERP rej clean , ERP wICA clean and ERP PM clean . Since the ERP waveforms differ from channel to channel, the channels were not averaged to calculate group statistics. Instead, subjects were selected individually then a statistical test was performed using the 64 channel-ERPs as sample population for pairwise comparison or the removal methods. The two-sample t-tests for each subject produced the results shown in Table 2. The proposed method performed significantly better than the wICA or rej ICA methods for each subject. Similar outcome is obtained for the λ and ∆SNR performance metrics. The distribution of the results for each method is shown in Figure 21.  In the second test, all epochs of the cleaned signals were used to compute the ERP signal and compared to the pure reference ERP. This shows the ability of the method to recreate the pure ERP waveforms after artifact removal. Figure 22.b shows the ERP distortion under real conditions, in presence of EOG artifacts. The solid blue line indicates the contaminated P300 signal. The proposed method follows the pure ERP curve with the smallest error, see inset in Figure 22.b. The distortion of the removal methods was tested by two ways. First, the pure ERP signal was compared to the cleaned ERP signals averaged from the same epochs as the pure ERP (artifactual epochs excluded). This shows the distortion of each method operating on artifact-free data ( Figure 22). The rej ICA and wICA introduce larger distortions, since the entire signal is affected by EOG removal, even if only artifact-free epochs are averaged afterwards. By using the proposed method, however, artifact-free sections of the signal are unaffected and the averaged clean epochs are nearly identical to the reference signal. See inset in Figure 22a. In the second test, all epochs of the cleaned signals were used to compute the ERP signal and compared to the pure reference ERP. This shows the ability of the method to recreate the pure ERP waveforms after artifact removal. Figure 22.b shows the ERP distortion under real conditions, in presence of EOG artifacts. The solid blue line indicates the contaminated P300 signal. The proposed method follows the pure ERP curve with the smallest error, see inset in Figure 22.b.

Discussion
EOG artifacts are random, high-amplitude distortions in EEG recordings that, if appear frequently, can make entire measurements unusable. Due to the unpredictable nature of artifacts, traditional artifact removal is based on manually data inspection and rejection of contaminated data segments. This process is both time-consuming and prone to human errors. The introduction of independent component analysis for artifact removal [11] revolutionized the field, first by providing a theoretical framework for separating artifacts, then secondly, by paving the way to automatic, intervention-free implementations. Unfortunately, the strong statistical independence assumption of ICA does not always hold in practice, resulting in neural data leaking into artifact components. In these cases, independent component rejection based artifact removal methods lose valuable neural activity information. In the second test, all epochs of the cleaned signals were used to compute the ERP signal and compared to the pure reference ERP. This shows the ability of the method to recreate the pure ERP waveforms after artifact removal. Figure 22b shows the ERP distortion under real conditions, in presence of EOG artifacts. The solid blue line indicates the contaminated P300 signal. The proposed method follows the pure ERP curve with the smallest error, see inset in Figure 22b.

Discussion
EOG artifacts are random, high-amplitude distortions in EEG recordings that, if appear frequently, can make entire measurements unusable. Due to the unpredictable nature of artifacts, traditional artifact removal is based on manually data inspection and rejection of contaminated data segments. This process is both time-consuming and prone to human errors. The introduction of independent component analysis for artifact removal [11] revolutionized the field, first by providing a theoretical framework for separating artifacts, then secondly, by paving the way to automatic, intervention-free implementations. Unfortunately, the strong statistical independence assumption of ICA does not always hold in practice, resulting in neural data leaking into artifact components. In these cases, independent component rejection based artifact removal methods lose valuable neural activity information.
The wavelet-enhanced ICA (wICA) [6] method showed that independent components can be cleaned from artifacts if they are not rejected entirely. As shown in the Results section, correcting components this way not only preserves information, but also reduces distortions that rejection ICA methods introduce in the time and frequency domain. Distortions in the frequency domain, for instance, can corrupt EEG-based connectivity analyses [6].
The novelty of the method proposed in this paper was that component artifact correction was only performed in EOG contaminated sections of the component, ensuring that non-EOG contaminated sections were left untouched. The statistical analysis of the artifact removal performance metrics confirmed that while wICA outperformed rejection ICA methods in most performance parameters our proposed method significantly outperformed the quality of both the wICA and the rejection ICA EOG cleaning methods, both in the time and spectral domains, resulting in close-to-ideal pure EEG signals.

Conclusions
This paper described an improved wavelet-based ICA method for removing EOG components from EEG measurements. The method operated on independent components produced by an independent component analysis, automatically selected EOG-contaminated components for subsequent wavelet decomposition. EOG peaks were detected in the selected components, then the wavelet components representing EOG artifact waveforms were removed in windows placed around the EOG peaks. The component was then reconstructed from the remaining wavelet coefficients and used in signal reconstruction using the inverse ICA process. This partial component cleaning approach significantly outperformed the popular wICA and rejection ICA based artifact removal methods in all key artifact removal performance metrics. In addition, our method was fully automatic; it did not require manual component and artifact inspection, which could simplify and speed up high-quality artifact removal processes. Our future research will focus on implementing this method using high-performance computing techniques to support very fast and potentially online artifact cleaning.