IEICE TRANS. INF. & SYST., VOL.E99–D, NO.6 JUNE 2016

1437

PAPER

Special Section on Human Cognition and Behavioral Science and Technology

Enhancing Event-Related Potentials Based on Maximum a
Posteriori Estimation with a Spatial Correlation Prior
Hayato MAKI†a) , Nonmember, Tomoki TODA†,†† , Sakriani SAKTI† , Members,
Graham NEUBIG† , Nonmember, and Satoshi NAKAMURA† , Member

SUMMARY
In this paper a new method for noise removal from singletrial event-related potentials recorded with a multi-channel electroencephalogram is addressed. An observed signal is separated into multiple
signals with a multi-channel Wiener filter whose coeﬃcients are estimatedbased on parameter estimation of a probabilistic generative model that locally models the amplitude of each separated signal in the time-frequency
domain. Eﬀectiveness of using prior information about covariance matrices
to estimate model parameters and frequency dependent covariance matrices
were shown through an experiment with a simulated event-related potential
data set.
key words: electroencephalogram (EEG), event-related potential (ERP),
generative model, independent component analysis (ICA), Wiener filter,
noise removal, Wishart distribution, spatial correlation prior

1.

Introduction

Electroencephalogram (EEG) is one of the useful tools
to investigate human brain activity non-invasively, alongwith functional magnetic resonance imaging (fMRI), nearinfra-red spectroscopy (NIRS), and magnetoencephalography (MEG). Compared to other methods, EEG has advantages such as relatively low cost, high temporal resolution,and experimental flexibility. Event-related potential (ERP)
is a brain response to an internal or external event such as
visual stimuli or cognitive processes. ERP is a useful toolfor cognitive neuroscience, and often used for manipulationof a brain-machine interface (BMI). However, analysis of
ERP data usually suﬀers from its low signal-to-noise ratio(SNR) due to superimposition of task-unrelated brain activities as well as non-neural artifacts such as eye blinks, muscle
movements, and heart beats.
Signal averaging in the time domain is often employed
in order to improve SNR of ERP by attenuating taskunrelated background EEGs. In a usual ERP data analysis,
the whole EEG recording during an experiment is cut into
small intervals containing a single stimulus. Each of theseintervals is called a trial or an epoch. ERP components are
assumed to be identical over each trial, while background
EEGs are assumed to be arbitrary phase with regards to thestimuli. Therefore, taking the average of all trials aligned toManuscript received October 31, 2015.
Manuscript revised March 11, 2016.
Manuscript publicized April 1, 2016.
†
The authors are with Nara Institute of Science and Technology, Ikoma-shi, 630–0101 Japan.
††
The author is with Nagoya University, Nagoya-shi, 464–8601
Japan.
a) E-mail: maki.hayato.lt3@is.naist.jp
DOI: 10.1587/transinf.2015CBP0008

the onsets of stimuli of each trial will weaken the amplitudeof background EEGs while preserving that of ERPs.
However, variabilities of ERP across trials exist in latency, amplitude, and scalp distribution [1]. This acrosstrial averaging procedure sacrifices all the information concerning across-trial variability of brain response. Moreover,to perform reliable analysis at least 20 trials for each experimental condition are required to be averaged, resulting in long experimental time and subject fatigue. Therefore, single-trial analysis is an important research topic tostudy event-related brain dynamics more precisely [2]–[4].
Throughout this paper, we refer to the ERPs as the signal topreserve and task-unrelated neural activities as well as nonneural artifacts as noise to be removed.
With regards to external artifact removal, especially
eye blink removal, there is a large amount of previous research, and several eﬀective methods [5]. However, background EEG removal is less studied, although there is someresearch proposing methods based on independent component analysis (ICA) [6], [7].
In this paper, a new approach to remove background
activities from single-trial ERPs is addressed. The proposed
method is an extension of a framework that has emerged recently in under-determined sound source separation [8], [9].
While ICA aims to separate an observed signal into its
sources, this framework separates into signals elicited byeach event such as ERP or eye blinks. We denote the groupof source indexes that contributes to the k-th event by Ek anddenote the EEG signal that is elicited by Ek and captured
by ck (t) and k-th event signal. In the time-frequency domain, short-time Fourier transform (STFT) coeﬃcients ofthe event signal ck (n, f ) are locally modeled by a multivariate complex Gaussian distribution whose parameters are afunction of (n, f ), where n is the index of the time frame
and f is the index of the frequency bin. One of the parameters is a full rank time-invariant covariance matrix thatis called a spatial correlation matrix because it encodes
the spatial spread of the event signal. The model parameters are estimated by the maximum likelihood criterion andused to separate an observed multi-channel signal into someevent signal by a multi-channel Wiener filter. This framework was applied to EEG signal separation where it wasused to remove external ocular artifacts from alpha bandrhythmic activity [10] and background activities from an

c 2016 The Institute of Electronics, Information and Communication Engineers
Copyright IEICE TRANS. INF. & SYST., VOL.E99–D, NO.6 JUNE 2016

1438

ERP. We denote the group of source indexes corresponding
to the ERP by E. Then, signal reconstruction is done asfollows:⎧L⎪⎪⎨0 (i  E)
c = Ãŝ =ãi ŝi (t), ãi = ⎪(3)⎪⎩âi (i ∈ E)
i=1

Fig. 1

Overview of the proposed method.

