NeuroImage 51 (2010) 629–641

Contents lists available at ScienceDirect

NeuroImagej o u r n a l h o m e p a g e : w w w. e l s e v i e r. c o m / l o c a t e / y n i m g

A generalized method to estimate waveforms common across trials from EEGs
Yusuke Takeda a,⁎, Masa-aki Sato a, Kentaro Yamanaka b, Daichi Nozaki c, Yoshiharu Yamamoto ca

ATR Computational Neuroscience Laboratories, 2-2-2 Hikaridai, Keihanna Science City, Kyoto 619-0288, Japan
Department of Health Design, Showa Women's University, 1-7 Taishido, Setagaya-ku, Tokyo 154-8533, Japanc
Educational Physiology Laboratory, Graduate School of Education, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-0033, Japanba r t i c l e

i n f oArticle history:
Received 29 April 2009
Revised 24 December 2009
Accepted 2 February 2010
Available online 10 February 2010
Keywords:
Electroencephalography (EEG)
EEG waveforms common across trials
Go/NoGo taska b s t r a c t
We propose a generalized method to estimate waveforms common across trials from electroencephalographic (EEG) data. From single/multi-channel EEGs, the proposed method estimates the number ofwaveforms common across trials, their delays in individual trials, and all of the waveforms. After verifying
the performance of this method by a number of simulation tests with artiﬁcial EEGs, we apply it to EEGs
during a Go/NoGo task. This method can be used in general situations where the number and the delays of
EEG waveforms common across trials are unknown.
© 2010 Elsevier Inc. All rights reserved.

Introduction
In many electroencephalographic (EEG) studies, EEG waveforms
common across trials have been estimated by averaging EEG epochs
across trials. For example, evoked potentials, such as visual evoked
potentials (VEPs), are estimated by averaging EEG epochs that are
triggered on visual stimulus onsets. When waveforms temporally
overlap, however, the averaging procedure cannot estimate exact
waveforms because they are mutually contaminated (Kok, 1988;
Takeda et al., 2008a; Verleger, 1988). Further, when the delays ofwaveforms are variable and unknown, the averaging procedure cannot
be used (Takeda et al., 2008b; Tallon-Baudry and Bertrand, 1999).
Recently, we proposed two methods that overcome thoselimitations (Takeda et al., 2008a; Takeda et al., 2008b). These two
methods assume that two waveforms common across trials exist in an
EEG epoch and estimate them from single-channel EEG epochs. One
method is used when the delays of two waveforms are given (Takeda
et al., 2008a). By this method, we obtain pure waveforms that are notcontaminated with each other from EEGs during stimulus–responsetasks, in which the delays of two waveforms are given from stimulus
and response onsets. The other method is used when the delays of a
waveform are not given (Takeda et al., 2008b). By this method, we canobtain the delays as well as the pure waveforms from EEGs during
covert response tasks, such as decision-making tasks.
However, these methods lack versatility. In particular, the validityof the assumption, the existence of two EEG waveforms common
across trials, is not always guaranteed. Generally, the total number ofwaveforms is unknown. It is possible that three or more waveforms

exist. In fact, the study of Verleger et al. (2005) is suggestive of three
waveforms in EEGs during Go trials of a Go/NoGo task: a stimuluslocked waveform, a response-locked waveform, and a waveform timelocked to neither stimulus nor motor response onsets. In such a case,applying the above methods is inappropriate. On the other hand, theprevious methods use only single-channel EEGs. While this property isan advantage when only single-channel EEGs are available, it
sometimes becomes a disadvantage when multi-channel EEGs areavailable because using all available EEGs may provide more
information than just using a small portion of it. To investigate various
types of EEGs in detail, we need a more general method that can deal
with an unknown number of waveforms and multi-channel EEGs.
In this paper, we propose a generalized method to estimate EEG
waveforms common across trials. From single/multi-channel EEGs,
the proposed method estimates the number of waveforms common
across trials, their delays in individual trials, and all of the waveforms.
We examine the performance of this method by a number ofsimulation tests with artiﬁcial EEGs. Then as an example, we apply
this method to EEGs during a Go/NoGo task.
Methods
Proposed method to estimate waveforms common across trials
Assumption and purpose
An EEG epoch of a channel, which is assumed to consist of
waveforms common across trials and noise, is expressed byðchÞyn ðt Þ =⁎ Corresponding author. Fax: +81 774 95 1259.
E-mail address: takeda@atr.jp (Y. Takeda).
1053-8119/$ – see front matter © 2010 Elsevier Inc. All rights reserved.
doi:10.1016/j.neuroimage.2010.02.002

K
X

k=1

ðchÞ

sk

!
"
ðchÞ
t − τn;k + vn ðt Þ

ðt = 0; : : : ; T − 1; n = 1;: : :; N; ch = 1; : : : ; CHÞ;

ð1Þ

630

Y. Takeda et al. / NeuroImage 51 (2010) 629–641

The Evaluation procedure is based on the results of Delay estimationand Waveform estimation. The procedure of Delay estimation is based onthat of Waveform estimation. Therefore, we describe these steps in thefollowing order: Waveform estimation, Delay estimation, and Evaluation.
Waveform estimation
In this step, we estimate S(ch)(ω) from Y(ch)(ω) when τ and K are
given.
In the least squares method, estimated S(ch)(ω) is expressed byŜðchÞ ðωÞ = argmin j jY
ðchÞ

S

ðchÞ

ðωÞ

ðchÞ

ðωÞ − Eðω; τÞS2ðωÞj j ;where argmin : : : represents x that minimizes ⋯ and || ⋯ || representsχ

the norm of ⋯. By solving
Fig. 1. Flowchart to estimate number of waveforms common across trials, their delays
in individual trials, and all of the waveforms.

Aj jYðchÞ ðωÞ − Eðω; τÞSðchÞ ðωÞj j 2= 0;ASðchÞ ðωÞ

we obtain(ch)where y(ch)(t): kn (t): observed EEG epoch of channel ch in trial n, sk(ch)th waveform of channel ch, τn,k: delay of sk (t) in trial n, v(ch)n (t):
noise of channel ch in trial n, and K: number of waveforms. Noise v(ch)n(t) is assumed to be a stationary process.
For simplicity, we rewrite Eq. (1) as a matrix representation in the
Fourier domain as below. By taking the discrete Fourier transform ofEq. (1), we obtainðchÞ

Yn ðωÞ=

K
X

k=1

!
"
ðchÞ
ðchÞ
exp − i2πωτn;k = T Sk ðωÞ+ Vn ðωÞðω =0; : : : ; T − 1Þ;

ðchÞ

ðchÞ

ðωÞ = Eðω; τÞSwhereðωÞ + V

ðchÞ

ðωÞ;

h
iT
ðchÞ
ðchÞ
YðchÞ ðωÞ = Y1 ðωÞ; : : : ; YN ðωÞ ;
2
!
"
:::
exp − i2πωτ1;1 = T
66Eðω; τ Þ = 6
v
O
6
!
"
4
:
::
exp − i2πωτ = T
N;1

ð3Þ

!
"3
exp − i2πωτ1;K = T
7
7
7;
v
7
!
"5
exp − i2πωτ N;K = T

τ = τ n;k ðn = 1; : : : ; N; k = 1; : : : ; K Þ;
h
iT
ðchÞ
ðchÞ
SðchÞ ðωÞ = S1 ðωÞ;: : : ; SK ðωÞ ;

and
V

ðchÞ

ðωÞ =

h

i

T
ðchÞ
ðchÞ
V1 ðωÞ;: : :; VN ðωÞ :The purpose of the proposed method is to obtain Ŝ (ch)(ω)τ̂ and K̂ only
from Y (ch)(ω), where Ŝ (ch)(ω): the estimated S(ch)(ω),τ̂: the estimated τ, and K̂: the estimated K.
Overview of procedure for obtaining Ŝ (ch)(ω),τ̂ and K̂
The proposed method consists of three steps: Delay estimation,
Waveform estimation, and Evaluation, which are consecutively repeatedas shown in Fig. 1. At the beginning of the procedure, we set K to 1. In
the Delay estimation step, we obtain τ̂ using preset K. In the Waveform
estimation step, we obtain Ŝ (ch)(ω) using τ̂ and preset K. In the Evaluation step, we evaluate whether preset K is true using Ŝ (ch)(ω) and τ̂.
If preset K is evaluated as wrong, we return to the Delay estimation stepby increasing preset K by 1. By repeating this procedure until preset K
is evaluated as true, we can simultaneously obtain Ŝ (ch)(ω),τ̂, and K̂.

