An Evaluation of EEG Ocular Artifact Removal with a Multi-channel Wiener
Filter Based on Probabilistic Generative Model
Hayato Maki1 , Tomoki Toda2 , Sakriani Sakti2 , Graham Neubig2 , and Satoshi Nakamura2

Abstract— Data contamination by ocular artifacts such as eye blinks
and eye movements is a major barrier that must be overcome whenattempting to analyze electroencephalogram (EEG) and event-related
potential (ERP) data. To handle this problem, a number of artifact
removal methods has been proposed. Specifically, we focus on a method
using a multi-channel Wiener filters based on a probabilistic generativemodel. This method assumes that the observed signal is the sumof multiple signals elicited by psychological or physical events, andseparates the observed signal into each event signal using estimatedmodel parameters. Based on this scheme, we have proposed a model
parameter estimation method using prior information of each event
signal. In this paper, we examine the potential of this model to deal
with highly contaminated signals by collecting EEG data intentionallycontaminated by eye blinks and relatively clean ERP data, and using
them as prior information of each event signal. We conducted anexperimental evaluation using a classical attention task. The results
showed the proposed method effectively enhances the target ERP
component while reducing the contamination caused by eye blinks.

I. INTRODUCTION
Due to its high temporal resolution, electroencephalogram (EEG)
has been one of the major techniques to investigate human brain
dynamics [1]. However, it is very sensitive to external artifacts suchas eye blinks, body movements, and so forth, which is a majorcommon problem to all EEG studies.
In an experimental paradigm using event-related potential (ERP),
the whole recording is cut into small intervals containing a single
stimulus. Each of interval is called a trial. Trial signals are usually
averaged synchronously to attenuate the effect of background EEGs
while preserving the amplitude of ERPs.
However, this synchronous averaging is not effective for removing external artifacts, because artifacts have far larger potentials,
several tens or hundreds times as large as neural sources. Consequently, external artifacts distort the averaged signal because thenumber of available trials is rarely enough to drown out the effect of
artifacts. Therefore, in order to avoid such distortion, it is common
usually discard trials in which the potential exceeds some predefined
threshold either in specific electrodes or in all electrodes, thenaverage remaining trials. Nevertheless, this procedure is problematicor even unacceptable when, as is often with the case with ERP
studies, the number of available trials is small or artifacts happen
frequently, as the number of remaining trials is too small to performreliable analysis. In addition, this procedure makes it impossible
to conduct an experimental paradigm in which the occurrence ofartifacts is inevitable such as visual tracking experiments. Especially, ocular artifacts are the most troublesome among various
artifacts since they are likely to happen frequently, especially duringvisual experiments. In fact, Small [2] reported that almost all trials*This work was not supported by any organization.
*The Ethical Review Board of Nara Institute of Science and Technology(NAIST) approved all experimental procedures involving human subjects.
1 Hayato Maki is with Graduate School of Information Science, NAIST,Japan maki.hayato.lt3@is.naist.jp
2 Tomoki Toda, Sakriani Sakti, Graham Neubig and Satoshi Nakamura
are Faculty of Graduate School of Information Science, NAIST, Japan