ERP signal [11]† .
This approach is extended using prior information.
This paper provides two main contributions. First, the effectiveness of using prior information to estimate covariance matrix parameters is demonstrated. Second, the eﬀectiveness of frequency dependent spatial correlation matrices
are examined, in contrast to [10], which uses frequency independent ones. The overview of our proposed method isshown in Fig. 1.
2.

Existing Techniques for Signal Enhancement

In this section, we describe two diﬀerent types of signal
separation techniques, independent component analysis inSect. 2.1 and the maximum likelihood approach in Sect. 2.2.
2.1

Independent Component AnalysisM

ai si (t),Lλ (xk ) : xk −→ (1 − λ)xk + λx̄, λ ∈ [0, 1](4)where xk is the k-th trial, x̄ is the averaged signal across trials and λ is the regularization parameter. They applied the
ICA algorithm to Lλ (xk ) instead of xk to estimate the mixing matrix, then decomposed each xk using the inverse of
the estimated matrix. Another problem of ICA is permutation ambiguity. The ICA algorithm doesn’t tell us which
ICs correspond to the signal of interest, which decreases thepractical usability of ICA. The component classification isrequired to reconstruct a noise free signal that can be performed by visual inspection [3], [7], predefined [13], [14] or
adaptive [15] thresholding.
2.2 Maximum Likelihood Approach

ICA is the most commonly used technique in the context of
EEG signal separation. ICA generally models an observed
EEG signal x(t) at time t as an instantaneous linear combination of M sources that is captured by L sensors as follows:
x(t) = As(t) =In general, it is assumed that the number of sensors isequal to that of sources because the linear mixture is invertible in such a case. However, there is no reason to believe thenumber of EEG sources is equal to that of the channels of the
EEG recording equipment. In fact, EEG source separation
is often treated as an under-determined problem, implying
the condition M > L [7], [12]. In such a case, additional assumptions or prior knowledge about the underlying sources
are usually exploited. For example, in [7], they used the following function Lλ (x) to incorporate the prior information
about an averaged ERP:

(1)This approach has been proposed originally in the field ofsound source separation [9] but its eﬀectiveness for EEG
signals has not been explored extensively. This approach
has a virtue that, in contrast to ICA, it doesn’t assume thenumber of sources because it doesn’t estimate sources but
contributions of each event signal to each EEG electrode.

i=1

where A is a L × M mixing matrix whose i-th column is
ai . ICA is an algorithm that estimates a mixing matrix A
given only observation x such that sources s are as staticallyindependent as possible. An estimation of sources s is obtained by multiplying the inverse of the mixing matrix withthe observed signal as follows:−1ŝ(t) = Â x(t).

(2)

We refer to estimated sources ŝ = ( ŝ1 , · · · , ŝL ) as independent components (ICs) and Â is an estimated mixing matrix
whose i-th column is âi . In order to obtain a noise-free ERP
signal c, it is necessary to select ICs corresponding to the†
In this paper, we present further details of the proposed methods, more discussions, and more evaluations than those in our previous work [11].

2.2.1

Observation Model

The STFT coeﬃcients of the k-th event signal in timefrequency slot (n, f ) are expressed as ck (n, f ):hl sl (n, f ),(5)
ck (n, f ) =l∈Ek

ck (n, f ) = [ck,1 (n, f ), · · · , ck,L (n, f )] ,hl = [h1l , · · · , hLl ] ,(6)

where hil is the transfer function from the l-th source to thei-th channel assuming that 0 ≤ hil ≤ 1. The observed multichannel EEG signal x(n, f ) is expressed asx(n, f ) = x1 (n, f ), · · · , xL (n, f )  =ck (n, f ),
K

k=1

(7)

MAKI et al.: ENHANCING ERPS BASED ON MAP ESTIMATION WITH A SPATIAL CORRELATION PRIOR

1439

where K is the total number of events that should be given
as a hyper-parameter.
Recall that in ICA’s observation model, observed signals are expressed as the sum of the product of time-variant
scalars and time-invariant vectors in (1). Here, instead theobserved signals are expressed as the sum of time-variant
vectors, which means that signal separation of event signalsthat take this form can not be achieved by ICA.
2.2.2

Generative Model

To model the generation of observed signals probabilistically, we make the following two assumptions:
Assumption1
The amplitude of the l-th source that contributes to thek-th event in each slot (n, f ) follows a complex normal
distribution with mean 0. The variance is given by the
product of the degree of the event activity vk (n, f ) andthe source activity λl :=vk (n, f )λl hl hll∈Ek= vk (n, f )Rk ,where {·}H is the complex conjugate transpose. We call Rk
the spatial correlation matrix of the k-th event signal.
2.2.3

Event SparsenessTo simplify the generative model, we assume that the eventsin each time frequency slot are sparse. In other words, we
assume that we observe only one event signal in each slot.
In order to express this sparseness mathematically, we introduce the latent variables zk (n, f ), which take a 1-of-K representation, in which one particular element zk (n, f ) is equal
to 1 and all other elements are equal to 0. In other words,
the value of zk (n, f ) satisfies
zk (n, f ) ∈ {0, 1},