ð4ÞWhen ω = 0, all of the values in E(ω,τ) become 1, and E(ω,τ)TE(ω,τ)
becomes singular. This corresponds to the fact that arbitrary constants
can be added to ŝ k(ch)(t). Therefore, we set the time average of ŝ (ch)(t)
k
to be 0 byŜðchÞ ðωÞ =

(

ð2Þ

(ch)
(ω), and Vn(ch)(ω): the discrete Fourier transforms
where Y(ch)n (ω), Sk(ch)of y(ch)
(t),
s
(t),
and
v(ch)
n
k
n (t), respectively. Eq. (2) is rewritten as

Y

h
i−1
T
T ðchÞ
ŜðchÞ ðωÞ = Eðω; τÞ Eðω; τÞEðω; τ Þ Y ðωÞ:

h0
i−1
T
T ðchÞ
Eðω; τ Þ Eðω; τ ÞEðω; τ Þ Y ðωÞω=0ω ≠ 0:

ð5Þ

We can obtain ŝ k(ch)(t) by taking the inverse discrete Fourier transformof Ŝ (ch)(ω).
Delay estimation
In this step, we estimate τ from Y(ch)(ω) when K is given.
In the least squares method, estimated τ is expressed byτ̂ = argmin
τ

T =2
CH
X
X

ch = 1 ω = 1

j jY

ðchÞ

ðchÞ

ðωÞ − Eðω; τ ÞS2ðωÞj j :ð6Þ

In this equation, since S(ch)(ω) is neither known nor determined by τ,
we replace it by Ŝ (ch)(ω) in Eq. (5), which is the least squares solutionof S(ch)(ω) determined by τ. Then Eq. (6) is rewritten asτ̂ = argmin
τ

T =2
CH
X
X

ch = 1 ω = 1
T ðchÞ

× Eðω; τ Þ Y

j jY2ðωÞ j j :

ðchÞ

h
i−1
T
ðωÞ − Eðω; τ Þ Eðω; τÞ Eðω; τÞð7Þ

By solving Eq. (7), we can estimate τ. However, because E(ω,τ)[E(ω,τ)TE(ω,τ)]-1E(ω,τ)T in Eq. (7) is nonlinear with respect to τ, wecannot solve Eq. (7) analytically. We solve it by a hybrid optimization
algorithm that consists of two consecutive stages: global search andlocal search. First, we obtain an approximate solution by global searchand then obtain the optimal solution by local search. In global search,
we conduct a random search (Zhigljavsky, 1991) (see Appendix A)
modiﬁed from that in our previous study (Takeda et al., 2008b) for M(=50) times with a different initial τ. Then we obtainM sets of τ iandh−1
PCH
PT = 2ðchÞoτ( =ðωÞ − Eðω; τÞ Eðω; τÞT Eðω; τÞch = 1
ω = 1 j jY
T ðchÞ
2
Eðω; τ Þ YðωÞj j ) and select the τ that minimizes oτ. In local
search, we conduct a grid search (see Appendix A) by setting the τ
selected in global search as initial τ.

631

Y. Takeda et al. / NeuroImage 51 (2010) 629–641

After optimization, we adjusted the averages across n of obtained τ.
This adjustment is required because the averages of τ vary dependingon the time points deﬁned as the onsets of the waveforms, and theonsets are arbitrarily determined in optimization. For example, we
adjusted the average of τn,2 obtained from the EEGs during a NoGo task
so that a peak in the estimated waveform-2 became its onsets(described in Data analysis section). Adjusted τ is referred to as τ̂.
Evaluation
In this step, we evaluate whether preset K is true.
As described above, once K is set, we can obtain τ̂ and Ŝ (ch)(ω) by
the steps of Delay estimation and Waveform estimation, respectively,
and we can obtain the residual error between observed and reconstructed EEGs by
"
#
K
X
ðchÞ
ðchÞ
ðchÞ
ren ðt Þ = IDFT Yn ðωÞ −expð − i2πω τ̂n;k = TÞ Ŝk ðωÞ ; ð8Þ
k=1

where IDFT[⋯]: inverse discrete Fourier transform of ⋯. We evaluatepreset K by examining re(ch)n (t). When preset K is smaller than true,
(ch)
s(ch)
(t) should remain in re(ch)kn (t). As a result, ren (t) should be anonstationary process. In contrast, when preset K is true, s(ch)
(t)
k
(ch)
should disappear from re(ch)n (t). As a result, ren (t) should be astationary process. Therefore, we evaluate preset K by examining
whether re(ch)n (t) is a stationary process. We examine whether the
distribution of re(ch)n (t) differs before and after stimulus onsets. We
divide all re(ch)(t) (n = 1,⋯, N;ch = 1,⋯,CH) into two samples
ncorresponding to before and after stimulus onsets, and conduct a
two-tailed Two-Sample Kolmogorov–Smirnov test to test the null
hypothesis that the two samples are drawn from the same
distribution. The probability of p b 0.05 is accepted as signiﬁcant. If
the null hypothesis is rejected, we regard preset K as wrong. If the null
hypothesis comes to be accepted, we regard preset K as true.
Simulation tests
To examine the performance of the proposed method, we
conducted simulation tests for Waveform estimation, Delay estimation,
and Evaluation.
In these simulation tests, we generated simulated data y(ch)n (t) asfollows. We generated ﬁve waveforms, s1(t),⋯s5(t), by an exponentialfunction, a cosine function, a rectangular function, a sawtooth wave,
and a triangular pulse, respectively. All of the waveforms wereidentical regardless of channels ch. We set the variance across t of
the waveforms to 1. We set the delays of s1(t) to 1 and generated thedelays of the other waveforms by Gaussian random numbers[mean = 18, standard deviation (SD) = 5]. We used white noise as
(ch)
(ch)
noise v(ch)n (t). Then we generated yn (t) from sk(t), τ and vn (t),based on Eq. (1).
Simulation tests for Waveform estimation
To evaluate the quality of Waveform estimation, we scrutinized theresidual errors between estimated waveforms ŝ k(t) and original
waveforms sk(t).
First, we examined whether the residual errors ﬂuctuated randomly
around 0. We generated simulated data yn(ch)(t) with the following
parameters: number of waveforms K = 2, number of channels CH= 1,
number of trials N = 100, and signal-to-noise ratio (SNR)= 6. SNR was
deﬁned as 10log10(Var[sk(t)]/Var[vn(t)]), where Var[sk(t)]= 1. We
repeated the estimation of sk(t) 500 times from different sets of yn(ch)(t) and obtained 500 sets of residual errors. Because the averages
across time of ŝk(t) are inherently indeﬁnite, we adjusted the averages
across time of the residual errors to 0. We then plotted the time coursesof the means and SDs of the 500 residual errors.
Second, we examined whether the magnitudes of the residual
errors became smaller as the number of trials N increased. We

generated y(ch)n (t) with the following parameters: K = 2, CH = 1,
N = 100,⋯,1000, and SNR = 6. For each N, we repeated the extractionof sk(t) 500 times from different sets of y(ch)n (t) and obtained 500 setsof residual errors. The magnitudes of the residual errors were
quantiﬁed by the variance across time of the residual errors:
500 X
K
h
i
X
1
ðp;N Þðt Þ − sk ðt Þ ;Var ŝk
500 × K p = 1 k = 1

where ŝ k (p,N)(t) represents the p-th estimated waveform-k obtained
from the simulated data consisting of N trials. The variance was ﬁtted
by function y = a/N by the least squares method.
Finally, we compared error coefﬁcients a of ﬁtting function y = a/N obtained in different situations. a represents the magnitude of theresidual errors for all N, and thus a smaller value of a is better. We
generated y(ch)n (t) with the following parameters: K = 1,⋯,5, CH = 1,
N = 100,⋯,1000, and SNR = 6. For each K, we obtained a. Because theaveraging procedure can be used when K = 1, we also obtained a for
K = 1 using the averaging procedure. Because our previous method(Takeda et al., 2008a) can be used when K = 2, we also obtained a for
K = 2 using that method.
Simulation tests for Delay estimation
To evaluate the quality of Delay estimation, we examined thenormalized root mean squared errors (RMSEs) between sets of delays τ
obtained in the Delay estimation step and true τ. RMSE was calculated by
K
1X
K k=1
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi!
ffiffiffiffi!
ffiffiffiffiffi
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi!
u N
2
N
N
u1 X
1X
1X
t
true
true
×
τ n;k −
τ n;k − τn;k −
τn;k
;
N n=1
N n=1
N n=1

RMSE =

true
represent delay and true delay, respectively.
where τn,k and τn,kFirst, we examined the validity of Eq. (7) and whether τ becamecloser to the true one as optimization proceeded. We generated
simulated data y(ch)n (t) with the following parameters: number ofwaveforms K = 3, number of channels CH = 1, number of trials
N = 100, and SNR = 6. From y(ch)n (t) by global search, we obtained M(=50) sets of oτ and RMSEs, and by local search, we obtained one setof those. Then we made a scatter plot of oτ and the RMSEs.
Second, we examined the relation between estimation accuracy
and SNR. We generated y(ch)n (t) with the following parameters:
K = 2,3, CH = 1, N = 100, and SNR = -20,-15,-10,-5,0. For each K and
SNR, we estimated τ and calculated RMSE. The estimation was
repeated 10 times from different sets of y(ch)n (t).
Finally, we examined the relation between estimation accuracy
and CH. We generated y(ch)n (t) with the following parameters: K = 3,
CH = 1,2,3, N = 100, and SNR = -5. For each CH, we estimated τ andcalculated RMSE. The estimation was repeated 10 times from different
sets of y(ch)n (t). To examine the effect of adding meaningless channels,
we also generated y(ch)n (t) byðchÞyn ðt Þ =

8
K
!
"
X
>
ðchÞ
>
<
sk t − τ n;k + vn ðt Þ>>:

ch = 1k=1vðnchÞ ðt Þ;

ch N 1and calculated RMSE in the same way. In this simulation, we refer tochannels including and not including sk(t) as meaningful andmeaningless channels, respectively.
Simulation tests for Evaluation
In the simulation tests for Evaluation, we ﬁrst checked the rationaleof the procedure for evaluating the number of waveforms K. We
examined the time courses of the means and SDs across trials of the632

Y. Takeda et al. / NeuroImage 51 (2010) 629–641

Experimental procedure
The experimental population was comprised of nine healthy adults(age 28.4 ± 3.7 years), all of whom gave informed consents. The local
ethics committee approved the experimental procedure.
The subjects were comfortably seated on a chair in a dimly lit,
electrically shielded room. About 50 cm in front of their eyes, red andgreen light-emitting diodes (LEDs) for imperative signals werevertically arrayed 1.5 cm apart on a black panel. The subjects were
instructed to perform two tasks in the following order: a Go/NoGo task
and a passive viewing task. In the Go/NoGo task, they participated infour experimental blocks, each consisting of 50 trials. They were
instructed to push a button immediately after a “Go” signal (green LED)and to not push it after a “NoGo” signal (red LED). The green and red
LEDs were illuminated in random order at almost equal probability. In
two blocks, the subjects had to respond with their right index ﬁnger andin the other two blocks with their left index ﬁnger. In off-line analysis,
the data from the blocks of the left and right ﬁngers were mixed. In thepassive viewing task, the subjects participated in two experimental
blocks, each consisting of 50 trials. They were instructed to passively
view the same stimulus as in the Go/NoGo task. In both tasks, each trial
began with a warning signal (a beep), followed, after a variable delay of
1.8–2.2 s, by the imperative signals (duration: 500 ms). Inter-trial
intervals were randomized from 3.5 to 7.5 s.
During the tasks, surface EEGs were recorded from 19-ch tin
electrodes, mounted in a cap (Electro-Cap International, Inc., Eaton,
Ohio, USA) based on the international 10–20 system. A tin electrodeplaced on AFz was used as a ground. The EEGs were ampliﬁed on a
Nihon Kohden EEG-1100 with a time constant of 0.3 s. Because we
expected that large EEG activity related to the task execution would
not appear around the earlobes, we placed Ag/AgCl electrodes onboth earlobes and recorded their potentials separately. Their averaged
potentials were subtracted from the EEGs off-line. For monitoring eye
movements, an electrooculogram (EOG) was recorded with a pair of
Ag/AgCl electrodes placed above and below the left eye. The sampling
rate of the EEGs and EOGs was 1000 Hz.

at 50 Hz). Then the ﬁltered EEGs were segmented into 2 s epochs from-500 to 1500 ms after stimulus onsets.
If waveforms common across trials exist in EEGs, their mean and
SD across trials should change with time. The time courses of themeans and SDs of a subject's EEGs during the Go/NoGo task at Fz, C3,
C4, Cz, and Pz clearly showed transient changes, indicating theexistence of such waveforms in these EEGs. Therefore, we applied the
proposed method to them.
The reaction times (RTs) of the Go trials were deﬁned as the intervals
between the stimulus and the button-push signal (response) onsets. We
excluded the Go trials with an RT either shorter than 100 ms or longer
than 400 ms, and the NoGo trials with any response. An artifact criterionof ±50 μV was used for the EEGs and EOGs to reject trials with excess
ocular artifacts or measurement noise. We obtained 96 Go trials and 100
NoGo trials. The RTs were 275.63±46.24 (mean± SD) ms.
In the Delay estimation step, we estimated the delays of the EEG
waveforms common across trials using priori knowledge about thedelays. In the Go/NoGo task, the approximate delays of two waveforms
can be given from the stimulus and the response onsets; using thedelays simpliﬁes the optimization problem [Eq. (7)]. We assumed thatthe delays of waveform-1s slightly ﬂuctuate after the stimulus onsets
and the delays of waveform-2s in the Go trials slightly ﬂuctuate aroundthe response onsets. Therefore, based on reports that examined thevariable latencies of visual evoked potentials (Mihaylova et al., 1999;
Vassilev et al., 2002; Vaughan et al., 1966), we searched for the delays ofwaveform-1s from 0 to 50 ms after the stimulus onsets and the delays ofwaveform-2s in the Go trials from -25 to 25 ms after the response
onsets. We searched for the delays of the other waveforms setting theinitial delays to Gaussian random numbers [mean = 180 (ms), SD= 50(ms)]. We set the delay ranges at a width of 500 ms because the EEGs'
SDs, which are indications of the variable delays (Takeda et al., 2008b),
were clearly greater than the pre-stimulus level for about 500 ms after
the stimulus onset. After the optimization of Eq. (7), we adjusted thedelays of the estimated waveforms. The delays of waveform-1s were
adjusted so that the minimum values of the estimated delays became
0 ms. The delays of waveform-2s in the Go trials were adjusted so thatthe average of the estimated delays was identical with that of the RTs.
The delays of the other waveforms were adjusted so that the estimateddelays represented the latencies of the maximum positive peak in eachof the estimated waveforms at Cz. In the Waveform estimation step, weestimated the waveforms using the estimated delays. For comparison
with the estimated waveform-1s, we obtained stimulus-triggered
average EEGs by averaging the EEGs triggered on the stimulus onsets.
For comparison with the estimated waveform-2s in the Go trials, we also
obtained response-triggered average EEGs by averaging the EEGs during
the Go trials triggered on the response onsets. In the Evaluation step, we
evaluated whether the preset number of waveforms was true. We
compared the distributions of the residual errors at Fz, C3, C4, Cz, and Pz
from -0.5 to 0 s with those from 0 to 0.7 s after the stimulus onsets.
After the estimation, using the estimated delays, we estimated the
waveforms from the 19-channel EEGs by the Waveform estimation
procedure and obtained the scalp distributions of the variance acrosstime of the estimated waveforms. Finally, because the estimatedwaveforms seemed to have large oscillatory components, we calculated
the amplitude spectra by taking the discrete Fourier transform of h(t)ŝk(ch)(t), where h(t) was the Hanning window. We calculated theamplitude spectra at Fz, Cz, and Pz and averaged them across thechannels. For comparison, we also calculated the amplitude spectra of
the stimulus- and response-triggered average EEGs in the same way.