were contaminated by ocular artifacts in his visual ERP experimentconducted on children with autistic disorders.
For the reasons above, methods to remove the effect of artifacts
from recorded EEG data are necessary. A recent successful approach to the problem of noise reduction of EEG signals is independent component analysis (ICA), which decomposes a multi-channel
signal in a set of sources with maximally independent components(ICs). However, it has limitation on the number of separable ICs, N
ICs from N electrodes, which makes the decomposition imperfectas the number of EEG sources is much higher than the number ofICs [3], [4], [5]. In addition, ICA-based methods require subjective
decision making [6] or arbitrary tuning of the thresholds [7] todistinguish artifacted ICs from non-artifacted ICs, which is known
as the permutation problem.
Another approach [8], [9] using a multi-channel Wiener filter designed based on a probabilistic generative model has been proposed
to address the limitation of ICA. This approach refers to the multichannel signal related to the k-th event simply as the “k-th event
signal” and separates an observed signal not into signals generated
by individual sources but into event signals. This method has beensuccessfully applied to unsupervised EEG event signal separation
using multi-channel EEG signals [10]. Furthermore, in our previous
paper [11], we extended this method by using prior information of
the target event signal, allowing us to specifically enhance the target
event signal. This approach was shown to improve enhancement
performance and effectively remove background EEGs from singletrial ERPs. Moreover, in contrast to BSS techniques including
ICA, this method has no problems of permutation, which is often
troublesome in EEG/ERP research. However, the effectiveness of
our approach was verified only for removing background EEGs andit is still not clear it can remove more severe noise such as ocular
artifacts.
In this paper, we extend this work to ocular artifact removal from
ERP data by considering prior information of event signals elicitedby eye blinks. This procedure is essential because the frequency
spectra of potentials given by ocular artifacts and ERP components
are similar. For an experimental evaluation, we use ERPs intentionally contaminated by eye blinks and show the proposed method caneffectively remove the effect of the ocular artifacts and enhance thetarget event signal of ERPs.
II. P ROPOSED M ETHOD
A. Observation Model
Given I channels, the k-th event signal ck (n, f ) =[]⊤ck,1 (n, f ), · · · , ck,I (n, f )at time n and frequency f in thetime-frequency domain is expressed asck (n, f )=∑ hl sl (n, f ),(1)

[h1l , · · · , hIl ]⊤ ,

(2)

l∈Ek

hl

=

where {·}⊤ is the transpose, sl (n, f ) is the l-th source signal fromvarious sources, such as synapse, muscle, and so on, given by acomplex value, Ek is a set of the sources activated in the k-th event,and hil is the transfer function from the l-th source to the i-thchannel assuming that 0 ≤ hil ≤ 1. The observed multi-channel EEG
signal is x(n, f ) = [x1 (n, f ), · · · , xI (n, f )]⊤ expressed asK

x(n, f ) = ∑ ck (n, f ),

(3)

k=1

where K is the number of events.
B. Probabilistic Generative Model
We assume that the probability density function of the source
signal sl (n, f ) is modeled by the following zero-mean complex
Gaussian distribution,p (sl (n, f )) = Nc (sl (n, f ); 0, vk (n, f )) ,l ∈ Ek ,(4)

because EEG signals always contain multiple event signals, we candesign an experiment that elicits a specific ERP component, andwe can record EEG signals intentionally contaminated by a specificartifact e.g., eye blinks. Thus we can record EEG signals in which atarget event signal is supposed to be observed and also record EEG
signals which contain a non-target signal to establish a contrast
between them.
Based on these recorded signals, we calculate the hyper parameters iteratively as follows:
1: {Initialization}vk (n, f ) ← 1

2: repeat3:{Iterative update}(where the variance vk (n, f ) varies in the time-frequency domain depending on the k-th event. Further assuming that the source signals
are non-correlated to each other, the probability density function of
the k-th event signal ck (n, f ) is modeled by a multivariate complex
Gaussian distribution as follows:(5)=

Nc (ck (n, f ); 0, Rck (n, f )) ,[]E ck (n, f )ck (n, f )H=vk (n, f )Rk ,(7)∑ hl h⊤l ,(8)p (ck (n, f ))=Rck (n, f )Rk=where {·}H is the complex conjugate transpose. The spatial covariance matrix Rck (n, f ) is factorized into the time-frequency invariant
spatial covariance matrix Rk and the time-frequency variant variance component vk (n, f ).
We also assume that only one event signal is active in each timefrequency slot as follows:(9)

where z(n, f ) is the index of the active event signal. Consequently,
the probability density function of the observation signals x in thetime-frequency domain is modeled by a Gaussian mixture model
as follows:p(x|θ )=∏ p(x(n, f )|θ )
n, f

K

=

∏ ∑ αk Nc (x(n, f ); 0, vk (n, f )Rk ) ,

)
/NF