p (sl (n, f )) = Nc (sl (n, f ); 0, vk (n, f )λl )=(8)where l ∈ Ek and the group Ek contains source indexes
associated with the k-th event signal.
Assumption2
Sources are not correlated with each other within eachtime-frequency slot (n, f ):E[sl1 (n, f )s∗l2 (n, f )] = 0 for l1  l2 ,exp −ck (n, f )πL det(Rck (n, f ))=l1 ∈Ek l2 ∈Ek

=

l∈Ek

E[|sl (n, f )|2 ]hl hlKzk (n, f )ck (n, f ) = cl (n, f ),

(13)

k=1

where zl (n, f ) = 1.
2.2.4

Observation Likelihood

From all assumptions above, finally we derive the likelihoodof the observed signals x in the time-frequency slot (n, f ),which is modeled by a Gaussian mixture model as follows:p(x|θ) =p(x(n, f )|θ)
n, f

.

(10)

For simplicity, we restrict the covariance matrix Rck (n, f )to the product of the time-invariant covariance matrix Rkthat encodes spatial spread of the k-th event signal andthe time-frequency variant scalar vk (n, f ) that encodes timefrequency power of the signal as follows:Rck (n, f ) = E ck (n, f )ck (n, f )H
⎤
⎡
H⎥

⎢⎢⎢ 
⎥⎥
hl1 sl1 (n, f )hl2 sl2 (n, f ) ⎥⎥⎥⎦
= E ⎢⎢⎢⎣
l1 ∈Ekl2 ∈Ek 
E[sl1 (n, f )s∗l2 (n, f )]hl1 hl2 ]
=


(12)

If we observe the l-th event signal in the slot (n, f ), we can
express the observed signal as:
x(n, f ) = x1 (n, f ), · · · , xL (n, f ) p (ck (n, f )) = Nc ck (n, f ); 0, Rck (n, f )=zk (n, f ) = 1.

(9)where l1 , l2 ∈ Ek and {·}∗ indicates the complex conjugate ofsl2 (n, f ).
From the assumptions above, the probability density
function that each event signal follows is derived, namely
the multivariate complex normal distribution with zero
mean.
R−1
ck (n, f )ck (n, f )

K

k=1

|sl (n, f )|21exp,
πvk (n, f )λlvk (n, f )λl

H

(11)

=

K


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

n, f k=1

p(zk (n, f ) = 1) = αk for all (n, f ),(15)

where αk is the prior probability that the k-th event signal
is observed assuming that their values are shared among allslots. The model parameter set θ consists of αk , vk (n, f ), and
Rk of each mixture component.
2.3 Event Separation
The optimal model parameter set θ is estimated by maximizing the observation likelihood described in Eq. (14). Thismaximization process can also be eﬀectively solved with the
expectation-maximization (EM) algorithm. θ is chosen by
maximizing the following auxiliary function, which is theconditional expectation of the complete-data log likelihood:

IEICE TRANS. INF. & SYST., VOL.E99–D, NO.6 JUNE 2016

1440

Q=



mk (n, f ) log(αk ) − L log(vk (n, f ))n, f,k+ log(|R−1k |) −

1x(n, f )H R−1
k x(n, f ) + Const,
vk (n, f )(16)where in the E-step, the posterior probability mk (n, f ) is calculated at each time-frequency slot:
αk Nc (x(n, f ); 0, vk (n, f )Rk ).
k =1 αk Nc (x(n, f ); 0, vk (n, f )Rk )mk (n, f ) = K

In the M-step, α̂k and v̂k (n, f ) are updated as follows:n, f mk (n, f )α̂k = ,n, f,k mk (n, f )1v̂k (n, f ) = x(n, f )H R−1
k x(n, f ),L mk (n, f )Hn, f v̂k (n, f ) x(n, f )x(n, f )

R̂k =
,
n, f mk (n, f )

(17)

Mk

(18)
(19)

(21)

(22)

where Rab = E[ab ] and Ra = E[aa ]. Assuming event
signals are independent from each other, Rck (n, f )x (n, f ) canbe expressed as:
Rck (n, f )x (n, f ) = [ck (n, f )x(n, f ) ]
= ck (n, f )ck (n, f )

3.

q−p−1
2

exp −

Tr(W−1 A)
2

where W is a p × p positive definite matrix and q is called
the number of degrees of freedom and is restricted to q ≥ p,
and B(W, p, q) is a normalizing factor defined byqpq

B(W, p, q) = |W|− 2 2 2 πp(p−1)
4

where Γ(·) is the Gamma function.
The Wishart distribution is the conjugate prior distribution of a multivariate Gaussian random variable’s covariancematrix. Assume that we observe b-dimensional random vectors that follow the normal distribution with the known mean
vector μ and the unknown precision matrix Λ, and want toinfer Λ. In such a case, the conjugate prior distribution forΛ is the Wishart distribution with a hyper parameter of b × b
matrix [16].
3.2 Maximum a Posteriori Estimation of Spatial Correlation Matrices
From the Bayes’ theorem, the posterior probability density
function of the time-frequency invariant spatial covariancematrices R = {R1 , · · · , RK } is proportional to the product of
the likelihood of the observed data and their prior as follows:

n, f

(23)

(24)

p
 q + 1 − i  −1

Γ
.
2
i=1

(26)

∝

Consequently, an event signal is estimated as follows:ĉk (n, f ) = mk (n, f )v̂k (n, f )R̂k R−1
x (n, f )x(n, f )W p (W, q) = B(W, p, q)|A|p(R|x,θ\R , Ψk , qk , wk )= Rck (n, f )= mk (n, f )v̂k (n, f )R̂k .

