Artifacts in EEG and how to remove them: ATAR, ICA
Electroencephalogram (EEG) is an approach to record the electrical activities of the brain using sensors placed on the scalp. While it is very fascinating to work with EEG signals, since they represents the brain activity during any task, it is also difficult to interpret them just by looking at them. With experience, you could learn to see a few things in the given EEG signal but not very complex phenomenon. The most common issue working EEG signal is to remove the artifacts. Artifacts in EEG are the signals that are caused by non-brain activity, such as eye movements, muscle contraction, body movement, heart beating etc. To reliably extract and interpret the EEG events, it is therefore, crucial to remove the artifacts. There are many approaches that have been proposed in the literature to remove the selective artifacts. In this post we only discuss two approaches (1) ICA based (2) ATAR. Most of the details are from my published work: Automatic and tunable algorithm for EEG artifact removal using wavelet decomposition [1]
The objective of this post is show how to remove artifact from EEG using python library — spkit [2]. We will go into a very brief idea about both approaches. For more details, the reference [1] is suggested.
ICA Based:
One of the widely used approach is Blind Source Separation (BSS), as name suggests, given multi channel EEG, BSS recovers the source signals without having any ground truth, hence the word 'blind'. It is done using Independent Component Analysis (ICA), which decompose multi channel observed signals into multi channel source signals, such that source signals are statistically independent. In simple words, given N-channel observed (recorded) EEG signals, ICA will transform them to different set of N-channel signals, assumed to be source signals, which of them share no information (hence statistically independent). After that, any source signal that have attributes of artifacts are removed (by multiplying with 0) and transform back to EEG.
Though, ICA is effective approach to separate sources, however there are a few limitations with it. (1) Need of multi-channel: ICA can work only on multi-channel, for given N-channel, it will compute N-sources, (2) Effective convergence of ICA: Since this ICA works on minimizing MI among source signals, this is an optimization problem, there is always an issue with convergence. There might not be independent source signals. (3) Aggressive: For EEG Artifact, a source component which is assumed to be artifact, is removed and EEG is computed back. This makes ICA approach aggressive (for details check explanation in [1]. (4) While removing any components, it induces extra frequency components, which can be observed in spectrum (periodogram plot). In time domain, these extra frequency components are like some jerky effect, something like sudden change. We will observe these all things in the coding section.
Aggressiveness of any artifact removal algorithm for EEG is very destructive for the predictive modeling. Since we do not know, if algorithm is really removing the artifactual component or some important neural information along with assumed component that might be more useful for the predictive modeling.
ATAR Algorithm — Wavelet Based
Due to characteristics of wavelet, capturing short temporal patterns, there have been many artifact removal algorithms based on wavelet. One major benefit is that most of these approaches operates on single channel, do not require multi-channel EEG data. And there is no question of convergence.
Though most of approaches based on wavelet have fixed statistical threshold, which makes them aggressive to remove components, however, there is a way to compute threshold on different parameters, and tune them to control aggressiveness of the algorithm. ATAR is precisely based on that idea (from my published work). ATAR is based on Wavelet Package Decomposition, and Automatic and Tuneble Artifacts Removal Algorithm, where you have control over the aggressiveness of the algorithm. ATAR algorithm comes with three different modes of operations — Soft Thresholding, Linear Attenuation and Elimination and a controlling parameter β. All three modes also help. While Elimination mode is like any other algorithm — wavelet component is completely removed if it fits the criteria, on the other hand, Linear Attenuation, attenuate the magnitude of presumed artefactual, and Soft Thresholding rather than removing it decreases (clips) the magnitude of wavelet component to get out of the criteria. An example of tuning the β is showing in figure below as animation for three different modes, as β increases, the presumed artefactual components are removed more aggressively.
Uniqueness of ATAR lies in three aspects; most importantly (1) It’s ability to tune the aggressiveness, which helps in efficient predictive modelling, (2) controlling the adaptiveness of threshold on each window and mode of operation, and (3) it is based on wavelet, which capture more localized pattern and does not harm entire segment of signal.
Coding: How to remove using spkit
In this section, we will show how to remove artifacts with both, ICA and ATAR algorithm. For that you would need python library — spkit, along with numpy and matplotlib. To install spkit use pip and make sure you have spkit with version greater than equal to 0.0.9.4:
pip install spkitLet’s first import libraries and get an example of 16 sec EEG signal
import numpy as np
import matplotlib.pyplot as pltimport spkit as sp
from spkit.data import load_data
sp.__version__
It should be 0.0.9.4
X,ch_names = load_data.eegSample()
fs = 128Xf = sp.filter_X(X,band=[0.5], btype='highpass',fs=fs,verbose=0)
print(Xf.shape)t = np.arange(Xf.shape[0])/fs
plt.figure(figsize=(12,4))
plt.plot(t,Xf+np.arange(-7,7)*200)
plt.xlim([t[0],t[-1]])
plt.xlabel('time (sec)')
plt.yticks(np.arange(-7,7)*200,ch_names)
plt.grid()
plt.title('Xf: 14 channel - EEG Signal (filtered)')
plt.show()
ICA Based
Let’s apply ICA based approach:
Following code apply ICA based algorithm on each window (128 samples = 1sec) with 50% overlapping to reconstruct 16 sec signals with overlap add method.
XR = sp.eeg.ICA_filtering(Xf.copy(),verbose=1,winsize=128)
print(shape: XR.shape)ICA Artifact Removal : extended-infomax 100.05%|#############################|2112\2113| Done…
shape: (2048,14)
Let’s plot as see the results:
plt.figure(figsize=(15,10))
plt.subplot(221)
plt.plot(t,Xf+np.arange(-7,7)*200)
plt.xlim([t[0],t[-1]])
#plt.xlabel('time (sec)')
plt.yticks(np.arange(-7,7)*200,ch_names)
plt.grid()
plt.title('X: Filtered signal',fontsize=16)
plt.subplot(222)
plt.plot(t,XR+np.arange(-7,7)*200)
plt.xlim([t[0],t[-1]])
#plt.xlabel('time (sec)')
plt.yticks(np.arange(-7,7)*200,ch_names)
plt.grid()
plt.title('XR: Corrected Signal',fontsize=16)
#plt.show()
#plt.figure(figsize=(12,5))
plt.subplot(223)
plt.plot(t,(Xf-XR)+np.arange(-7,7)*200)
plt.xlim([t[0],t[-1]])
plt.xlabel('time (s)')
plt.yticks(np.arange(-7,7)*200,ch_names)
plt.grid()
plt.title('Xf - XR: Difference \n(removed signal)',fontsize=16)
# #plt.show()
plt.subplot(224)
plt.plot(t,Xf[:,0],label='Xf')
plt.plot(t,XR[:,0],label='XR')
#plt.plot(t,Xf[:,0]-XR[:,0])
plt.xlim([t[0],t[-1]])
plt.xlabel('time (s)')
plt.grid()
plt.legend()
plt.tight_layout()
plt.show()
It can be seen that ICA is able to remove that very peaky component, however, it leaves a lot of noise around that peak in corrected signal.
Let’s check Power Spectrum of first channel
Pr1 = sp.Periodogram(Xf[:,0])
Pr2 = sp.Periodogram(XR[:,0])
frq = (fs/2)*np.arange(len(Pr1))/(len(Pr1)-1)plt.plot(frq, np.log10(Pr1),label='Xf')
plt.plot(frq, np.log10(Pr2),label='XR')
plt.legend()
plt.grid()
plt.xlim([0,frq[-1]])
plt.ylabel('Power Spectrum')
plt.xlabel('Frquency (Hz)')
plt.show()
As explained that ICA introduces extra frequency components, that can be seen here. Even around 40Hz, some frequency components are amplified and after 50Hz, more frequency components with as good magnitude as 45Hz. These high frequency components are associated with that noise left behind by ICA.
You could change the parameters of ICA Algorithm, as spkit allows you to control window size, ICA method ICA_method(can be FastICA, Infomax or Extended Infomax), threshold on Kurtosis kur_thr(2 default) and correlation threshold corr_thr(0.8 default) Check help(sp.eeg.ICA_filtering)
ICA_method='extended-infomax', kur_thr=2, corr_thr=0.8ATAR Algorithm
Same as in ICA, following code apply ATAR algorithm on each window (128 samples = 1sec) with 50% overlapping to reconstruct 16 sec signals with overlap add method.
Parameters are using db4 wavelet, β = 0.1, threhold method = ‘ipr’, Operating Mode = ‘soft’. Check ref [1] for more details on parameters or check help(sp.eeg.ATAR)
XR1 = sp.eeg.ATAR(Xf.copy(),wv='db4', winsize=128, beta=0.1,thr_method='ipr',OptMode='soft', verbose=1)WPD Artifact Removal
WPD: True Wavelet: db4 , Method: ipr , OptMode: soft
IPR= [25, 75] , Beta: 0.1 , [k1,k2]= [10, 100]
Reconstruction Method: custom , Window: [‘hamming’, True] , (Win,Overlap)= (128, 64)
Let’s plot and see
plt.figure(figsize=(15,10))
plt.subplot(221)
plt.plot(t,Xf+np.arange(-7,7)*200)
plt.xlim([t[0],t[-1]])
#plt.xlabel('time (sec)')
plt.yticks(np.arange(-7,7)*200,ch_names)
plt.grid()
plt.title('X: Filtered signal',fontsize=16)
plt.subplot(222)
plt.plot(t,XR1+np.arange(-7,7)*200)
plt.xlim([t[0],t[-1]])
#plt.xlabel('time (sec)')
plt.yticks(np.arange(-7,7)*200,ch_names)
plt.grid()
plt.title('XR: Corrected Signal',fontsize=16)
#plt.show()
#plt.figure(figsize=(12,5))
plt.subplot(223)
plt.plot(t,(Xf-XR1)+np.arange(-7,7)*200)
plt.xlim([t[0],t[-1]])
plt.xlabel('time (s)')
plt.yticks(np.arange(-7,7)*200,ch_names)
plt.grid()
plt.title('Xf - XR: Difference \n(removed signal)',fontsize=16)
# #plt.show()
plt.subplot(224)
plt.plot(t,Xf[:,0],label='Xf')
plt.plot(t,XR1[:,0],label='XR')
#plt.plot(t,Xf[:,0]-XR[:,0])
plt.xlim([t[0],t[-1]])
plt.xlabel('time (s)')
plt.grid()
plt.legend()
plt.tight_layout()
plt.show()
Much better!!!
First thing that can be noticed is removed part from signal (bottom-left figure), a lot of part of removed signal is zero, which alone shows that ATAR is not removing more than necessary part of signal, unlike ICA. Second thing to notice is left out part in EEG around the peaky part (at 10 sec). Unlike ICA, ATAR do not leave extra noise around it, which can further be verified with periodogram.
Let’s see power spectrum, using periodogram.
Pr1 = sp.Periodogram(Xf[:,0])
Pr2 = sp.Periodogram(XR1[:,0])
frq = (fs/2)*np.arange(len(Pr1))/(len(Pr1)-1)plt.plot(frq, np.log10(Pr1),label='Xf')
plt.plot(frq, np.log10(Pr2),label='XR')
plt.legend()
plt.grid()
plt.xlim([0,frq[-1]])
plt.ylabel('Power Spectrum')
plt.xlabel('Frquency (Hz)')
plt.show()
Again, comparing to ICA, it is clearly a better one.
Let’s Try with Elimination Mode OptMode='elim' but less aggressive beta=0.01
XR2 = sp.eeg.ATAR(Xf.copy(),wv='db4', winsize=128, beta=0.01,thr_method='ipr',OptMode='elim')# Plotsplt.figure(figsize=(15,10))
plt.subplot(221)
plt.plot(t,Xf+np.arange(-7,7)*200)
plt.xlim([t[0],t[-1]])
#plt.xlabel('time (sec)')
plt.yticks(np.arange(-7,7)*200,ch_names)
plt.grid()
plt.title('X: Filtered signal',fontsize=16)
plt.subplot(222)
plt.plot(t,XR2+np.arange(-7,7)*200)
plt.xlim([t[0],t[-1]])
#plt.xlabel('time (sec)')
plt.yticks(np.arange(-7,7)*200,ch_names)
plt.grid()
plt.title('XR: Corrected Signal',fontsize=16)
#plt.show()
#plt.figure(figsize=(12,5))
plt.subplot(223)
plt.plot(t,(Xf-XR2)+np.arange(-7,7)*200)
plt.xlim([t[0],t[-1]])
plt.xlabel('time (s)')
plt.yticks(np.arange(-7,7)*200,ch_names)
plt.grid()
plt.title('Xf - XR: Difference \n(removed signal)',fontsize=16)
# #plt.show()
plt.subplot(224)
plt.plot(t,Xf[:,0],label='Xf')
plt.plot(t,XR2[:,0],label='XR')
#plt.plot(t,Xf[:,0]-XR[:,0])
plt.xlim([t[0],t[-1]])
plt.xlabel('time (s)')
plt.grid()
plt.legend()
plt.tight_layout()
plt.show()
Looks better, You could almost see all the other activities in Corrected EEG with a hint of absence of peak at 10 sec.
As name suggest, ATAR is a tunable algorithm, you might want to tune parameters like β with different mode of operation and wavelet which helps your prediction modeling.
Comparison & Conclusion
Artifacts in EEG are most disruptive for the prediction modelling. While working on project including EEG signal, removing artifacts is the first job. While there are many algorithms available (ICA and wavelet based), most of them are very aggressive to remove presumed artifactual components. ATAR algorithm helps to control the aggressiveness of algorithm and provide tuning parameters such as β and different mode of operation. ATAR is based on Wavelet Packet Decomposition and can be applying to single channel EEG.
As seen above figures, it can be clearly seen that ATAR performs better than ICA. Comparing time domain signals alone, it is very clear that ICA is very aggressive and induce extra frequency components. Comparing residual of both algorithms in figure below, it can be seen that ATAR is less aggressive.
Comparing Power Spectrum of first channel with both algorithms, it is further clear that no extra frequency components are produced by ATAR.
Tunable ability of ATAR makes it unique for predictive modeling.
Which is an important aspect, since with EEG, we are never sure if any component removed had any useful neural information for predictive task or not.
Here is a comparison Table to summaries the differences
Code Download
- Download entire code of this blog from here
- For further on tuning, check Jupyter-notebook here
References:
[1] Bajaj, Nikesh, et al. “Automatic and tunable algorithm for EEG artifact removal using wavelet decomposition with applications in predictive modeling during auditory tasks.” Biomedical Signal Processing and Control 55 (2020): 101624.
[2] Bajaj, Nikesh, Requena Carrión, Jesús “SpKit: Signal Processing Toolkit”, Python Package, https://spkit.github.io, version=0.0.9.4
