Removing Noise from Event-Related Potentials using a Probabilistic Generative
Model with Grouped Covariance Matrices
Hayato Maki1 , Tomoki Toda2,3 , Sakriani Sakti2 , Graham Neubig2 , and Satoshi Nakamura2

Abstract— Analysis of electroencephalograms (EEG) usually suffersfrom a variety of noises. In this paper, we propose a new method
for background noise removal from single-trial event-related potentials(ERPs) recorded with a multi-channel EEG. An observed signal isseparated into multiple signals with a multi-channel Wiener filter, whose
coefficients are estimated based on a probabilistic generative modelin the time-frequency domain. The main contribution is a method toestimate covariance matrices for each frequency bins of short-time
Fourier transform (STFT) representing different spatial spread of amulti-channel EEG signal according to frequencies. An experimentusing a pseudo-ERP data set demonstrates the effectiveness of our
proposed method.

I. INTRODUCTION
Event-related potentials (ERPs) are a brain response to a specific stimulus recorded by electroencephalograms (EEGs). ERPs
are usually averaged in time domain to weaken task-unrelated
background activities. However, a number of variabilities exist
across trials exist in latency, amplitude, and scalp distribution [1].
This across-trial averaging procedure conceals all the information
concerning across-trial variability of brain response. Moreover, toperform reliable analysis, at least 20 trials for each experimental
condition are generally averaged, resulting in long experimental
time and subject fatigue. Therefore, methods that allow for singletrial analysis without averaging are an important research topic
with the potential of making it possible to study event-related brain
dynamics more precisely [2]. Throughout this paper, we refer to the
ERPs as the signal to preserve and task-unrelated neural activities
as well as non-neural artifacts as noise to be removed.
Within the paradigm of EEG signal processing, there is a large
amount of research regarding external artifact removal, especially
eye blink removal [3]. However, removal of background EEGs isless studied, although there is some research proposing methods
based on independent component analysis (ICA) [4], [5].
In this paper, we propose a new approach to remove background
activities from single-trial ERPs based on an approach that has
been proposed originally in the field of sound source separation[6] and further applied to EEG signals [7]. In this method, model
parameters are estimated by the maximum likelihood criterion
without prior knowledge of target signal. We have extended thisapproach, exploiting prior knowledge of the target signal to further
improve the effectiveness of model parameter estimation[8]. Thisapproach has a virtue that, in contrast to ICA, it doesn’t assume
the number of sources because it does not estimate sources, but
contributions, of each event signal to each EEG electrode. While*Part of this work was supported by JSPS KAKENHI Grant Number
26540117.
*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
3 Tomoki Toda is also Faculty of the Information Technology Center,
Nagoya University, Japan

ICA aims to separate an observed signal into its sources, our method
separates it into signals elicited by each event such as ERP or aneye blinks.
We denote the group of source indexes that contributes to thek-th event by Ek and denote the EEG signal that is elicited by Ek
and captured at the scalp by ck (t) and call it k-th event signal. In
the time-frequency domain, short-time Fourier transform (STFT)
coefficients of the event signal ck (n, f ) are locally modeled by amultivariate complex Gaussian distribution whose parameters are afunction of (n, f ) where n is the index of the time frame and fis the index of the frequency bin. A multi-channel Wiener filter
is constructed using estimated model parameters to separate anobserved EEG signal into some event signal.
Model parameters includes full rank time-invariant covariancematrices called spatial correlation matrices encoding the spatial
spread of the event signals. In the previous works [7], [8], spatial
correlation matrices were assumed independent from frequency.
In this paper, we propose a generative model with frequency
dependent spatial correlation matrices because EEG signals usually
have different spatial spread according to frequencies.
II. E XISTING T ECHNIQUES FOR S IGNAL E NHANCEMENT
1) Observation Model: The STFT coefficients of the k-th event
signal in time-frequency slot (n, f ) are expressed as ck (n, f ):ck (n, f ) = ∑ hl sl (n, f ) = [ck,1 (n, f ), · · · , ck,L (n, f )]⊤

(1)

l∈Ek

where hl is the transfer function vector and sl (n, f ) is the l-th source
activity in the slot (n, f ) . The observed multi-channel EEG signalx(n, f ) is expressed asK