The Wishart distribution is a type of probability distribution,
whose probability density function is given by(20)

where ĉk (n, f ) is the k-th event signal separated from x(n, f ).
Solving this minimization problem, we obtainM̂k (n, f ) = Rck x (n, f )R−1
x (n, f ),3.1 Spatial Correlation Prior(25)

where v̂k (n, f ) and R̂k ( f ) are iteratively updated because
they depend on each other.
Finally, the target event signal is extracted from the observed EEG signals using a multi-channel Wiener filter. A
multi-channel Wiener filter M̂k ∈ RL×L is a linear filter thatminimizes the square errors as follows:M̂k = arg min Mk x(n, f ) − ĉk (n, f ) 2

described in the previous section. Hence, we introduce priorinformations of event signals to estimate the spatial correlation matrices utilizing the property of the Wishart distribution.

p (x(n, f )|θ)

(27)
K



wk
−1
WL R−1(28)
k |Ψk , qk

k=1

where x is the M-dimensional observation signal, θ\R is themodel parameter set except for R, and wk ≥ 0 is the weighting parameter that expresses our belief about the strength ofthe prior. If we have little confidence in the prior of Rk , weset a small value for wk and vice versa. When wk = 0, thisformula reduces to maximum likelihood (ML) estimation.

Proposed Method
3.3 Hyper Parameter Learning

We often know which event signal we would like to enhance, such as ERP, and we can record EEG signals related
to the target event beforehand and use them as prior knowledge for enhancement. In such a case, we need to enhance
the target event signal from the observed EEG signal, rather
than blindly separating it into multiple EEG event signals

We denote the j-th pre-recorded signal that contains the kth event at time t by c(kj) (t), although we can’t observe thesingle event signal itself because EEG signals always contain multiple event signals. However, we can design an experiment that elicits a specific ERP component, and we canMAKI et al.: ENHANCING ERPS BASED ON MAP ESTIMATION WITH A SPATIAL CORRELATION PRIOR

1441

record EEG signals intentionally contaminated by a specificartifact e.g., eye blinks. Thus we can record EEG signals
in which a specific event signal is supposed to be observed
and also record EEG signals that are not supposed to be observed to establish a contrast between them. In addition, we
can assume that we can obtain a noise-reduced event signal
by synchronous averaging if the event signal is time-lockede.g., ERPs.
The hyper parameters are calculated by following
steps:
1: {Initialization}v̄k (n, f ) ← 1(29)

2: repeat3:{Iterative update}Fig. 2 Scalp topographies within the diﬀerent frequency bins. Each line
stands for an EEG channel.

⎞⎛⎜⎜⎜ c̄k (n, f )c̄H (n, f ) ⎟⎟⎟k⎟⎟⎟ /NΨk ( f ) ← ⎜⎜⎜⎝
⎠
v̄
(n,
f
)
k
n, f


v̂k (n, f ) ← c̄k (n, f )H Ψ−1k c̄k (n, f ) /L

4: until convergence
where c̄k (n, f ) is a signal that contains the k-th event signal,
N is the total number of time frames, v̄k (n, f ) is an estimationof the scaling parameter vk (n, f ). If the k-th event signal istime-locked e.g., ERP, we use the averaged signal E j [ĉ(kj) (t)]as the template of the k-th event signal to calculate hyper
parameters where ĉ(kj) (t) is the j-th trial signal that contains
the k-th event signal:c̄k (t) = E j [ĉ(kj) (t)].

(30)

If the k-th event signal is not time-locked e.g., background
EEGs or eye blinks, we use the concatenated signal as follows:c̄k = [ĉ(1) (1) · · · ĉ(1) (t), · · · , ĉ( j) (1) · · · ĉ( j) (t)] .
3.4

(31)

Event Separation

The optimal model parameter set θ is estimated by maximizing the posterior probability described in Eq. (27) insteadof observation likelihood described in Eq. (14). This maximization process can also be eﬀectively solved with the EM
algorithm. The following auxiliary function is maximized
with respect to θ.
Q=n, f,k

+

K


wk

k=1

+ Const.

where v̂k (n, f ) and R̂k each other. Finally, the target event
signal is extracted from the observed EEG signals using a multi-channel Wiener filter constructed according toEqs. (21), (22), and (23).
3.5 Frequency Dependent Spatial Correlation Matrix
EEG signals usually have diﬀerent spatial spread over thescalp according to the frequency as shown in Fig. 2. Hence,
another generative model with frequency dependent spatial
correlation matrices Rk ( f ) and Ψk ( f ) instead of Rk and Ψkis proposed. In this manner, the estimation of spatial covariance matrices in the EM algorithm is done as follows: (n, f )HΨk ( f ) + 2 n mv̂kk(n,f ) x(n, f )x(n, f ).
(34)R̂k ( f ) =(qk ( f ) − L − 1) + 2 n mk (n, f )Note that the range of summation is over only time frame n.
4.

Experiments

To validate the performance of the proposed methods, weapplied the methods to pseudo-ERP data, which were madeby the superposition of measured background EEGs (noise)and the measured and averaged ERP of P300.