(
)
v̂k (n, f ) ← c̄k (n, f )H Ψ−1k c̄k (n, f ) /I
4: until convergence
5: Ψk ← NFΨk

(6)

l∈Ek

x(n, f ) = cz(n, f ) (n, f ),Ψk ←c̄k (n, f )c̄H (n, f )∑ v̄k (n,k f )
n, f

(10)

n, f k=1

where αk is a prior probability that the k-th event signal is active.
The model parameter set θ consists of αk , vk (n, f ), and Rk of each
mixture component.
C. Spatial Correlation Prior
The Wishart distribution is known as the conjugate prior distribution of the precision matrix of a multivariate Gaussian distribution
with known mean vectors. The prior distributions of time-frequency
invariant spatial covariance matrices are designed as follows:
(
(
) 1
])
1 [
−1 −1
−1 q−I−1−12
p Rk |Ψk , q = |Rk |exp − Tr Ψk Rk, (11)
Z
2
where Z is the normalizing constant, Ψk is a I-by-I symmetricpositive definite matrix, and qk is the degrees of freedom.
How to learn hyper parameters
Let us denote a pre-recorded signal that contains the k-th eventat time t by ĉk (t). While we can’t observe single event signal itself

where ĉk (n, f ) is the pre-recorded signal that contains the k-th
event at time-frequency slot (n, f ), N is the total number of timeframes, F is that of frequency bins, and v̄k (n, f ) is an estimationof the variance component vk (n, f ). The other hyper parameter qk
is calculated as qk = NF.
D. EEG Signal Enhancement with Spatial Correlation Priors
Given the observation signals x, the model parameter set θ isestimated by maximizing the posterior probability density functionof the time-frequency invariant spatial covariance matrices, R ={R1 , · · · , RK }, as follows:θ̂=

arg max p(R|x, θ\R , Ψk , q)

(12)

=

(
)
K
arg max ∏ p (x(n, f )|θ ) ∏ p R−1
k |Ψk , q ,

(13)

θ
θ

n, f

k=1

where θ\R is the model parameter set except for R. This maximization process can also be effectively solved with EM algorithm. The
following auxiliary function, which is the conditional expectationof the complete-data log likelihood evaluated for the posterior
probability in the equation (13) and some parameter set θ is
maximized with respect to θ ,(Q = ∑ mk (n, f ) log(αk ) − I log(vk (n, f ))n, f ,k)1
H −1+ log(|R−1k |) − v (n, f ) x(n, f ) Rk x(n, f )
k
])
K (
1 [
qk − I − 1−1
log |R−1
|
−
Tr
Ψ
R
+
k k
k
2
2
k=1

∑

+ Const.

(14)

In the E-step, the posterior probability mk (n, f ) is calculated at eachtime-frequency slot as follows:
mk (n, f ) =αk Nc (x(n, f ); 0, vk (n, f )Rk )∑K
k′ =1 αk′ Nc (x(n, f ); 0, vk′ (n, f )Rk′ )

.

(15)

=

v̂k (n, f )

=

R̂k

=

(16)
(17)

(18)

ĉk (n, f )=R̂ck (n, f )R̂−1
x (n, f )x(n, f ),(19)

R̂ck (n, f )=mk (n, f )v̂k (n, f )R̂k ,

(20)

K

=

∑ R̂c (n, f ),
k

10
10
0
0

(b)

0
0

200
200

400
400

600
600

100
100

where v̂k (n, f ) and R̂k are iteratively updated as they dependon each other. Finally, the target event signal is extracted fromthe observed EEG signals using a multi-channel Wiener filter asfollows:R̂x (n, f )

20
20

−10
−10

Amplitude (µV)α̂k∑n, f mk (n, f ),∑n, f ,k′ mk′ (n, f )1x(n, f )H R−1
k x(n, f ),
I
(
1
1
Ψ
qk − I − 1
2 k+ ∑ mk (n, f )2n, f)mk (n, f )H+∑x(n, f )x(n, f ),n, f v̂k (n, f )(a)Amplitude (µV)

In the M-step, model parameters α̂k , v̂k (n, f ), Rk are updated asfollows:

50
50
0
0
−50
−50

2000
2000

(21)

III. E XPERIMENTAL E VALUATION
While we showed in [11] that the previously described method iseffective in separating ERPs from background EEGs, it is not clear
whether these results generalize to ocular artifacts. This section
reports the result of an experiment in which we apply the proposed
method to an ERP data set contaminated by eye blinks to evaluate
its performance.
All EEG signals were recorded from 25 scalp electrodes at
locations based on a modified International 10-20 system, digitally
sampled at 1000 Hz, and downsampled to 200 Hz. A subject, who
was a right-handed 22 years old woman without any neurologicaldisorders participated in all experiments.
A. Prior Information Acquisition
In order to collect prior information of event signals that corresponded to P3 and eye blinks respectively, we conducted thefollowing EEG data recording.
First, we conducted an oddball paradigm experiment to obtain atemplate of P3. A random sequence of auditory stimuli including
2 kHz and 1 kHz sine waves was presented to the subject. 1 kHz
sounds were presented 200 times as Non-target stimuli and 2 kHzsine wave were presented 50 times as Target stimuli. The subject
was told to count the number of Target stimuli without any body
movement. Single-trial ERPs of each condition were averaged andthe result is shown in Fig. 1 (a). As expected, Target stimuli elicitedlarger P3 components. We used the averaged signal elicited byTarget stimuli as a template of P3 wave and used it to calculate
hyper parameters Ψ1 and q1 . Similarly, Ψ2 and q2 were calculatedusing the averaged signal elicited by Non-target stimuli.

6000
6000

8000
8000

Time (ms)k=1where ĉk (n, f ) is the k-th event signal separated from x(n, f ).
In the proposed enhancement method, we don’t have to selectthe target event signal from the separated event signals as in theconventional blind source separation methods, because the target
signal corresponds to the mixture component with the target priordistribution. We may also deal with event signals not modeled withthe prior distributions by just using additional mixture components
without the prior distributions.

4000
4000

Fig. 1. (a) The result of the oddball experiment. The red line is for the
averaged signal Target stimulus and the blue line is for Non-target stimulus.
(b) An EEG signal recorded while the subject was making eye blinks every
1.5 second. Peaks that have large potential are effects of eye blinks. All
plotted signals were recorded at channel Cz.

Next, we recorded an EEG signal for 2 minutes while the subject
was making an eye blink every 1.5 seconds, which was used to
calculate hyper parameters Ψ3 and q3 and shown in Fig. 1 (b). The
amplitude peaks is considered to be generated by the eye blinks.
B. Test Data Acquisition
To collect the data for evaluation, we used a modified oddballparadigm where the subject was told to make an eye blink every
time she heard the either type of stimulus and count the number ofpresence of Target stimuli only.
C. Result
We averaged the Target trial signal and Non-target signal of thetest data respectively. As expected, both of the averaged signals haveextremely large potentials (over 100 µ V) and P3 component wereburied under the effect of eye blinks, which makes the potential of
the Non-target averaged signal larger than Target averaged signal(Fig. 2 (a)). We applied the proposed method to the test data set
to remove the effect of the eye blinks setting the total number ofthe events K = 3, and averaged the signals for each condition. The
result is shown in Fig. 2 (b). The proposed method removed thelarge potentials made by eye blinks in both conditions, and madethe difference between them clear.
Furthermore, we show scalp maps of the test data before andafter the proposed method was applied (Fig. 3). Many studies havereported that P3 wave activates about at the central area of thescalp [12]. Before the proposed method is applied, we can see test
data is obviously contaminated by eye blinks as the area near theeyes is more highly activated instead of the central area. After its
application, we can see the activation of the central area, and seethat the potentials elicited by eye blinks have been removed.

Amplitude (µV)

(a)

(a)

L

100

3

50

1

0

0

0

200

400

600

800

1000

]

(b)

]