x(n, f ) = [x1 (n, f ), · · · , xL (n, f )]⊤ = ∑ ck (n, f ),

(2)

k=1

where K is a hyper parameter for the total number of events.
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 :p (sl (n, f )) = Nc (sl (n, f ); 0, vk (n, f )λl )()|sl (n, f )|2
1
exp
,
=
π vk (n, f )λlvk (n, f )λl(3)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 eachIn the M-step, α̂k and v̂k (n, f ) are updated as follows:

time-frequency slot (n, f ):E[sl1 (n, f )s∗l2 (n, f )] = 0 for l1 ̸= l2 ,(4)where l1 , l2 ∈ Ek and {·}∗ indicates the complex conjugate ofsl2 (n, f ).
From the assumptions above, the probability density functionthat each event signal follows is the multivariate complex normal
distribution with zero mean,()
p (ck (n, f )) = Nc ck (n, f ); 0, Rck (n, f ))(exp −ck (n, f )H R−1
ck (n, f )ck (n, f )=.
π I det(Rck (n, f ))(5)

For simplicity, we restrict the covariance matrix Rck (n, f ) to the
product of the time-invariant covariance matrix Rk that encodes
spatial spread of the k-th event signal, and the time-frequency
variant scalar vk (n, f ) that encodes time-frequency power of thesignal:[]Rck (n, f ) = E ck (n, f )ck (n, f )H()H 
= E  ∑ hl1 sl1 (n, f )l1 ∈Ek∑ hl sl (n, f )l2 ∈Ek2

2= ∑ E[|sl (n, f )|2 ]hl h⊤ll∈Ek= ∑ vk (n, f )λl hl h⊤l = vk (n, f )Rk ,

(6)

l∈Ek

where {·}H is the complex conjugate transpose. We call Rk the
spatial correlation matrix of the k-th event signal.
To simplify the generative model, we assume that the events ineach 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 particularelement zk (n, f ) is equal to 1 and all other elements are equal to 0.
If we observe the l-th event signal in the slot (n, f ), we can express
the observed signal assuming zl (n, f ) = 1 as follows:Kx(n, f ) = ∑ zk (n, f )ck (n, f ) = cl (n, f ),

(7)

k=1

From all assumptions above, the likelihood of the observed
signals x in the time-frequency slot (n, f ) is expressed by a Gaussian
mixture model as follows:Kp(x(n, f )|θ ) = ∑ αk Nc (x(n, f ); 0, vk (n, f )Rk )

(8)

k=1

where αk is the prior probability that the k-th event signal activates.
The model parameter set θ consists of αk , vk (n, f ), and Rk of each
mixture component.
The optimal model parameter set θ is estimated by maximizing
the observation likelihood described in eq. (8). This maximization
process can be effectively solved with the EM algorithm.
In the E-step, the posterior probability mk (n, f ) is calculated ateach time-frequency slot:
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′ ).

(9)∑n, f mk (n, f ),∑n, f ,k′ mk′ (n, f )1v̂k (n, f ) = x(n, f )H R−1
k x(n, f ),Lm (n, f )∑n, f v̂kk(n, f ) x(n, f )x(n, f )HR̂k =,∑n, f mk (n, f )

α̂k =

(10)
(11)
(12)

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 as follows:ĉk (n, f ) = mk (n, f )v̂k (n, f )R̂k R−1
x (n, f )x(n, f ).

(13)