mk (n, f ) log(αk ) − L log(vk (n, f ))+ log(|R−1k |) −In the E-step, the posterior probability mk (n, f ) is calculated
at each time-frequency slot as shown in Eq. (18). In the Mstep, α̂k and v̂k (n, f ) are updated as shown in Eqs. (18) and(19). The spatial correlation matrices are updated as follows:(n, f )
Hwk Ψk + 2 n, f mv̂kk(n,f ) x(n, f )x(n, f )

.
(33)
R̂k =
wk (qk − L − 1) + 2 n, f mk (n, f )1x(n, f )H R−1
k x(n, f )vk (n, f )qk − L − 1
1−1
log |R−1
k | − Tr Ψk Rk
2
2

4.1 Subjects and EEG Recording(32)

EEG data were collected from three healthy male volunteers
aged from 23 to 26 years. All participants gave written informed consent and the Ethical Review Board of the Nara

IEICE TRANS. INF. & SYST., VOL.E99–D, NO.6 JUNE 2016

1442

Institute of Science and Technology approved all experimental procedures.
All EEG signals were recorded from 22 scalp electrodes at locations based on a modified International 10-20
system [17], digitally sampled at 1000 Hz, downsampled to
250 Hz referred to the average of the both sides of ears, andfiltered between 0.01 Hz and 30 Hz.
4.2

Experimental Paradigm

We conducted two sessions for each subject. In the first session, we conducted an auditory oddball paradigm experiment [18] using 1 kHz and 2 kHz sound. A random sequenceof auditory stimuli including 2 kHz and 1 kHz sine waves
was presented to the subject with an earphone. All the stimuli had the duration of 0.2 seconds and were presented atintervals of 1.4 seconds. 1 kHz sounds were presented 200
times as non-target stimuli and 2 kHz sine wave were presented 50 times as target stimuli. Subjects sat in the front of
a computer display and were told to look at the fixation markin the display and count the number of target stimuli without
any body movement during the session. At the onset of eachstimulus, the two small areas corresponding to two types ofstimuli in the corner of the display flashed. The photo sensor sent a trigger marker that marked the type of stimulus
presented and its onset time to the EEG recording computerwhen it caught the flash. The flashes were hidden from thesubjects using a thick paper to avoid distracting them.
In the second session, the subjects were told to relax
without any body movement to record resting state EEG for
two minutes.
4.3

Parameter Settings

We implemented regularized ICA based on [7]. We decomposed each raw single-trial multi-channel EEG signal
into two event signals with time-frequency modeling methods. Hyper-parameters of the prior distribution for the firstevent signal was calculated from the signal that were given
by averaging the validation dataset as described in Sect. 3.3.
Hyper-parameters for the second event signal were calculated from the resting state EEG signal recorded in the second session. All hyper-parameters are calculated individually for each subject.
4.4

Fig. 3 The target trial (left, solid blue line) and non-target trial (left,dashed red line) ERPs that are embedded into each of resting state signal(right), respectively.
Table 1
Conventional method
Proposed method 1
Proposed method 2

Prior
No
Yes
Yes

Frequency dependency
No
No
Yes

trial ERP was also added to resting state signals as sameas the target trial ERP to make non-target trial pseudo-ERP
data set. The averaged and embedded ERP for each type ofstimulus and an example of resting state trial were shown in
Fig. 3.
4.5 Evaluation
We compared five denoising methods, regularized ICA [7]and time-frequency modeling methods that have a spatial
correlation matrix with or without prior information and frequency dependency (See Table 1).
We evaluate these methods for each subject using amplitude deviation (AD), latency deviation (LD) and root
mean square errors (RMSE) defined as follows because
ERPs are often quantified by their peak amplitude and peaklatency,1|x̄(t¯peak ) − x̂( j) (t peak )|
J j=1
J

AD =