high

0
−5

0

200

400

600

800

−99.4

−49.7

0

5

49.7

10
99.4

Amplitude (µV)

15

low

1000

(b)

Time (ms)Fig. 2. The test data given by the modified oddball paradigm. (a) The redsolid line is for the averaged signal Target-stimuli and the blue dot line isfor Non-target stimuli. (b) The test data processed by the proposed method.

12

6

IV. CONCLUSIONSACKNOWLEDGMENT

ROW
ms

from

−6

−12

500
LatencyIn this paper, we addressed the problem of ocular artifact removalfrom EEG data. We assumed a probabilistic generative modelfor each event signal and estimated model parameters using priorinformation obtained from pre-recorded signals, which enabled us
to enhance only the target signal without the permutation problem.
Experimental evaluation showed that the proposed method caneffectively remove the effect of eye blinks in recorded ERP data.

0

Student Version of MATLAB

Fig. 3. The scalp maps given by (a) test data and (b) test data processedby the proposed method. Both of them were drawn at the peak latency of
their own averaged signals elicited by Target-stimuli.

Part of this work was supported by JSPS KAKENHI Grant
Number 26540117.
R EFERENCES[1] Ernst Niedermeyer and FH Lopes da Silva. Electroencephalography:
basic principles, clinical applications, and related fields. 2005.
[2] Joyce G Small. Sensory evoked responses of autistic children. Infantile
Autism, 1971.
[3] Sylvain Baillet, John C Mosher, and Richard M Leahy. Electromagnetic brain mapping. IEEE Signal Processing Magazine, 2001.
[4] Arthur K Liu, Anders M Dale, and John W Belliveau. Monte Carlosimulation studies of EEG and MEG localization accuracy. Humanbrain mapping, 2002.
[5] David M Groppe, Scott Makeig, and Marta Kutas. Identifying reliableindependent components via split-half comparisons. NeuroImage,
2009.
[6] Tzyy-Ping Jung, Scott Makeig, Marissa Westerfield, Jeanne Townsend,
Eric Courchesne, and Terrence J Sejnowski. Removal of eye activity
artifacts from visual event-related potentials in normal and clinicalsubjects. Clinical Neurophysiology, 2000.
[7] Carrie A Joyce, Irina F Gorodnitsky, and Marta Kutas. Automatic
removal of eye movement and blink artifacts from EEG data using
blind component separation. Psychophysiology, 2004.
[8] Ngoc QK Duong, Emmanuel Vincent, and Rémi Gribonval. Underdetermined reverberant audio source separation using a full-rank
spatial covariance model. IEEE Transactions on Audio, Speech, andLanguage Processing, 2010.

[9] Ryutaro Sakanashi, Shigeki Miyabe, Takeshi Yamada, and Shoji
Makino. Comparison of superimposition and sparse models in blindsource separation by multichannel Wiener filter. Proceedings of
the Annual Summit and Conference of the Asia-Pacific Signal andInformation Processing Association (APSIPA ASC), 2012.
[10] Yusuke Kurihana, Shigeki Miyabe, Tomasz M. Rutkowski, Yoshihiro
Matsumoto, Takeshi Yamada, and Shoji Makino. (in Japanese) Signal
Separation of EEG using Multivariate Probabilistic Model. Technical
Report of IEICE, MBE, 2013.
[11] Hayato Maki, Tomoki Toda, Sakriani Sakti, Graham Neubig, andSatoshi Nakamura. EEG Signal Enhancement using Multichannel
Wiener Filter with a Spatial Correlation Prior. Proceedings of the IEEE
International Conference on Acoustics, Speech, and Signal Processing(ICASSP), 2015.
[12] David Friedman, Yael M Cycowicz, and Helen Gaeta. The novelty P3:
an event-related brain potential (ERP) sign of the brain’s evaluationof novelty. Neuroscience & Biobehavioral Reviews, 2001.