Utilizing Spatial Correlation Prior and MAP Estimaition
We often know which event signal we would like to enhance,such as ERP, and we can record EEG signals related to thetarget event beforehand and use them as prior knowledge forenhancement. 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 described in the previous section.
Hence, we have proposed utilizing prior information of event
signals to estimate spatial correlation matrices [8] utilizing Wishart
distributions, whose probability density function is given by
)
(
q−p−1
Tr(W−1 R)
W p (R|W, q) = B(W, p, q)|R| 2 exp −(14)2
where W is a p × p positive definite matrix, q is called the numberof degrees of freedom and is restricted to q ≥ p, and B(W, p, q)
is a normalizing factor. The Wishart distribution is the conjugate
prior distribution of a p-dimensional multivariate Gaussian randomvariable’s precision matrix with a known mean vector.
From Bayes’ theorem, the posterior probability density functionof the time-frequency invariant spatial covariance matrices R ={R1 , · · · , RK } is proportional to the product of the likelihood of
the observed data and their prior as follows:()
K−1
p(R|x, θ\R , Ψk , qk ) ∝ ∏ p (x(n, f )|θ ) ∏ WL R−1
|Ψ
,
q
(15)
k
k
k
n, f

k=1

where θ\R is the model parameter set except for R.
The optimal model parameter set θ is estimated by maximizing
the posterior probability described in eq. (15) instead of observation
likelihood described in eq. (8). This maximization process can also
be effectively solved with the EM algorithm.
In the E-step, the posterior probability mk (n, f ) is calculated ateach time-frequency slot as shown in eq. (10). In the M-step, α̂kand v̂k (n, f ) are updated as shown in eq. (10) and eq. (11). The
spatial correlation matrices are updated as follows:)(m (n, f )
Ψk + 2 ∑n, f v̂ k(n, f ) x(n, f )x(n, f )H
k
R̂k =
.
(16)
(qk − L − 1) + 2 ∑n, f mk (n, f )III. P ROPOSED M ETHOD
A. Frequency Dependent Spatial Correlation Matrix
Spatial covariance matrices, which represent the spatial spread ofmulti-channel EEG signals in the generative model of the existing
work, are independent from frequency. However, EEG signals
usually have different spatial spread over the scalp according to thefrequency, as shown in Fig. 2, so a probabilistic model considering

6.1

10.0

22.0 Hz+

6.1 Hzthe frequency difference has the potential to be more effective. In
this paper, we propose a generative model with frequency dependent
spatial correlation matrices Rk ( f ) instead of frequency independent
matrix Rk of the previous works. In this manner, the estimation of
spatial covariance matrices in the EM algorithm is done as follows:)(m (n, f )Ψk ( f ) + 2 ∑n v̂ k(n, f ) x(n, f )x(n, f )HkR̂k ( f ) =.
(17)(qk ( f ) − L − 1) + 2 ∑n mk (n, f )22.0 Hz
22.0
22.0Hz
Hz

High
+

2020

--

-

Low

1010
20
20 20

00

Power 10*log 10 (µV22 /Hz)Power 10*log 10(µV /Hz)Power 10*log 10 (µV2 /Hz)Power 10*log 10 (µV2 /Hz)Power 10*log10(µV2/Hz)-

B. Grouped Spatial Correlation Matrix

10
10 10

-10
-10
00

-20
-20

0

-10
-10 -10
-30

-30

-20
-20 -20

-40
-40

-30
-30 -30

-50
-50

5

10

5

-40
-40 -40

10

15

55

5

10
10 10

20 2
0
20

15

Frequency (Hz)Frequency (Hz)-50
-50 -50

15
15 15

20
20 20

Frequency
Frequency
Frequency(Hz)(Hz)(Hz)

Fig. 1. Scalp topographies within the different frequency bins. Each line
stands for an EEG channel.
Differences of normalized ℓ1 norm

In the frequency dependent covariance model described in theprevious section, the number of samples x(n, f ) to estimate a
covariance matrix is reduced to 1/F compared to the conventional
method described in Section II, as the samples that belong to only
one frequency bin are involved in estimating the corresponding
covariance matrix as shown in eq. (17). This has an adverse effecton robust covariance matrix estimation although it can effectively
represent a frequency dependent spatial spread of a multi-channel
EEG signal, while all samples x(n, f ) in the time-frequency gridare involved in the conventional method. Therefore, there is atradeoff between frequency specific covariance matrix modelingand robust covariance matrix estimation. To solve this problem,
we group covariance matrices from multiple frequency bins together to share their statistics. The difference of the normalizedℓ1 norm of covariance matrices for the adjacent frequencies ofa multi-channel EEG signal are shown in Fig. 2 to show how
much covariance matrices of frequencies differ. We can see thatcovariance matrices corresponding to higher frequencies, except for
the Nyquist frequency have a relatively small difference between
adjacent matrices. Based on this, we estimate one covariance matrix
corresponding to higher frequencies while individually estimating
covariance matrices corresponding to lower frequencies. In thismanner, the frequency bins are grouped as lower and higher groups
as:22.0 Hz

10.0
10.0 10.0++

Note that the range of summation is over only time frame n.

19
18
17
16
15
14
13
12
11
10
0

10

20

30

40

50

60

Indexes of Frequency

Fig. 2. The differences of normalized ℓ 1 norm of adjacent covariancematrices. The last point stands for the Nyquist frequency.

third session, he was told to relax without any body movement torecord resting state EEG for two minutes.
B. Creation of Validation Dataset{ f1 , · · · , fF } = {{ f1 , · · · , fg }, { fg+1 , · · · , fF }}= {glow , ghigh }10.0 Hz

6.1
6.1 6.1(18)

where fg is an index of a frequency bin. The estimation of spatial
covariance matrices in the EM algorithm is done as follows:()
m (n, f )Ψk ( f ) + 2 ∑n ∑ f ∈ghigh v̂ k(n, f ) x(n, f )x(n, f )HkR̂k ( f ) =(19)(qk ( f ) − L − 1) + 2 ∑n ∑ f ∈ghigh mk (n, f )if f ∈ ghigh , otherwise Rk ( f ) is updated by eq. (17).
IV. E XPERIMENTAL E VALUATION
A. Experimental Paradigm
EEG data were collected from a healthy male volunteer aged
24 years without any neurological disorders. All EEG signals wererecorded from 20 scalp electrodes at locations based on a modified
International 10-20 system [9], digitally sampled at 250 Hz.
We conducted three sessions. In the first and second sessions, we
conducted auditory oddball paradigm experiments [10] using 1kHz
and 2kHz sound. A random sequence of auditory stimuli including
2kHz and 1kHz sine waves was presented to the subject. 1kHz
sounds were presented 200 times as non-target stimuli and 2kHz
sine wave were presented 50 times as target stimuli for each session.
The subjects were told to count the number of target stimuli without
any body movement. It is widely known that the ERP componentof P300 will be elicited when a subject hears target stimuli. In the

We created a pseudo-ERP dataset in order to evaluate noise
removal methods objectively. First, we cut the multi-channel EEG
data obtained from session 2 for each type of stimulus into trialsof 1.2 second length, then averaged them across the trials for eachtype of stimulus. Next, we cut the resting state EEG obtained fromsession 3 into 40 signals of 1.2 second length, and then we added
the averaged ERP of each corresponding channel that stands fortrue ERP to the 20 resting state trials for each type of stimuli thatstand for noise with randomly shifted phase from -60ms to 60ms.
C. Evaluation
We compared three denoising methods, ICA [4] and two timefrequency modeling methods that have spatial correlation matrices
dependent on individual frequencies and on grouped frequencies
where a covariance matrix corresponding to frequencies higher than
25th frequency bins (about 48 Hz) is shared. All twenty trials
for each type of stimuli are concatenated as one signal that weredenoised by each compared method, and then re-epoched. An ERP
signal in time domain x(t) is often quantified by its peak amplitudex(t peak ) and latency t peak , the former is its maximum value of
amplitude within the given time window and the latter is the elapsed
time from the onset of stimuli to reach the peak amplitude. We
quantified the P300 components with the time window from 200ms
to 500ms from the onset of stimuli. After quantification of the eachsingle-trial ERP (true ERP + noise) and true ERPs or averaged

The Ratio of Samples Correctly Classified

30 30

ICA

30 30

Peak Amplitude (µV)

Peak Amplitude (µV)

Without Process

20 20
20 20
30 30
30 30
10 10
10 10
20 20
20 20
0 0
0 0
10 10
10 10
-10 -10-10 -100 300300 400400 500500 200
0 300300 400400 500500
0 200
0 200
200
Peak Latency (ms)Peak Latency (ms)

-10 -10-10 -10200Filter
200 300
300 400400 500500
30
30
30
200
300300 400400 50050030
200
Wiener
WienerFilter+ Grouped Covariance

Peak Amplitude (μV)

Peak Amplitude (µV)

+ Individual Covariance

20 20
20 20
30 30
30 30
10 10
10 10
20 20
20 20
0 0
0 0
10 10
10 10
-10 -10-10 -100 300300 400400 500500 200
0 300300 400400 500500
0 200
0 200
200

-10 -10-10 -10
200200 300300 400400 500500 200200 300300 400400 500500
Peak Latency (ms)Peak Latency (ms): True target ERP: True non-target ERP: Single-trial target ERPs

: Single-trial non-target ERPs

Fig. 3. Plots of peak amplitudes and latencies. Red circles stands for singletrial ERPs of target trials and blue triangles for non-target. Two filled points
stands for true ERPs for each type of stimuli.

ERP signals of each type of stimuli obtained from session 2 andadded to the resting state EEGs, we calculated amplitude (AD) andlatency deviation (LD) for each single-trial from the true ERPs as
follows:
AD( j) = x̄(t¯peak ) − x̂( j) (t peak )

(20)

( j)
LD( j) = t¯peak − tˆpeak(21)

where x̄ is the ’true’ ERP and x̂ is the j-th single-trial ERP. Using
these deviation measures, we calculated Mahalanobis distance fromeach single-trial point to two true points to classify each singletrial point to the nearer true point, and then counted the number ofsingle-trials that were correctly classified.
D. Parameter Settings
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 first event signal
was calculated from the signal that were given by the session 1.
Hyper parameters for the second event signal were calculated fromthe resting state EEG signal recorded in the session 3.
V. R ESULTS
The peak amplitudes and latencies of the true ERP and of theeach single-trial are plotted in Fig. 3. When single-trial ERPs arewithout processing, the points in the figure are interspersed (upper
left). By the denoising effect of proposed methods (lower), thescattered points gathered around their true point for each type ofstimuli.
The ratio of single-trials correctly classified is shown in Fig. 4.
When single-trial ERPs are without processing, only 60% of them

1
0.9
0.8
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

1

2

3

4

1: Without Process
2: Independent Component Analysis
3: Wiener Filter + Individual Covariance
4: Wiener Filter + Grouped Covariance

Fig. 4.

The ratio of forty single-trials that are correctly classified.

are correctly classified. After denoising by the proposed method(Wiener filter + Grouped Covariance) 82.5% are correctly classified.
VI. C ONCLUSION
We developed a new noise removal method from ERP data
extending a conventional method that models generation of each
event signal in the time frequency domain with frequency dependent
covariance matrices with prior information of ERP. Both of the
proposed methods, especially with the grouped covariance model
showed effectiveness in removing background EEG signals fromsingle-trial ERPs.
We used EEG data to obtain prior information of each event
signal recorded in the same day and the same subject. Ideally, priorinformation collected from other days or subjects should be used
for a ready-to-use system in the future work.
R EFERENCES[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] T-P. Jung, S. Makeig, M. Westerfield, J. Townsend, E. Courchesne, andT. J Sejnowski. Analysis and visualization of single-trial event-related
potentials. Human brain mapping, Vol. 14, No. 3, pp. 166–185, 2001.
[3] J. A. Urigüen and B. Garcia-Zapirain. EEG artifact removal - Stateof-the-art and guidelines. Journal of Neural Engineering, 2015.
[4] D. M Groppe, S. Makeig, and M. Kutas. Independent component
analysis of event-related potentials. Cognitive science online, Vol. 53,
No. 4, pp. 601–607, 2008.
[5] S. Lemm, G. Curio, Y. Hlushchuk, and K. R. Müller. Enhancing thesignal-to-noise ratio of ICA-based extracted ERPs. IEEE Transactions
on Biomedical Engineering, Vol. 53, No. 4, pp. 601–607, 2006.
[6] N. QK Duong, E. Vincent, and R. Gribonval. Under-determined
reverberant audio source separation using a full-rank spatial covariance model. IEEE Transactions on Audio, Speech, and Language
Processing, Vol. 18, No. 7, pp. 1830–1840, 2010.
[7] Y. Kurihana, S. Miyabe, T. M. Rutkowski, Y. Matsumoto, T. Yamada,and S. Makino. (in Japanese) Signal Separation of EEG using
Multivariate Probabilistic Model. Technical Report of IEICE, MBE,
Vol. 112, No. 479, pp. 161–166, 2013.
[8] 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 ICASSP, 2015.
[9] American Clinical Neurophysiology Society. Guidelines for standard
electrode position nomenclature. Clinical Neurophysiology, 1994.
[10] S. J Luck. An introduction to the event-related potential technique.
MIT press, 2014.