1
j)
|t¯peak − tˆ(peak
|
J j=1
!
"
J
T


1
RMSE =
x̄(t) − x̂(t)
J j=1
t

(35)

J

LD =

Creation of Validation Dataset

The multi-channel EEG data obtained from session 1 ofeach type of stimulus were cut into trials of 1.2 second
length, then averaged the trials across the trials for eachchannel, subject and type of stimulus. Next, the resting state
EEG obtained from session 2 was also cut into 40 signalsof 1.2 second length. The averaged target trial ERP of eachcorresponding channel that stands for signal was added tothe each of 40 resting state trial that stands for noise withrandomly shifted phase from −60 ms to 60 ms to maketarget-trial pseudo-ERP data set. The averaged non-target

Compared methods.

(36)

(37)

where x̄ is the averaged ERP signal obtained from session 1
added to the resting state EEGs and, x̂ is the j-th single-trial
denoised EEG signal, J is the total number of trials, T is thesignal length, and t¯peak is the time when the averaged ERPj)
is the time when
signal x̄ has the largest amplitude and t¯(peakthe j-th single-trial denoised EEG signal x̂( j) does.
5.

Results

Figure 4 shows the plots of a single-trial EEG signal of

MAKI et al.: ENHANCING ERPS BASED ON MAP ESTIMATION WITH A SPATIAL CORRELATION PRIOR

1443
Table 2

The result of the multiple comparison test (Target trials).

Method 1
Raw
Raw
Raw
Raw
ICA
ICA
ICA
Conv
Conv
Pro 1

Method 2
ICA
Conv
Pro 1
Pro 2
Conv
Pro 1
Pro 2
Pro 1
Pro 2
Pro 2

Table 3
Method 1
Raw
Raw
Raw
Raw
ICA
ICA
ICA
Conv
Conv
Pro 1

Fig. 4 Plots of a raw signal, a denoised signal, and the embedded averaged ERP: Comparison of regularized ICA (upper left), Winer Filter (Conventional Method, upper right), Winer Filter + Spatial Prior (Proposed
Method 1, lower left), and Winer Filter + Spatial Prior + Frequency Dependent please describe the (Proposed Method 2, lower right).

the validation dataset, a single-trial denoised signal, and theaveraged ERP signal that is embedded into the validation
dataset. The raw signal has two main peaks about 300ms
and 380ms. Compared to the averaged ERP, the former peakis eliminated by averaging while the latter peak still exists inthe averaged ERP. Hence, we can think that the former peakis made by chance due to background EEG to be removed by
the denoising methods while the latter peak is due to embedded ERP to be preserved. First, looking at the ICA results,
we can see that it removed the ERP component as well as
the background EEG. The conventional method, Wiener filtering without spatial correlation prior, on the other hand,preserved the noise as well as the ERP component. Both ofthe proposed method 1, which uses prior information andfrequency independent spatial correlation matrices, and the
proposed method 2, which uses prior information and frequency dependent spatial correlation matrices removed thenoise from the ERP. However, noise removal done by the
proposed method 1 was not enough to make the first peaksmaller than the second, resulting in worse estimation of thepeak latency. The proposed method 2 managed to remove

p-val (AD)**0.000
0.952**0.004**0.000**0.000**0.001*0.011*0.039**0.0040.960

p-val (LD)**0.000
1.000
0.449
1.000
**0.000
0.075**0.000
0.410
1.000
0.533

p-val (RMSE)
0.994
0.955
0.252
*0.0450.796
0.105*0.013
0.679
0.238
0.947

The result of the multiple comparison test (Non-target trials).
Method 2
ICA
Conv
Pro 1
Pro 2
Conv
Pro 1
Pro 2
Pro 1
Pro 2
Pro 2

p-val (AD)*0.0190.102**0.000
0.227
0.975
0.105
0.870
*0.0200.996**0.006

p-val (LD)**0.000
1.000*0.0470.594**0.000**0.000**0.000*0.0280.479
0.691

p-val (RMSE)**0.000
0.349**0.000
0.864**0.000
0.122**0.000*0.0310.911**0.002

the first peak by an amount suﬃcient to make the first peaksmaller than the second, allowing for better estimation of
the peak amplitude and latency, although it was not able toremove the first peak entirely.
In addition to this qualitative analysis, we also objectively measure the eﬀectiveness of the proposed method
quantitatively by measuring AD, LD, and RMSE, respectively. Figure 4 shows the averages and standard deviationsof these values over all trials for each type of stimulus. We
performed ANOVA to analyze diﬀerence among comparedmethods and we found significant diﬀerence between methods in either of AD, LD and RMSE of target and non-target
trials (p < 10−13 for AD, p < 10−4 for LD, and p < 0.01 for
RMSE of target trials, and p < 10−5 for AD, p < 10−21 for
LD and p < 10−12 for RMSE of non-target trials). We also
performed multiple comparison test with Tukey’s honestly
significant diﬀerence (HSD) criterion [19] and the result isshown in Table 2 for target trials and Table 3 for non-target
trials where Pro stands for the proposed method and Conv
stands for the conventional method and the asterisks * and** indicate statically significant diﬀerence at the 5 and 1 percent level, respectively. With regards to target trials, we cansee that regularized ICA resulted in significantly smaller AD
than raw signals and all the other denoising methods. However, it had also significantly larger LD than others. ICA
constantly outputted a denoised signal with small amplitude
regardless of how large amplitude the raw signal had. Bothof the proposed methods made AD significantly smaller than
the raw signal and the conventional method in target trials.
The proposed method 2 made RMSE significantly smaller
than raw signals and regularized ICA. With regards to nontarget trials, the proposed method 1 had significantly smaller
AD and RMSE compared to raw signals and the conven-

IEICE TRANS. INF. & SYST., VOL.E99–D, NO.6 JUNE 2016

1444

nent with the target prior distribution.
6.2 Frequency Dependency
We could not find the proposed method 2 (frequency dependent covariance model) always significantly better thanthe proposed method 1 (frequency independent covariancemodel) although the former could have been able to more
eﬀectively represent a frequency dependent spatial spread ofa multi-channel EEG signal as stated in Sect. 3.5. In the frequency dependent covariance model, the number of samples
x(n, f ) to estimate a covariance matrix is reduced to 1/F
compared to the frequency independent covariance model,
as the samples that belong to only one frequency bin are
involved in estimating the corresponding covariance matrix
as shown in Eq. (34) while all samples x(n, f ) in the timefrequency grid are involved in the frequency independent
model as shown in Eq. (33). This has an adverse eﬀect onrobust covariance matrix estimation and could deteriorate
signal separation performance placing a trade oﬀ between
frequency specific covariance matrix modeling and robust
covariance matrix estimation.
6.3 Target and Non-Target Trials

Fig. 5 Amplitude deviation of target trials (upper left), latency deviationof target trials (middle right), root mean square errors of target trials (lower
right), amplitude deviation of target trials (upper right), latency deviationof target trials (middle right), and root mean square errors of non-target
trials (lower right). The bar length are the means and the lines are standard
deviations of 40 trials for each type of stimuli.

tional method. The result of AD and RMSE was consistentin each of type of stimulus.
6.

Discussion

6.1

Eﬀectiveness of Using Prior Information

All the compared denoising methods are statistical methods that require suﬃciently sized data. However, we facea shortage of data in the context of estimation of single trial
ERPs because trial length is usually only about 1 second or
less. In fact, we tried to compare one more time-frequency
modeling method that does not use prior information withfrequency dependent spatial correlation matrices. However,the number of time frames of STFT was only eight, which
was too small to estimate 22 × 22 matrices. As a result, theEM algorithm never finished properly, forcing us to discard
this method. Utilizing prior information eﬀectively compensated for the shortage of data size and led to the better denoising performance of the proposed methods, as well as allows us to avoid the permutation ambiguity of event signals
because the target signal corresponds to the mixture compo-

Denoising performance of all compared methods on nontarget trials were mostly worse than target trials. As shownin Fig. 3, the amplitude of P300 that was embedded into resting state EEGs was smaller in non-target trials than target
trials as is often the case with an oddball paradigm experiment [18], which made SNR worse in non-target trials andcould make this noise removal problem harder than in targettrials.
7.

Conclusion

We developed a new noise removal method from ERP data
extending a conventional method that models generation ofeach event signal in the time frequency domain using priorinformation of ERP. Both of the proposed methods showed
eﬀectiveness in removing background EEG signals fromsingle-trial EEG signals to quantify ERP data more precisely
and lead us to better understanding of event-related brain dynamics.
There are a number of avenues for future work. First,
in this study we used EEG data to obtain prior informationof each event signal recorded in the same session and thesame subject. Ideally, we would like to use prior information collected from other days or subjects, which is a futurechallenge. Second, in this paper we used a carefully controlled experimental paradigm to avoid severe noise such as
body movements. Hence, we hope to further validate the effectiveness of our proposed methods to remove such severe
noise. Moreover, we used prior information about spatial
spread over the scalp by setting a prior distribution of spatial
correlation matrix Rk . However, we did not use any prior information about a timing of event signals. We can set a prior

MAKI et al.: ENHANCING ERPS BASED ON MAP ESTIMATION WITH A SPATIAL CORRELATION PRIOR

1445

distribution of the scaling parameter vk (n, f ) to exploit a timing prior, which may be eﬀective to do better estimation of
single-trial ERP because ERP components have characteristic timing of their amplitude shifts.
Acknowledgements
Part of this work was supported by JSPS KAKENHI Grant
Number 26540117.

[16] C.M. Bishop, Pattern recognition and machine learning, Springer
New York, 2006.
[17] American Clinical Neurophysiology Society, “Guideline Thirteen:
Guidelines for standard electrode position nomenclature,” J. Clinical
Neurophysiology, vol.11, no.1, pp.111–113, 1994.
[18] S.J. Luck, An introduction to the event-related potential technique,
MIT Press, 2014.
[19] J.W. Tukey, “Comparing individual means in the analysis of variance,” Biometrics, vol.5, no.2, pp.99–114, 1949.

References[1] M. Kutas, G. McCarthy, and E. Donchin, “Augmenting mental
chronometry: The P300 as a measure of stimulus evaluation time,”
Science, vol.197, no.4305, pp.792–795, 1977.
[2] B. Blankertz, R. Tomioka, S. Lemm, M. Kawanabe, and K.R.
Müller, “Optimizing spatial filters for robust EEG single-trial analysis,” IEEE Signal Process. Mag., vol.25, no.1, pp.41–56, 2008.
[3] T.-P. Jung, S. Makeig, M. Westerfield, J. Townsend, E. Courchesne,
and T.J. Sejnowski, “Analysis and visualization of single-trial
event-related potentials,” Human Brain Mapping, vol.14, no.3,
pp.166–185, 2001.
[4] B. Blankertz, S. Lemm, M. Treder, S. Haufe, and K.R. Müller,“Single-trial analysis and classification of ERP components — Atutorial,” NeuroImage, vol.56, no.2, pp.814–825, 2011.
[5] J.A. Urigüen and B. Garcia-Zapirain, “EEG artifact removal—
State-of-the-art and guidelines,” Journal of Neural Engineering,vol.12, no.3, 031001, 2015[6] D.M. Groppe, S. Makeig, and M. Kutas, “Independent component analysis of event-related potentials,” Cognitive Science Online,
2008.
[7] S. Lemm, G. Curio, Y. Hlushchuk, and K.R. Müller, “Enhancing
the signal-to-noise ratio of ICA-based extracted ERPs,” IEEE Trans.
Biomed. Eng., vol.53, no.4, pp.601–607, 2006.
[8] C. Févotte and J.F. Cardoso, “Maximum likelihood approach forblind audio source separation using time-frequency Gaussian source
models,” IEEE Workshop on Applications of Signal Processing to
Audio and Acoustics, pp.78–81, 2005.
[9] N.Q. Duong, E. Vincent, and R. Gribonval, “Under-determined reverberant audio source separation using a full-rank spatial covariance model,” IEEE Trans. Audio Speech Language Process., vol.18,
no.7, pp.1830–1840, 2010.
[10] Y. Kurihana, S. Miyabe, T.M. Rutkowski, Y. Matsumoto, T. Yamada,and S. Makino, “ Signal separation of EEG using multivariate probabilistic model,” IEICE Technical Report, MBE2012-119, 2013 (inJapanese).
[11] H. Maki, T. Toda, S. Sakti, G. Neubig, and S. Nakamura, “EEG
signal enhancement using multi-channel Wiener filter with a spatial
correlation prior,” Proc. IEEE International Conference on Acoustics, Speech and Signal Processing, pp.2639–2643, 2015.
[12] M. Zibulevsky and B.A. Pearlmutter, “Blind source separationby sparse decomposition in a signal dictionary,” Neural Comput.,
vol.13, no.4, pp.863–882, 2001.
[13] Y. Okada, J. Jung, and T. Kobayashi, “An automatic identification
and removal method for eye-blink artifacts in event-related magnetoencephalographic measurements,” Physiological Measurement,
vol.28, no.12, pp.1523–1532, 2007.
[14] F.C. Viola, J. Thorne, B. Edmonds, T. Schneider, T. Eichele, and
S. Debener, “Semi-automatic identification of independent components representing EEG artifact,” J. Clinical Neurophysiology,
vol.120, no.5, pp.868–877, 2009.
[15] A. Mognon, J. Jovicich, L. Bruzzone, and M. Buiatti, “ADJUST: An
automatic EEG artifact detector based on the joint use of spatial andtemporal features,” Psychophysiology, vol.48, no.2, pp.229–240,
2011.

Hayato Maki
received his B.S. from University of Tsukuba, Ibaraki, Japan in 2009, andrecieved his M.E. from Nara Institute of Scienceand Technology (NAIST), Nara, Japan in 2015.
He is currently a Ph.D. student at NAIST. He
is a graduate school member of IEEE SPS and
BME. His research interests include cognitive
science and signal processing with a focus on
biological signal processing.

Tomoki Toda
earned his B.E. degree from
Nagoya University in 1999 and his M.E. and
D.E. degrees from Nara Institute of Science andTechnology (NAIST) in 2001 and 2003, respectively. He was a Research Fellow of JSPS from2003 to 2005. He was an Assistant Professor
at NAIST from 2005 to 2011, and an Associate
Professor from 2011 to 2015. He is currently
a Professor at Nagoya University. His researchinterests include statistical approaches to speechprocessing. He received more than 10 paper
awards including the IEEE SPS 2009 Young Author Best Paper Awardand the 2013 EURASIP-ISCA Best Paper Award (Speech Communication
Journal).

Sakriani Sakti
received her B.E. degree
in Informatics (cum laude) from Bandung Institute of Technology, Indonesia, in 1999, andher M.Sc. degree from University of Ulm, Germany in 2002. While working in Japan at ATR(2003–2009) and NICT (2006–2011), she continued her study (2005–2008) with Dialog Systems Group University of Ulm, Germany, andreceived her Ph.D. degree in 2008. Currently,
she is an assistant professor of Nara Institute of
Science and Technology, Japan. Her researchinterests include statistical pattern recognition, speech recognition, spoken language translation, cognitive communication, and graphical modeling framework.

IEICE TRANS. INF. & SYST., VOL.E99–D, NO.6 JUNE 2016

1446

Graham Neubig
received his B.E.
from University of Illinois, Urbana-Champaign,
U.S.A, in 2005, and his M.E. and Ph.D. in informatics from Kyoto University, Kyoto, Japan in2010 and 2012 respectively. He is currently an
assistant professor at the Nara Institute of Science an Technology, Nara, Japan. His researchinterests include speech and natural language
processing, with a focus on machine learning
approaches for applications such as machinetranslation, speech recognition, and spoken dialog.

Satoshi Nakamura
is Professor of Graduate School of Information Science, Nara Institute of Science and Technology, Japan, Honorarprofessor of Karlsruhe Institute of Technology, Germany, and ATR Fellow. He received his
B.S. from Kyoto Institute of Technology in 1981
and Ph.D. from Kyoto University in 1992. He
was Associate Professor of Graduate School ofInformation Science at Nara Institute of Scienceand Technology in 1994-2000. He was Directorof ATR Spoken Language Communication Research Laboratories in 2000-2008 and Vice president of ATR in 2007-2008.
He was Director General of Keihanna Research Laboratories and the Executive Director of Knowledge Creating Communication Research Center,
National Institute of Information and Communications Technology, Japan
in 2009-2010. He is currently Director of Augmented Human Communication laboratory and a full professor of Graduate School of Information Science at Nara Institute of Science and Technology. He is interested in modeling and systems of speech-to-speech translation and speech recognition.
He is one of the leaders of speech-to-speech translation research and has
been serving for various speech-to-speech translation research projects inthe world including C-STAR, IWSLT and A-STAR. He received Yamashita
Research Award, Kiyasu Award from the Information Processing Societyof Japan, Telecom System Award, AAMT Nagao Award, Docomo Mobile
Science Award in 2007, ASJ Award for Distinguished Achievements in
Acoustics. He received the Commendation for Science and Technology by
the Minister of Education, Science and Technology, and the Commendationfor Science and Technology by the Minister of Internal Aﬀair and Communications. He also received LREC Antonio Zampoli Award 2012. He has
been Elected Board Member of International Speech Communication Association, ISCA, since June 2011, IEEE Signal Processing Magazine Editorial Board Member since April 2012, IEEE SPS Speech and Language
Technical Committee Member since 2013, and IEEE Fellow since 2016.