Data analysis

Results

In the off-line analysis, we resampled the EEGs at a rate of 100 Hz.
They were ﬁltered with a bandpass of 2–40 Hz using two kinds ofﬁnite impulse response (FIR) ﬁlters: a high-pass ﬁlter of 2 Hz (300point, -26 dB at 1 Hz) and a low-pass ﬁlter of 40 Hz (15-point, -45 dB

Simulation tests

residual errors between the original and reconstructed simulated
data. We generated simulated data y(ch)(t) with the followingnparameters: K = 3, number of channels CH = 1, number of trials
N = 100, and SNR = 6. From y(ch)n (t), we obtained three sets of residual
errors using a preset K of 1,2,3.
Then we examined the relation between the reliability of Evaluation and SNR. In the Evaluation step, we calculate the p-value of thenull hypothesis that re(ch)n (t) before and after stimulus onsets are
drawn from the same distribution. Correct K is estimated if the p-value
for K smaller than true is larger than 0.05 and the p-value for true K issmaller than 0.05. We obtained the probability of satisfying thiscondition in the following way. We generated y(ch)n (t) with thefollowing parameters: K = 3, CH = 3, N = 100, and SNR = -20,-15,10,-5,0. For each SNR, we conducted Delay estimation and Waveform
estimation with preset K = 2 (smaller than true) and K = 3 (true),obtained two sets of re(ch)n (t), and calculated the two p-values. We
repeated this procedure 10 times and obtained the probability ofsatisfying the condition.
Applications to EEGs during Go/NoGo task
As an example, we applied the proposed method to the EEGs
during a Go/NoGo task.

Figs. 2A–D shows original waveforms sk(t) (k = 1,2,3) and noise
(ch)
v(ch)
n (t). Fig. 2E shows simulated data yn (t) generated with the

Y. Takeda et al. / NeuroImage 51 (2010) 629–641

633

(ch)
(ch)
Fig. 2. Simulation tests for overall procedure. (A–D) Original waveforms sk(t) and noise v(ch)n (t). (A) s1(t). (B) s2(t). (C) s3(t). (D) vn (t). (E) Simulated data yn (t) generated froms1(t), s2(t), s3(t), and v(ch)
(t).
(F–H)
Estimatedwaveforms
ŝ
(t).
(F)
ŝ
(t).
(G)
ŝ
(t).
(H)
ŝ
(t).
n
k
1
2
3

following parameters: number of waveforms K = 3, number ofchannels CH = 1, number of trials N = 100, and SNR = 6. Before
proceeding to the simulation tests for the individual steps, we
examined whether the overall procedure (Fig. 1) worked well. Fromthe simulated data (Fig. 2E), we estimated K, the delays of waveforms
τ, and the sk(t) (Figs. 2F–H). Estimated K is correctly 3. Estimated sk(t)are highly correlated with the original ones; all of the correlation
coefﬁcients between the estimated and original sk(t) are higher than
0.99 (Figs. 2F–H). The RSMEs between the estimated and original τ

are 0, indicating that the estimation is completely accurate. These
results indicate that the proposed method, as a whole, works well forthe simulated data.
Simulation tests for Waveform estimation
Fig. 3 shows the results of the simulation tests for Waveform
estimation. Figs. 3A and B shows the time courses of the means and SDsof the residual errors between the estimated and original waveforms.
These time courses are nearly constant (Figs. 3A and B), indicating that

Fig. 3. Simulation tests for Waveform estimation. (A, B) Means (solid lines) and means ± SDs (dotted lines) of residual errors between estimated and original waveforms across 500
repeated simulations. (A) Result for estimated waveform-1 ŝ1(t). (B) Result for estimated waveform-2 ŝ2(t). (C) Variance of residual errors as function of number of trials N(diamonds). Solid lines represent ﬁtting curves in form of y = a/N. (D) Error coefﬁcient a of ﬁtting function y = a/N for each number of waveforms K.

634

Y. Takeda et al. / NeuroImage 51 (2010) 629–641

Fig. 4. Simulation tests for Delay estimation. (A) Scatter plots of values of objective function oτ and RMSEs between estimated and original delays τ. Open circle indicates selected τ inglobal search. Filled dots indicate non-selected τ in global search. Diamond indicates converged τ in following local search. (B) RMSEs between estimated and original τ for each SNR.
Solid line represents results for number of waveforms K = 3. Dotted line represents results for K = 2. (C) RMSEs between estimated and original τ for each number of channels CH.
Solid line represents results for adding meaningful channels. Dotted line represents results for adding meaningless channels. In (B) and (C), diamonds and error bars respectively
represent means and SDs of RMSEs.

the residual errors have no temporal modulation patterns andﬂuctuate randomly. Fig. 3C shows how the variance of the residual
errors changes as the number of trials N increases. The variance isinversely proportional to N (Fig. 3C), indicating that the residual errors
become smaller as the number of trials increases. Fig. 3D shows how
error coefﬁcients a of ﬁtting function y = a/N changes as the numberof waveforms K increases. a becomes larger as K increases, indicatingthat the accuracy of the estimated waveforms becomes lower as thenumber of waveforms increases. When K = 1, the a of the averaging
procedure is 0.25, which is the same as that of the Waveform estimation
procedure. This indicates that the accuracy of the Waveform estimation
procedure is identical as that of the averaging procedure. When K = 2,
the a of our previous method (Takeda et al., 2008a) is 1.91, which islarger than that of the Waveform estimation procedure. This indicatesthat the accuracy of the Waveform estimation procedure is higher thanthat of our previous method (Takeda et al., 2008a).

non-selected τ (ﬁlled dots), indicating the validity of global search.
RMSE of τ̂ obtained by local search (diamond) is smaller than that of
the selected τ in global search (open circle), indicating the validity
and necessity of local search. Fig. 4B shows the relations between
RMSEs and SNR. RMSEs become smaller as SNR becomes higher,
indicating that the estimation accuracy becomes higher as SNR
becomes higher. For SNRs of -15,-10,-5,0, RMSEs for the number ofwaveforms K = 2 (dotted line) are signiﬁcantly smaller than those for
K = 3 (solid line) (p b 0.05, two-tailed Mann–Whitney test), indicatingthat the estimation accuracy for K = 2 is higher than that for K = 3.
Fig. 4C shows the relations between RMSEs and the number ofchannels CH. While RMSEs obtained by adding meaningful channels(solid line) become smaller as CH increases, those by adding
meaningless channels (dotted line) become larger. This indicatesthat adding meaningful channels increases estimation accuracy but
not adding meaningless channels.

Simulation tests for Delay estimation
Fig. 4 shows the results of the simulation tests for Delay estimation.
Fig. 4A is the scatter plot of RMSEs between the estimated and original
delays τ and the values of objective function oτ. RMSE becomes
smaller as oτ becomes smaller, suggesting the validity of Eq. (7). RMSE
of selected τ in global search (open circle) is smallest among those ofSimulation tests for Evaluation
Fig. 5 shows the results of the simulation tests for Evaluation. Fig.
5A shows the time course of the means and SDs of the residual errors
between the original and reconstructed simulated data. When thepreset number of waveforms K is smaller than true (K = 1,2), theresidual errors are nonstationary; the time courses of the means and/Fig. 5. Simulation tests for Evaluation. (A) Means (solid lines) and means ± SDs (dotted lines) across trials of residual errors between original and reconstructed simulated data. (B)Correct rate (%) of estimating true number of waveforms K for each SNR.

Y. Takeda et al. / NeuroImage 51 (2010) 629–641

635

Fig. 6. SDs across trials of residual errors obtained from EEGs at Fz, Cz, and Pz. (A) Go trials. (B) NoGo trials. Red lines represent SDs across trials of EEGs during Go/NoGo task. Blue,green, and black lines respectively represent SDs across trials of residual errors for preset K of 1, 2, and 3. Time 0 corresponds to stimulus onsets.

or SDs of the residual errors are not constant. In contrast, when preset
K is true (K = 3), the residual errors are stationary; the time courses ofthe means and SDs of the residual errors are constant. These results
indicate the rationale of the Evaluation procedure. Fig. 5B shows therelations between the correct rates (%) of estimating true K and SNR.
The correct rates become higher as SNR becomes higher, indicatingthat the reliability of Evaluation becomes higher as SNR becomes
higher.
Application to EEGs during Go trials
Figs. 6A, 7, 8, and 9 show the results of applying the proposed
method to the EEGs during the Go trials. Fig. 6A shows the SDs acrossthe trials of the residual errors between the original and reconstructed

EEGs. When 3 is the preset number of waveforms, the time courses ofthe means and the SDs of the residual errors become almost constantand the number of waveforms is estimated to be 3. Figs. 7, 8, and 9A
show the estimated waveforms at Fz, Cz, and Pz. The correlation
coefﬁcient between the estimated waveform-1 at Cz and thestimulus-triggered average EEG at Cz is 0.67, indicating that both
waveforms are similar (Fig. 7A, middle). Waveform-1s exhibit P300 as
well as the stimulus-triggered average EEGs (Fig. 7A). Waveform-1s
have large oscillatory components around 700 ms after the stimulus
onsets (Fig. 7A). The correlation coefﬁcient between the estimatedwaveform-2 at Cz and the response-triggered average EEG at Cz is
0.54, indicating that both waveforms are similar (Fig. 8A, middle).
Waveform-2s have large oscillatory components around 400 ms after
the response onsets. Figs. 7, 8, and 9B show the histograms of theFig. 7. Estimated waveform-1s from EEGs during Go trials. (A) Estimated waveform-1s at Fz, Cz, and Pz (black lines). Red lines represent stimulus-triggered average EEGs at Fz, Cz,
and Pz. Time 0 corresponds to stimulus onset. (B) Histograms of delays of estimated waveform-1s. Time 0 corresponds to stimulus onset. (C) Scalp distributions of variance acrosstime of estimated waveform-1s. (D) Average amplitude spectra of estimated waveform-1s at Fz, Cz, and Pz (black line). Red line represents average amplitude spectra of stimulustriggered average EEGs at Fz, Cz, and Pz.

636

Y. Takeda et al. / NeuroImage 51 (2010) 629–641

Fig. 8. Estimated waveform-2s from EEGs during Go trials. (A) Estimated waveform-2s at Fz, Cz, and Pz (black lines). Red lines represent response-triggered average EEGs at Fz, Cz,
and Pz. Time 0 corresponds to button-push signal onset. (B) Histograms of delays of estimated waveform-2s. Time 0 corresponds to stimulus onset. (C) Scalp distributions of varianceacross time of estimated waveform-2s. (D) Average amplitude spectra of estimated waveform-2s at Fz, Cz, and Pz (black line). Red line represents average amplitude spectra of
response-triggered average EEGs at Fz, Cz, and Pz.

delays of the estimated waveforms. The estimated delays are 18.54 ±14.87 ms for waveform-1s, 276.25 ± 47.95 ms for waveform-2s, and329.48 ± 125.58 ms for waveform-3s. Figs. 7, 8, and 9C show thescalp distributions of the variance across time of the estimatedwaveforms. All have large variance around Cz. Waveform-1s also
have large variance around O1 and O2. Figs. 7, 8, and 9D show theaverage amplitude spectra of the estimated waveforms at Fz, Cz,
and Pz. All have large amplitude spectra around 2–3 Hz, andwaveform-3s also have large amplitude spectra around 4 Hz. Wave

form-1s and waveform-2s have larger amplitude spectra around
10 Hz than the stimulus- and response-triggered average EEGs,
respectively.
Application to EEGs during NoGo trials
Figs. 6B, 10, and 11 show the results of applying the proposed
method to the EEGs during the NoGo trials. Fig. 6B shows the SDs
across the trials of the residual errors between the original andFig. 9. Estimated waveform-3s from EEGs during Go trials. (A) Estimated waveform-3s at Fz, Cz, and Pz. Horizontal axes represent relative time to deﬁned onsets of estimatedwaveform-3s. (B) Histograms of delays of estimated waveform-3s. Time 0 corresponds to stimulus onset. (C) Scalp distributions of variance across time of estimated waveform-3s.
(D) Average amplitude spectra of estimated waveform-3s at Fz, Cz, and Pz.

Y. Takeda et al. / NeuroImage 51 (2010) 629–641

637

Fig. 10. Estimated waveform-1s from EEGs during NoGo trials. (A) Estimated waveform-1s at Fz, Cz, and Pz (black lines). Red lines represent stimulus-triggered average EEGs at Fz,
Cz, and Pz. Time 0 corresponds to stimulus onset. (B) Histograms of delays of estimated waveform-1s. Time 0 corresponds to stimulus onset. (C) Scalp distributions of variance acrosstime of estimated waveform-1s. (D) Average amplitude spectra of estimated waveform-1s at Fz, Cz, and Pz (black line). Red line represents average amplitude spectra of stimulustriggered average EEGs at Fz, Cz, and Pz.

reconstructed EEGs. When 2 is the preset number of waveforms, thetime courses of the means and the SDs of the residual errors
become almost constant and the number of waveforms is estimatedto be 2. Figs. 10 and 11A show the estimated waveforms at Fz, Cz,
and Pz. The correlation coefﬁcient between the estimated waveform-1 at Cz and the stimulus-triggered average EEG at Cz is 0.95,
indicating that both waveforms are similar (Fig. 10A, middle).
Waveform-1s exhibit N200 and P300 as well as the stimulustriggered average EEGs (Fig. 10A). Figs. 10 and 11B show thehistograms of the delays of the estimated waveforms. The estimateddelays are 18.10 ± 12.53 ms for waveform-1s, and 336.30 ± 109.61
ms for waveform-2s. Figs. 10 and 11C show the scalp distributionsof the variance across time of the estimated waveforms. All havelarge variance around Cz. Waveform-1s also have large variancearound O1 and O2. Figs. 10 and 11D show the average amplitude
spectra of the estimated waveforms at Fz, Cz, and Pz. All have large
amplitude spectra around 2–3 Hz, and waveform-2s also have large
amplitude spectra around 5 Hz.

Fig. 11. Estimated waveform-2s from EEGs during NoGo trials. (A) Estimated waveform-2s at Fz, Cz, and Pz. Horizontal axes represent relative time to deﬁned onsets of estimatedwaveform-2s. (B) Histograms of delays of estimated waveform-2s. Time 0 corresponds to stimulus onset. (C) Scalp distributions of variance across time of estimated waveform-2s.
(D) Average amplitude spectra of estimated waveform-2s at Fz, Cz, and Pz.

638

Y. Takeda et al. / NeuroImage 51 (2010) 629–641

Discussion
In this study, we proposed a generalized method to estimate EEG
waveforms common across trials. From single/multi-channel EEGs,
the method estimates the number of waveforms common acrosstrials, the delays of waveforms, and all of the waveforms. The
performance of the algorithm was veriﬁed by a number of simulation
tests. We also applied this method to EEGs during a Go/NoGo task.
The purpose of the proposed method, estimating all EEG waveforms common across trials, is much more demanding than the tasksof our previous methods (Takeda et al., 2008a; Takeda et al., 2008b)
and the other methods for estimating a waveform whose delays arevariable and unknown (Biggins et al., 1997; Jaśkowski and Verleger,
1999, 2000; McGillem et al., 1985; Möcks et al., 1988; Pham et al.,
1987; Puce et al., 1994a,b; Woody, 1967). We achieved such ademanding goal with a step-by-step approach. We divided our
purpose into three easier subproblems (estimate waveforms, estimate
their delays, and estimate their number) and solved the subproblems
step-by-step, as shown in Fig. 1. This strategy is an essential point in
the versatility of the proposed method.
Another method for separating EEGs is spatial decomposition by
independent component analysis (ICA) (Jung et al., 2001; Makeig etal., 2004). ICA is a computational method for decomposing multichannel data into mutually independent components with different
scalp distributions. Therefore, to separate waveforms by ICA, the EEGs
need to satisfy the following three conditions: (1) the number ofchannels is larger than that of the waveforms, (2) the variability of thedelays of waveforms is adequately large, and (3) the waveforms havedifferent scalp distributions. Condition 1 is needed because ICA cannot
distinguish a number of signals larger than the channels. Condition 2
is needed because ICA cannot distinguish correlated signals. When the
delays of waveforms are within a small range, a waveform tends tooverlap on certain phases of the other waveforms, and their timeseries tend to be correlated with each other. In fact, in our simulation
test, ICA did not separate waveforms when the variability of the delays
was small (data not shown). Condition 3 is needed because ICA cannot
distinguish signals attributable to identical sources. In fact, in our
simulation test, ICA did not separate waveforms when the source of
the waveforms was the same (data not shown). We decomposed theEEGs during the Go/NoGo trials into independent components by ICA,
and applied the proposed method to the components. As a result, two
or more waveforms were estimated from some of the components(data not shown), indicating that the EEGs during the Go/NoGo trials
did not satisfy the three conditions. In contrast to ICA, the proposed
method does not need these three conditions. Therefore, from theviewpoint of estimating waveforms common across trials, weconsider our method better than ICA.
The proposed method consists of three steps: Delay estimation,
Waveform estimation, and Evaluation. Below, we discuss the threesteps individually.
Waveform estimation
In the Waveform estimation step, we estimate waveforms
common across trials when the number and the delays of waveforms
are given. The validity of the procedure is veriﬁed by the simulation
results (Fig. 3).
When the number of waveforms is 2, the purpose of Waveform
estimation is almost the same as that of our previous method (Takeda
et al., 2008a). However, their performances are different. A comparison of error coefﬁcients a of Waveform estimation and our previous
method (Takeda et al., 2008a) indicates the superiority of the Waveform estimation procedure. This is due to the differences between theirprocedures. In our previous method (Takeda et al., 2008a), waveforms(ch)were estimated by algebraically solving equation y(ch)n (t) = s1 (t)(ch)+ s(ch)
(t-τ
)
+
v
(t)
with
some
averagingprocedures.
Asaresult,2nn

our previous method has a limitation: slow waves (∼ 1 Hz) in noise
are ampliﬁed by the estimation. In the Waveform estimation
procedure, waveforms are estimated by the least squares method,which somehow eliminates the limitation. Therefore, it appears best
to adopt the Waveform estimation procedure rather than our previous
method (Takeda et al., 2008a), even when the number of waveforms
is known to be 2.
Delay estimation
In the Delay estimation step, we estimate the delays of waveforms
common across trials by solving the optimization problem of Eq. (7).
The validity of Eq. (7) and the procedure for solving it are veriﬁed by
the simulation results (Fig. 4A).
The simulation results in Fig. 4B show that estimation accuracy is
low when SNR is low. This is because, when SNR is low, wrong delays
tend to minimize the value of objective function oτ more than true
delays. Consequently, the validity of Eq. (7) becomes low when SNR islow. On the other hand, the simulation results in Fig. 4C indicate that,when noise is independent across channels, estimation accuracy
increases as the number of meaningful channels increases. Thissuggests a solution for improving estimation accuracy when SNR is
low: adding meaningful channels.
The simulation results in Fig. 4B also show that estimation accuracy
becomes lower when the number of waveforms increases. This is dueto the increased difﬁculty of optimization. As the number of waveforms increases, the number of delays to be searched for increases andthe difﬁculty of optimization increases. Therefore, when there seems
to be many waveforms whose delays are unknown, we need to spend
much time for optimization or to ﬁnd better optimization algorithms.
Evaluation
In the Evaluation step, we evaluate whether the preset number ofwaveforms common across trials is true by examining the residual
errors between original and reconstructed EEGs. The rationale of theprocedure for Evaluation is veriﬁed by the simulation results (Fig. 5A).
To select a criterion for evaluating the preset number of waveforms, by using simulated data (not shown), we tested various
criteria, such as Akaike's information criterion (AIC) (Akaike, 1974)
and the cross-validation method. Among the criteria we tested, theprocedure used in this study is the best from the viewpoint ofproviding stable and reasonable performance (Fig. 5B).
The simulation results in Fig. 5B show that the reliability of Evaluation is low when SNR is low. This is because, as SNR becomes lower,
residual errors become more noisy, and testing whether the
distributions of residual errors before and after stimulus onsets are
the same becomes more difﬁcult. Therefore, to increase the reliabilityof Evaluation, we need to increase SNR.
It is possible that the numbers of waveforms are different acrosschannels. For example, EEGs at distant electrodes may have a different
number of waveforms, e.g., 1 waveform for electrode-1 and 2
waveforms for electrode-2. In such cases, the proposed method
estimates the maximum number of waveforms, and some waveforms
estimated from EEGs that have the smaller number of waveforms
become ﬂat. In the above example, the estimated number ofwaveforms should be 2, and one of the waveforms estimated from
EEGs at electrode-1 should be ﬂat. Therefore, we can know theeffective number of waveforms for each channel by examining
estimated waveforms.
Some extended usages
We have described a basic usage of our method to estimate EEG
waveforms common across trials. In practice, we can use this method
in a variety of ways depending on our needs and situations. For

Y. Takeda et al. / NeuroImage 51 (2010) 629–641

example, only the Waveform estimation step is needed when we
estimate the approximate waveforms of stimulus- and responselocked components (Braun et al., 2002; Endo et al., 1999; Goodin et al.,
1986; Jung et al., 2001; Makeig et al., 2004) from EEGs during
stimulus–response tasks. Further, we can easily extend the method.
We describe some extended usages of the methods below.
Using a priori knowledge about importance of channels and frequencies
Sometimes, we have a priori knowledge about the importance ofchannels and/or frequencies. For example, EEGs at frontal channels
are sometimes contaminated with EOG and have low SNR. Further,EEGs at high frequencies (N50 Hz) are sometimes contaminated withelectromyographic activity. In such cases, using these channels andfrequencies as well as others may decrease the accuracy of theestimation, as shown in Fig. 4C. In the Delay estimation step, we canuse the knowledge about the importance of channels and frequencies
by replacing Eq. (7) withτ̂ = argmin
τCHXch = 1

W ch ðchÞ

T =2
X

ω=1

W ω ðωÞj jYðchÞðωÞ−Eðω; τÞhi−1TT ðchÞ
2
Eðω; τÞ Y ðωÞj j ;× Eðω; τÞ Eðω; τÞ

ð9Þ

where Wch(ch): a weight function of ch, and Wω(ω): a weightfunction of ω.
Using a priori knowledge about delays of waveforms
Sometimes we have a priori knowledge about the delays of EEG
waveforms common across trials. For example, in stimulus–responsetasks, the approximate delays of two waveforms can be given fromstimulus and response onsets. In such cases, using the given delays
simpliﬁes the optimization problem [Eq. (7)] and may increase
estimation accuracy. In the Delay estimation step, we can use theknowledge about the delays by restricting the delays' search space tothe neighborhood of the given delays.
Using a priori knowledge about waveforms
Sometimes we have a priori knowledge about EEG waveforms
common across trials. For example, an approximate waveform can begiven by averaging EEGs triggered on stimulus onsets. In such cases,
using known waveforms simpliﬁes the optimization problem [Eq. (7)]and may increase the estimation accuracy. In the Delay estimation step,
we can use the knowledge about waveforms by replacing Eq. (7) with
T =2
CH X
X
ðchÞ
ðchÞ
τ̂ = argmin
j jY ðωÞ − E1 ðω; τ 1 ÞS1 ðωÞ−E2 ðω; τ 2 Þτ
ch = 1 ω = 1
h
i−1
!
"
T
T
ðchÞ
ðchÞ
2
× E2 ðω; τ 2 Þ E2 ðω; τ 2 Þ E2 ðω; τ2 Þ Y ðωÞ−E1ðω; τ1 ÞS1 ðωÞ jj ;ð10Þwhere S1(ch)(ω): a L-by-1 matrix generated from known L waveforms, E1(ω,τ1): a N-by-L matrix generated from τ1 [=τn,k(n=1,⋯,N;k =1,⋯,L)],and E2(ω,τ2): a N-by-(K-L) matrix generated from τ2 [=τn,k(n=1,⋯,N;
k=L+1,⋯,K)].
Using waveform correlation across channels
In the case of low spatial resolution data, such as EEGs, waveforms
are correlated across channels. For example, waveforms at O1 may
resemble those at O2 because the two electrodes are close to eachother. In such cases, using the waveform correlation may increase theestimation accuracy. There are two ways to use the waveform
correlation: (1) by using spatial decomposition techniques before
applying our method, and (2) by assuming a model that incorporates
the waveform correlation. In the ﬁrst case, we ﬁrst decompose EEGs
by spatial decomposition techniques, such as ICA or principal
component analysis (PCA). Then, we apply the proposed method to639

the decomposed components. Because the decomposition techniques
separate signals from noises, the estimation accuracy increases. In thesecond case, we assume, for example, a model in which temporal
waveforms are the same but their amplitudes are different acrosschannels. In this model, an EEG can be expressed byðchÞyn ðt Þ =

K
X

k=1

!
"
ðchÞ
ak ðchÞsk t − τ n;k + vn ðt Þ;

ð11Þ

where y(ch)n (t): observed EEG epoch of channel ch in trial n, ak(ch):amplitude of k-th waveform of channel ch, sk(t): k-th waveform, τn,k:delay of sk(t) in trial n, v(ch)n (t): noise of channel ch in trial n, and K:
number of waveforms.
The unknown parameters ak(ch), sk(t), τn,k, and K can be estimatedin an iterative way. When Eq. (11) is valid, its estimation accuracy
would be higher than that of Eq. (1) because Eq. (11) has fewer
parameters than Eq. (1). In fact, in our preliminary simulation test, theestimation accuracy of Eq. (11) was higher (data not shown).
Target data
We focused on EEGs, but the proposed method can also be applied
to other kinds of brain imaging data, such as magnetoencephalography (MEG) data. Also, the method can be applied to preprocessed
EEG/MEG data. For example, the method can be applied to timefrequency data, such as a scalogram obtained by taking a wavelet
transform of EEG/MEG. In this case, y(ch)n (t) in Eq. (1) is regarded as
the value of time-frequency data at time t and frequency ch.
Limitation
Although the proposed method seems generally useful for wide
EEG analyses, it also has an inherent limitation: the validity of theassumption [Eq. (1)]. The proposed method assumes that noise is astationary process, and in the Evaluation step, nonstationary residual
errors are due to the remaining waveforms common across trials.
However, it is possible that the nonstationary residual errors are due
to other factors. Nonstationary background noise or the variability of
waveforms (Mihaylova et al., 1999; Vassilev et al., 2002; Vaughan etal., 1966) may be responsible for the nonstationary residual errors. In
such cases, the proposed method may extract nonexistent false
waveforms. To prevent this, before applying the proposed method, weneed to examine whether the nonstationary residual errors are
actually due to waveforms common across trials.
Applications to EEGs during Go/NoGo task
Before applying the proposed method to the EEGs during the Go/
NoGo task, we examined whether nonstationary EEGs are due to otherfactors than waveforms common across trials. We focused on the SDsof the EEGs across the trials and examined whether the increased SDs
after the stimulus onsets (Fig. 6, red lines) are due to other factors
than the variable delays of waveforms common across trials.
We examined whether the variability in the amplitudes ofwaveforms was responsible for the increased SDs. From the EEGs
during the Go trials, we estimated stimulus- and response-locked
waveforms by the Waveform estimation procedure using the RTs,
estimated the trial-to-trial variability of the amplitudes of theestimated waveforms by the least squares method, and obtained
residual errors. As a result, increases in the residual errors' SDs still
occurred (data not shown). From the EEGs during the NoGo trials, weestimated stimulus-locked waveforms by the stimulus-triggered
averaging procedure, estimated the trial-to-trial variability of theamplitudes of the estimated waveforms by the least squares method,and obtained the residual errors. As a result, increases in the residual
errors' SDs still occurred (data not shown). Therefore, we concluded
that the variability of the amplitudes of the waveforms was not fully

640

Y. Takeda et al. / NeuroImage 51 (2010) 629–641

responsible for the increased SDs during the Go/NoGo task. Furthermore, we examined whether the stimulus increased the amplitude of
the background noise and whether the increased background noise
was responsible for the increased SDs. We examined the SDs of theEEGs during the passive viewing task. Since the SDs did not show suchdrastic increases as the EEGs during the Go/NoGo task (see Takeda etal., 2008b), we concluded that the increased background noise by thestimulus was not fully responsible for the increased SDs during the
Go/NoGo task. Based on these preliminary examinations, we assumed
the validity of the assumption [Eq. (1)] for the EEGs during the Go/
NoGo tasks and thus applied the proposed method to the EEGs.
From the EEGs during the Go/NoGo tasks, for the ﬁrst time weestimated the numbers of waveforms common across trials, theirdelays, and all of the waveforms. As the preset number of waveforms
increases, the time courses of the SDs of the residual errors become
more constant (Fig. 6). This suggests that the estimated waveforms
and their delays are responsible for the increased SDs. The estimatedwaveforms are discussed below.
Waveforms time-locked to stimulus onsets
Waveform-1s in the Go and NoGo trials are stimulus-locked.
Therefore, waveform-1s are considered to reﬂect stimulus-related
brain processes, such as perception of the visual stimuli. This is
conﬁrmed by the scalp distributions of waveform-1s, which havelarge power at the occipital regions (Figs. 7 and 10C).
The early parts (0–300 ms) of waveform-1s in the Go trials
resemble those of the stimulus-triggered average EEGs during the Go
trials (Fig. 7A). This indicates that, in the early parts (0–300 ms) of thestimulus-triggered average EEGs during the Go trials, the effect of theoverlapping of the other waveforms, such as movement-related
potentials (MRPs), is small. The effect of the overlapping of MRPs on
stimulus-triggered average EEGs during Go trials has been debated(Kok, 1988; Smith et al., 2008; Verleger, 1988; Verleger et al., 2006).
With the proposed method, we extracted ourselves from thatproblem because we can obtain pure waveforms uncontaminatedwith MRPs.
Waveform-1s in the NoGo trials resemble the stimulus-triggered
average EEGs during the NoGo trials (Fig. 10A). This indicates that, inthe stimulus-triggered average EEGs during the NoGo trials, the effectof the overlapping of the other waveforms is small. Waveform-1s in
the NoGo trials exhibit N200 and P300 as well as the stimulustriggered average EEGs. The peaks in the waveform-1s in the NoGo
trials may reﬂect the motor inhibition or the detection of the response
conﬂict as discussed for those in the Go/NoGo literature (Bokura et al.,
2001; Donkers and van Boxtel, 2004; Falkenstein et al., 1999;
Ramautar et al., 2004).
Waveforms time-locked to response onsets
Waveform-2s in the Go trials are response-locked. Therefore, thesewaveforms are considered to reﬂect such response-related brain
processes as execution of button pushes.
Around 400 ms after the response onsets, waveform-2s in the Go
trials have large alpha band oscillations (Fig. 8A, black lines). This fact
is conﬁrmed by the amplitude spectra (Fig. 8D, black line) and thetime-frequency plots (data not shown). The alpha oscillation after the
response onsets may be regarded as the event-related synchronization (ERS) of the alpha band (Neuper and Pfurtscheller, 2001). The
ERS of the alpha band may represent a deactivated cortical area or
inhibited cortical networks (Neuper and Pfurtscheller, 2001). In
contrast to the waveform-2s in the Go trials, the response-triggered
average EEGs during the Go trials do not have the ERS of the alphaband (Fig. 8A, red lines). This may indicate that the proposed method
extracts the ERS of the alpha band but the response-triggered
averaging procedure does not.

Waveforms time-locked to neither stimulus nor response onsets
Waveform-3s in the Go trials and waveform-2s in the NoGo trials
are time-locked to neither the stimulus nor the response onsets. Such
waveforms have been previously hard to see because of their
unknown delays.
Verleger et al. (2005) examined EEGs during choice reaction time
tasks by stimulus- and response-triggered averaging procedures andshowed evidences that a component (P3b) is time-locked to neither
stimulus nor motor response onsets. Like P3b, the positive peaks inthe waveform-3s in the Go trials appear around 300 ms after the
stimulus onsets and are time-locked to neither the stimulus nor the
response onsets. Verleger et al. (2005) suggested that P3b reﬂects aprocess that mediates between perceptual analysis and response
initiation. It remains unclear whether their interpretation can beapplied to the positive peaks in the waveform-3s in the Go trials
because the delays of the positive peaks in the waveform-3s are
longer than the RTs in about half of the Go trials; if the peaks reﬂect
the mediating process, their delays should be between the stimulus
and the response onsets.
On the other hand, the properties (temporal waveforms, delays,
scalp distributions, and amplitude spectra) of waveform-3s in the Go
trials resemble those of waveform-2s in the NoGo trials (Figs. 9 and
11). This may indicate that these waveforms reﬂect the same brain
process in the Go and NoGo trials, such as monitoring task
performance. To elucidate the functional role of the waveforms, weneed further examinations, such as a comparison with the waveforms
estimated from EEGs during a variety of tasks and a source estimationof the waveforms.
Conclusion
We proposed a generalized method to estimate EEG waveforms
common across trials. The main achievement of the proposed method
is its ability to deal with an unknown number of multiple waveforms
and multi-channel EEGs. In general situations, this method can beused in a variety of ways.
Acknowledgments
This research was supported by a contract with the National
Institute of Information and Communications Technology entitled,“Multimodal integration for brain imaging measurements”.
Appendix A
The procedure for random search in global search is described in
MATLAB style as follows:
Generate τ by random numbers;
Obtain oτ;for iter = 1:20
for k = 1:K
for n = 1:N
Make τ′ by changing τn,k in τ randomly;
Obtain oτ′;
if oτ′ b oτoτ = oτ′;τ = τ′;
end
end
end
end.

Y. Takeda et al. / NeuroImage 51 (2010) 629–641

The procedure for grid search in local search is described in
MATLAB style as follows:Set τ selected in global search;
Obtain oτ;o′ = oτ + 1;
while oτ b o′o′ = oτ;for k = 1:K
for n = 1:N
for t = τmin:τstep:τmax
Make τ′ by changing τn,k in τ to t;
Obtain oτ′;
if oτ′ b oτoτ = oτ′;τ = τ′;
end
end
end
end
end,
where τmin: minimum value of τn,k, τstep: step of τn,k, and τmax:
maximum value of τn,k.
References
Akaike, H., 1974. A new look at the statistical model identiﬁcation. IEEE Trans. Automat.
Contr. AC-19 716–723.
Biggins, C.A., MacKay, S., Clark, W., Fein, G., 1997. Event-related potential evidence forfrontal cortex effects of chronic cocaine dependence. Biol. Psychiatry 42, 472–485.
Bokura, H., Yamaguchi, S., Kobayashi, S., 2001. Electrophysiological correlates forresponse inhibition in a Go/NoGo task. Clin. Neurophysiol. 112, 2224–2232.
Braun, C.M.J., Villeneuve, L., Gruzelier, J.H., 2002. Topographical analysis of stimulusrelated and response-related electrical scalp activity and interhemispheric
dynamics in normal humans. Int. J. Psychophysiol. 46, 109–122.
Donkers, F.C.L., van Boxtel, G.J.M., 2004. The N2 in go/no-go tasks reﬂects conﬂict
monitoring not response inhibition. Brain Cogn. 56, 165–176.
Endo, H., Kizuka, T., Masuda, T., Takeda, T., 1999. Automatic activation in the humanprimary motor cortex synchronized with movement preparation. Cogn. Brain Res.
3, 229–239.
Falkenstein, M., Hoormann, J., Hohnsbein, J., 1999. ERP components in Go/Nogo tasks
and their relation to inhibition. Acta Psychol. 101, 267–291.
Goodin, D.S., Aminoff, M.J., Mantle, M.M., 1986. Subclasses of event-related potentials:
response-locked and stimulus-locked components. Ann. Neurol. 20, 603–609.
Jaśkowski, P., Verleger, R., 1999. Amplitudes and latencies of single-trial ERP's
estimated by a maximum-likelihood method. IEEE Trans. Biomed. Eng. 46,
987–993.

641

Jaśkowski, P., Verleger, R., 2000. An evaluation of methods for single-trial estimation ofP3 latency. Psychophysiol 37, 153–162.
Jung, T.-P., Makeig, S., Westerﬁeld, M., Townsend, J., Courchesne, E., Sejnowski, T.J.,
2001. Analysis and visualization of single-trial event-related potentials. Hum. Brain
Mapp. 14, 166–185.
Kok, A., 1988. Overlap between P300 and movement-related-potentials: a response to
Verleger. Biol. Psychol. 27, 51–58.
Makeig, S., Delorme, A., Westerﬁeld, M., Jung, T.-P., Townsend, J., Courchesne, E.,
Sejnowski, T.J., 2004. Electroencephalographic brain dynamics following manually
responded visual targets. PLoS Biol. 2, 0747–0762.
McGillem, C.D., Aunon, J.I., Pomalaza, C.A., 1985. Improved waveform estimation
procedures for event-related potentials. IEEE Trans. Biomed. Eng. 32, 371–379.
Mihaylova, M., Stomonyakov, V., Vassilev, A., 1999. Peripheral and central delay inprocessing high spatial frequencies: reaction time and VEP latency studies. Vision
Res. 39, 699–705.
Möcks, J., Köhler, W., Gasser, T., Pham, D.T., 1988. Novel approaches to the problem oflatency jitter. Psychophysiol 25, 217–226.
Neuper, C., Pfurtscheller, G., 2001. Event-related dynamics of cortical rhythms:
frequency-speciﬁc features and functional correlates. Int. J. Psychophysiol. 43,
41–58.
Pham, D.T., Möcks, J., Köhler, W., Gasser, T., 1987. Variable latencies of noisy signals:
estimation and testing in brain potential data. Biometrika 74, 525–533.
Puce, A., Berkovic, S.F., Cadusch, P.J., Bladin, P.F., 1994a. P3 latency jitter assessed using 2techniques: I. Simulated data and surface recordings in normal subjects.
Electroencephalogr. Clin. Neurophysiol. 92, 352–364.
Puce, A., Berkovic, S.F., Cadusch, P.J., Bladin, P.F., 1994b. P3 latency jitter assessed using 2techniques: II. Surface and sphenoidal recordings in subjects with focal epilepsy.
Electroencephalogr. Clin. Neurophysiol. 92, 555–567.
Ramautar, J.R., Kok, A., Ridderinkhof, K.R., 2004. Effects of stop-signal probability in thestop-signal paradigm: the N2/P3 complex further validated. Brain Cogn. 56,
234–252.
Smith, J.L., Johnstone, S.J., Barry, R.J., 2008. Movement-related potentials in the Go/
NoGo task: the P3 reﬂects both cognitive and motor inhibition. Clin. Neurophysiol.
119, 704–714.
Takeda, Y., Yamanaka, K., Yamamoto, Y., 2008a. Temporal decomposition of EEG during
a simple reaction time task into stimulus- and response-locked components.
NeuroImage 39, 742–754.
Takeda, Y., Yamanaka, K., Nozaki, D., Yamamoto, Y., 2008b. Extracting a stimulusunlocked component from EEG during NoGo trials of a Go/NoGo task. NeuroImage41, 777–788.
Tallon-Baudry, C., Bertrand, O., 1999. Oscillatory gamma activity in humans and its rolein object representation. Trends Cogn. Sci. 3, 151–162.
Vassilev, A., Mihaylova, M., Bonnet, C., 2002. On the delay in processing high spatial
frequency visual information: reaction time and VEP latency study of the effect of
local intensity of stimulation. Vision Res. 42, 851–864.
Vaughan Jr., H.G., Costa, L.D., Gilden, L., 1966. The functional relation of visual evoked
response and reaction time to stimulus intensity. Vision Res. 6, 645–656.
Verleger, R., 1988. The true P3 is hard to see: some comments on Kok's (1986) paper on
degraded stimuli. Biol. Psychol. 27, 45–50.
Verleger, R., Jaśkowski, P., Wascher, E., 2005. Evidence for an integrative role of P3b inlinking reaction to perception. J. Psychophysiol. 19, 165–181.
Verleger, R., Paehge, T., Kolev, V., Yordanova, J., Jaśkowski, P., 2006. On the relation of
movement-related potentials to the go/no-go effect on P3. Biol. Psychol. 73,
298–313.
Woody, C.D., 1967. Characterization of an adaptive ﬁlter for the analysis of variable
latency neuroelectric signals. Med. Biol. Eng. 5, 539–553.
Zhigljavsky, A.A., 1991. Theory of Global Random Search. Kluwer Academic Publishers.