IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 46, NO. 5, MAY 1999

515

MEG Spatio-Temporal Analysis Using a CovarianceMatrix Calculated from Nonaveraged
Multiple-Epoch Data
Kensuke Sekihara,* Member, IEEE, David Poeppel, Alec Marantz, Hideaki Koizumi, and Yasushi Miyashita

Abstract— We propose a magnetoencephalographic (MEG)
spatio-temporal analysis in which the measurement-covariancematrix is calculated using nonaveraged multiple epoch data.
The proposed analysis has two advantages. First, a very narrow
time window can be used for the source estimation. Second,
accurate localization is possible even when the source activation
has a time jitter. Experiments using auditory evoked MEG dataclearly demonstrate these advantages.
Index Terms—Biomagnetics, biomedical electromagnetic imaging, biomedical signal processing, functional brain imaging, inverse problems.

I. INTRODUCTIONSPATIO-TEMPORAL modeling [1], [2] is known to reducethe degree of the ill-posedness in the biomagnetic inverse
problem [3]. In this modeling, a measurement-covariancematrix (or equivalently a spatio-temporal data matrix) is constructed, and source locations are estimated using this covariance matrix, assuming that the source locations are unchanged
during the measurement.
It is common in studies of evoked neuromagnetic fields
to measure many epochs by applying a stimulus repeatedly
and to analyze the data averaged over these epochs. When
applying the spatio-temporal analysis to such evoked field
data, the measurement-covariance matrix is calculated by using
this averaged data. This conventional way of calculating themeasurement-covariance matrix has two problems.
First, the time window for calculating the covariance matrix
must be large enough to contain a sufficient number oftime points to ensure the accuracy of the sample covariancematrix. The use of a large window, however, may result in
Manuscript received August 7, 1997; revised November 19, 1998. This
work was supported in part by the MIT-JST International Research Project“Mind Articulation.” Asterisk indicates corresponding author.
*K. Sekihara is with the Mind Articulation Project, Japan Science andTechnology Corporation (JST), Yushima 4-9-2, Bunkyo-ku, Tokyo 113-0034,
Japan (e-mail: sekihara@ma.jst.go.jp).
D. Poeppel is with the Department of Linguistics and Biology, Universityof Maryland, College Park, MD 20742 USA.
A. Marantz is with the Department of Linguistics and Philosophy, Massachusetts Institute of Technology, Cambridge, MA 02139 USA.
H. Koizumi is with the Mind Articulation Project, Japan Science andTechnology Corporation (JST), Tokyo 113-0034, Japan.
Y. Miyashita is with the Mind Articulation Project, Japan Science andTechnology Corporation (JST), Tokyo 113-0034, Japan. He is also with the
Department of Physiology, The University of Tokyo, School of Medicine,
Hongo, Tokyo 113-0033, Japan.
Publisher Item Identifier S 0018-9294(99)03114-6.

an erroneous source localization, if the source configuration
rapidly changes within this time window. Second, averaging
the multiple-epoch data obviously assumes that all the epochs
are time-locked to the stimulus. Thus, when a nontime-locked
source exists, the source estimation may also be erroneous.
In this paper, we propose to calculate the measurementcovariance matrix directly from nonaveraged multiple-epochdata. By calculating the measurement-covariance matrix fromnonaveraged data, the source localization can be carried out
using a small time window, which is advantageous when thesource configuration changes rapidly. Also, source localization
using this covariance matrix is not influenced by whether thesource activation is time-locked or not, and the location of asource whose activation onset has a time jitter in each epoch
can be accurately estimated.
One difficulty with the use of nonaveraged data is that thedata inevitably contains a large amount of background spontaneous activity that is usually averaged out. In this paper, thisdifficulty is circumvented by the use of a prewhitening technique that incorporates the noise covariance matrix obtainedusing the prestimulus parts of the nonaveraged data. Specifically, we use the previously proposed prewhitening MUSIC
algorithm [4], [5], as the source localization procedure. The
MUSIC algorithm has the distinct advantage of providing a
suboptimal solution by using only a three-dimensional search,regardless of the number of sources. The use of the MUSIC
algorithm, however, is not essential to implement the proposed
method. Any other source estimation method [6] can be used
together with the prewhitening technique.
Section II describes our proposed analysis, followed bya brief explanation of the prewhitening MUSIC algorithm.
Section III presents the results of our experiments whichshow the effectiveness of the proposed analysis. Throughout
this paper, plain italics indicate scalars, lower-case boldface
italics indicate vectors, and upper-case boldface italics indicate
matrices. The superscript indicates the matrix transpose. The
eigenvalues are numbered in decreasing order.
II. METHOD
A. Prewhitening MUSIC Algorithm
The prewhitening MUSIC algorithm was proposed to reduce the influence of the background activity by taking its

0018–9294/99$10.00  1999 IEEE

516

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 46, NO. 5, MAY 1999

covariance matrix into account [5]. It can provide accurate source estimation results even when a large amount ofbackground activity exists. Since the details of the prewhitening MUSIC algorithm have already been reported [5], it
is only briefly explained here for the reader’s convenience.
th
We define the magnetic-field data measured by the, and also define a vector
detector coil at time as
as a set of measured
is the total number of detector coils. The
data where
spherical homogeneous conductor [7] is assumed, so we
consider two tangential components, the and components,of the source moment and assume the radial component
is zero. For simplicity, we assume that all sources have
fixed orientations during measurement. The lead field vecand
components of the source ataretors for theanddefined as. We define the lead field matrix.
for the source at as
Let us denote the measurement-covariance matrix as , and. Let us define as an eigenthe noise covariance matrix asvector obtained by solving the generalized eigenvalue problem,
, and
as
where
is the number of sources. The prewhitening MUSIC algorithm
makes use of the fact that the noise subspace projector
andare orthogonal at the true source locations,represents the normal vector of thewheresource orientation that gives the maximum orthogonality. Thisorthogonality can be checked by calculating the following
localization function:(1)
indicates the minimum generalized eigenwhere
value of the matrix pair in the parenthesis. The prewhitening
MUSIC algorithm detects the source locations by calculating
this localizing function in a volume where sources can exist,
and by choosing each location where this localizing function
reaches a peak as one dipole source location. It is worth notingthat the localizing function for nonprewhitening MUSIC is [8]Let us decompose the nonaveraged measured magnetic field
into the four components:
, and
,
(3)
is the signal evoked-magnetic field that is timeHere,is the signal evoked-magnetic
locked to the stimulus,
field that is not time-locked to the stimulus and has a timeis the background field
jitter in each activation onset,
disturbance caused by external sources such as magnetocardiography (MCG) interference or spontaneous neuromagnetic
is the intrinsic sensor noise approximated by
activity, andthe white Gaussian noise with variance .
epochs are measured. We
We assume that a total of
epochs as.
denote the summation over these
Then,
and

the

,

averaged field
are defined as:

,

,

,
,

,
,

, and
. Here, we can assume thatand. Accordingly, we have(4)
The conventional way of obtaining the measurementcovariance matrix is to calculate the time average during acertain time period. Let us denote the measurement-covariance. The summation overmatrix obtained in this manner as
and the number of time points used intime is denoted as. Then, we havethis summation is denoted as
(5)

where

is the unit matrix andis obtained from
(6)

(2)
, and
is an eigenvector
whereobtained by solving the ordinary eigenvalue problem,. Here,is the number of sources;
may be different
because in nonprewhitening cases, sources for thefrombackground activity may need to be taken into account.
B. Proposed Spatio-Temporal Analysis
Using Nonaveraged Data
1) Conventional Covariance Matrix Calculation: Because
calculating the covariance matrix theoretically requires aninfinite number of measurements, some kind of approximationneeds to be introduced in practice. In this subsection, we
describe both the conventional and proposed methods ofcalculating the measurement-covariance matrix, explaining
the difference between them.

In deriving (5), we assume that
. If the nontime-locked evoked activity
does
is reduced to
not exist, the measurement-covariance matrix(7). Thus, in this case,where
can provide accurate information regarding the evoked activity
and accurate localizations of the sources forcan.
be obtained using thisHowever, two problems arise when utilizing (5). First, toobtain sufficient accuracy, a fairly large number of time points
are needed to calculate the sample covariance matrix. Second,cannot be ignored, thewhen the time-unlocked activity
may contain errors caused
source estimation obtained from.
by the activity

SEKIHARA et al.: MEG SPATIO-TEMPORAL ANALYSIS

517

2) Proposed Covariance Matrix Calculation: In this paper,
instead of using (5), we propose to calculate the measurementcovariance matrix using nonaveraged field data, such that(8)

andare the covariance matrices for evoked
whereactivities and for the background disturbance. These matrices
are obtained from
(a)

(9)
and
(10)
In deriving (8), we assume that the evoked activities areuncorrelated with the background activity. Equation (8) shows
contains information regardingthat the covariance matrix
as well as regarding
the signal source to be measured. However, if we obtain
the influence of the disturbance, which is equal to,
the noise covariance matrix
we can reduce the influence of the background disturbance byapplying a prewhitening technique. We show in the following
section that the previously proposed prewhitening MUSIC
algorithm [5] can almost completely remove this influence.

(b)
Fig. 1. Auditory evoked response used in the experiments: (a) one of the
256 raw epoch data and (b) the data averaged over 256 epochs. The data
portion from 300 to 500 ms is shown.

0III. EXPERIMENTS USING AUDITORY-EVOKED FIELD DATA
A. General Description
We performed source localization experiments using actual
auditory-evoked data to test the effectiveness of the proposed
analysis. The auditory-evoked fields elicited by a 1000-Hz puretone of 400-ms duration were measured using a 37-channel
Magnes magnetometer installed at the Biomagnetic Imaging
Laboratory, University of California, San Francisco. All measurements were done in a magnetically shielded room. The
stimulus was presented to the subject’s right ear, and the sensor
array was placed above the subject’s left hemisphere with theposition adjusted to optimally record the N1m auditory evoked
field.
Two-hundred-fifty-six epochs were measured with the average interstimulus interval of 2 s randomly varied between
1.75 and 2.25 s. The sampling frequency was set at 1 kHz,and the data was collected from 300 ms before the stimulus
onset to 700 ms after it. An on-line filter with the bandwidthfrom 1–400 Hz was used, and no additional digital filter wasapplied. One of the 256 epochs and the field data averaged over
all the epochs are shown in Fig. 1(a) and (b), respectively. In
Fig. 1(a), the single-epoch data contains a large amount of
spontaneous activity, but in (b), much of this activity was
suppressed through the averaging.
The , , and coordinates used to express the estimatedresults are depicted in Fig. 2. These coordinates are measured
wasin centimeters. The measurement-covariance matrix

Fig. 2. The x, y , and z coordinates used to express the estimated results.
The midpoint between the left and right preauricular points is defined as thecoordinate origin. The axis directed away from the origin toward the leftpreauricular point is defined as the y axis, and that from the origin to thenasion is the x axis. The z axis is defined as the axis that is perpendicular
to both these axes and directed from the origin to the vertex.

+

+

+

calculated, by using (5). The averaged field data shown in
Fig. 1(b) was used with the time window between 0 and
100 ms. The results of applying the nonprewhitening MUSIC
algorithm with this covariance matrix is shown in Fig. 3.
The MUSIC localizing function in (2) was calculated with,an interval of 0.5 cm over a volume defined as, and. The figure shows that thealgorithm detected the N1m source, which is believed to belocated in the primary auditory cortex area.
B. Detecting a Source Using a Small Time Window
with
We calculated the measurement-covariance matrix
the time window between 98 and 102 ms. This time window
was selected because it includes the N1m peak. Only six time

518

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 46, NO. 5, MAY 1999

Fig. 3. The results of plotting the nonprewhitening MUSIC localizer from(2) when applied to the averaged data shown in Fig. 1(b) with the time window
between 0 and 100 ms. The localizer was calculated with an interval of 0.5
7, and 3 z 11;
cm within a volume defined as 4 x 6, 0 y
the projections of the localizer values onto the transverse, coronal, and sagittal
planes are shown. The circles depicting a human head represent the projectionsof the sphere used to calculate the forward solution. L and R indicate thesubject’s left and right hemispheres.

0  

 

Fig. 4. The results of applying the nonprewhitening MUSIC algorithm tothe averaged data with a small time window between 98 and 102 ms.

 

points were used to calculate. The results of plotting thenonprewhitening MUSIC localization function (2) are shownin Fig. 4. In these results, the distortion is so severe that nosource can be detected due to the inaccuracy of the sample
covariance matrix.
using (8) withWe then calculated the covariance matrix
the same time window between 98 and 102 ms. The noise
covariance matrix
was calculated with the time window
between 300 and 0 ms. The nonprewhitening MUSIC algo. The results are shown
rithm was first applied using onlyin Fig. 5. In these results, the N1m source and another source
were detected; the second source was probably caused by the
spontaneous activity. The prewhitening MUSIC algorithm wasand, and the results are shown
then applied using bothin Fig. 6. In these results, only the single source located in the
primary auditory cortex area is clearly detected.
A comparison between Figs. 4 and 6 demonstrates the effectiveness of calculating the measurement-covariance matrix
with the nonaveraged multiple-epoch data when the time
window is very small. The comparison between Figs. 5 and6 shows the necessity of applying a prewhitening techniquewhen dealing with nonaveraged data.
Note that in theory the proposed covariance matrix is subject
to a stronger influence from the uncorrelated sensor noise thanthe conventional covariance matrix is. This can be understoodby noting that the sensor-noise influence is represented byin (5), whereas its influence is represented byin (8). However, when comparing Figs. 3 and 6, theFig. 5. The results of applying the nonprewhitening MUSIC algorithm to^ b was calculated by using (8)nonaveraged data. The covariance matrix R
with the time window between 98 and 102 ms.

greater influence of the sensor noise is not evident. This isprobably because the sensor noise in most modern biomagnetic
instruments, such as the one used in the current experiment, is
too low to influence the final source localization results, evenwhen nonaveraged data is used.
To further investigate this point, we check the eigenvalueand. The spectrum
spectra of the covariance matrices

SEKIHARA et al.: MEG SPATIO-TEMPORAL ANALYSIS

519

Fig. 8. The magnetic field averaged over 256 raw epoch data with anartificial time-jitter. The jitter was created with a uniform random number
distributed between 0 and 60 ms.

Fig. 6. The results of applying the prewhitening MUSIC algorithm tononaveraged data with the time window between 98 and 102 ms. The same
^ b as used to obtain the results in Fig. 5 was used. The
covariance matrix R^ n was calculated with the time window between
noise covariance matrix R
300 and 0 ms.

0in order to obtain the results shown in Fig. 3. In
Fig. 7(b), the noise-level eigenvalues are much larger than
, which was used to obtain thethose in Fig. 7(a) because
results in Fig. 6, contained the influence of the large sensor
noise from raw-epoch data. Nonetheless, we can still observe
two eigenvalues that are clearly larger than the noise-level
eigenvalues in this spectrum. Therefore, we were able toseparate the signal subspace from the noise subspace eventhough a strong influence of the large sensor noise wascontained in the covariance matrix. In Fig. 7(c), however, false
signal-level eigenvalues arose due to the inaccuracy in thesample covariance matrix, and as a result, the noise subspace
was inaccurately estimated.
C. Detecting a Source Whose Activation
Onset Has a Time Jitter

(a)

(b)

We next tested the effectiveness of the proposed analysis in
detecting a source whose activation onset varies from epoch
to epoch. To generate the data used for this experiment, we
artificially created a jitter for each epoch data. That is, each
was artificially time shifted to create
epoch measurement. Here, was a uniform random number generated in
a computer and it was distributed between 0 and 60 ms. The
field data averaged over all epochs were calculated using
(11)

(c)

~ b used to obtain the results shown in
Fig. 7. (a) Eigenvalue spectrum of R
^ b used to obtain the results shown in Fig. 6, and (c) that
Fig. 3, (b) that of R
~ b used to obtain the results shown in Fig. 4.
of R

in Fig. 7(a) relates to the results found in Fig. 3, the spectrum
in Fig. 7(b) relates to the results in Fig. 6, and the spectrum
in Fig. 7(c) relates to the results in Fig. 4. In Fig. 7(a), thefirst eigenvalue is distinctively large and the second one islarger than the other noise-level eigenvalues. Thus, we set

is shown in Fig. 8. ComThis averaged magnetic fieldparing this figure with Fig. 1(b), we can see that the peak at
approximately 50 ms of latency became less clearly defineddue to the artificially introduced time jitter . In addition, theappearance of the N1m peak was modified.
was calculated using this avThe covariance matrix
eraged field data with a time window between 30 and 130
ms. The results of the nonprewhitening MUSIC algorithm
are shown in Fig. 9. These results
obtained using thiscontain severe blurring and no accurate source estimation canbybe made. We then calculated the covariance matrix
using (8) with the same time window between 30 and 130 ms.
was calculated with the time
The noise covariance matrix
window between 240 and 0 ms. First, the nonprewhitening, andMUSIC algorithm (2) was performed only using
the results are shown in Fig. 10. Here, a source that was520

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 46, NO. 5, MAY 1999

Fig. 9. The results of applying the nonprewhitening MUSIC algorithm to the~ b was calculated
averaged field data shown in Fig. 8. The covariance matrix R
using (5) with the time window between 30 and 130 ms.

Fig. 11. The results of applying the prewhitening MUSIC algorithm tononaveraged data having the artificial time jitter. The same covariance matrix
R^ b as used to obtain the results in Fig. 10 was used. The noise covariance^ n was calculated with the time window between 240 and 0 ms.
matrix R

0

between Figs. 10 and 11 again demonstrates the necessity andeffectiveness of the prewhitening technique used here.
IV. CONCLUSION
A spatio-temporal analysis is proposed in which themeasurement-covariance matrix is calculated by using the
nonaveraged multiple-epoch evoked data. The proposed
analysis has two advantages. First, a very narrow time
window can be used. Second, the analysis can provide accurate
localization of a source whose activation onset has a time jitterin each epoch. Experiments using measured auditory evokeddata clearly showed that these advantages were obtained byusing the proposed covariance matrix.
ACKNOWLEDGMENT
The authors wish to thank S. Honma and Dr. T. Roberts fortheir help in performing auditory MEG measurements.
REFERENCES
Fig. 10. The results of applying the nonprewhitening MUSIC algorithm to^bnonaveraged data having the artificial time jitter. The covariance matrix R
was calculated using (8) with the time window between 30 and 130 ms.

probably caused by the spontaneous background activity was
detected, instead of the N1m source in the auditory cortex
area. Next, we applied the prewhitening MUSIC algorithm
and. The results areusing both covariance matrices
shown in Fig. 11. Here, the N1m source was clearly detected.
The comparison between Figs. 9 and 11 demonstrates theeffectiveness of the proposed analysis for detecting a sourcewith a time jitter at its activation onset. The comparison[1] M. Scherg and D. von Cramon, “Two bilateral sources of the late AEP
as identified by a spatio-temporal dipole model,” Electroencephalogra.
Clin. Neurophysiol., vol. 62, pp. 32–44, 1985.
[2], “A new interpretation of the generators of BAEP waves I–V:
Results of a spatio-temporal dipole model,” Electroencephalogra. Clin.
Neurophysiol., vol. 62, pp. 290–299, 1985.
[3] M. Hämäläainen, R. Hari, R. J. IImoniemi, J. Knuutila, and O. V.
Lounasmaa, “Magnetoencephalography-theory, instrumentation, and applications to noninvasive studies of the working human brain,” Rev.
Mod. Phys., vol. 65, pp. 413–497, 1993.
[4] K. Sekihara, S. Miyauchi, and H. Koizumi, “Covariance incorporated
MEG-MUSIC algorithm and its application to detect SI and SII whenlarge background brain activity exists,” Neuro. Image, vol. 3, no. 3, p.
92, June 1996.
[5] K. Sekihara, D. Poeppel, A. Marantz, H. Koizumi, and Y. Miyashita,“Noise covariance incorporated MEG-MUSIC algorithm: A method forSEKIHARA et al.: MEG SPATIO-TEMPORAL ANALYSIS

multiple-dipole estimation tolerant of the influence of background brain
activity,” IEEE Trans. Biomed. Eng., vol. 44, pp. 839–847, 1997.
[6] A. Achim, François Richer, and J.-M. Saint-Hilaire, “Methodologicalconsiderations for the evaluation of spatio-temporal source models,”
Electroencephalogra. Clin. Neurophysiol., vol. 79, pp. 227–240, 1991.
[7] J. Sarvas, “Basic mathematical and electromagnetic concepts of thebiomagnetic inverse problem,” Phys. Med. Biol., vol. 32, pp. 11–22,
1987.
[8] J. C. Mosher, P. S. Lewis, and R. M. Leahy, “Multiple dipole modelingand localization from spatio-temporal MEG data,” IEEE Trans. Biomed.
Eng., vol. 39, pp. 541–557, 1992.

Kensuke Sekihara (M’88) received the M.S. degreein 1976 and the Ph.D. degree in 1987, both from theTokyo Institute of Technology, Tokyo, Japan.
Since 1976, he has worked with Central ResearchLaboratory, Hitachi, Ltd., Tokyo, Japan. He wasa Visiting Research Scientist at Stanford University, Stanford, CA, from 1985 to 1986, and at
Basic Development, Siemens Medical Engineering,
Erlangen, Germany, from 1991 to 1992. He iscurrently working with the “Mind Articulation”
research project sponsored by the Massachusetts
Institute of Technology and the Japan Science and Technology Corporation.
His research interests include biomagnetic inverse problems and statisticalestimation theory, especially its application to noninvasive measurements ofbrain functions.
Dr. Sekihara is a member of the IEEE Medicine and Biology Society, theIEEE Signal Processing Society, and the International Society of Magnetic
Resonance in Medicine.

David Poeppel received the B.S. and Ph.D. degrees in cognitive neuroscienceat Massachusetts Institute of Technology (MIT), Cambridge, in 1990 and 1995,
respectively.
He is currently an Adjunct Assistant Professor in the Department ofRadiology, University of California, San Francisco, and Assistant Professor
in the Departments of Linguistics and Zoology at the University of Maryland
at College Park. His research uses the functional neuroimaging methods
magnetoencephalography (MEG) and functional magnetic resonance imaging(fMRI) to investigate the neural basis of speech and language processing.

Alec Marantz received the B.A. degree in psycholinguistics from Oberlin College, Oberlin, OH,
in 1978 and the Ph.D. degree in linguistics from
the Massachusetts Institute of Technology (MIT),
Cambridge, MA, in 1981.
He joined the faculty at MIT in 1990, where
he is currently a Professor of Linguistics in the
Department of Linguistics and Philosophy. His research interests include the syntax and morphology
of natural languages, linguistic universals, and theneurobiology of language. He is currently involved
in revising morphological theory within linguistics and in exploring MEG
techniques to uncover how the brain processes language.

521

Hideaki Koizumi was born in Tokyo, Japan, in1946. He graduated from the Department of Pureand Applied Sciences, the University of Tokyo,
Tokyo, Japan, in 1971. He received the doctoral
degree in physics from the University of Tokyoin 1976 for his thesis “Creating Polarized ZeemanEffect Atomic Absorption Spectrometry.”
He joined Hitachi, Ltd. Tokyo, Japan, in 1976.
He was a member of the Visiting Faculty at theLawrence Berkeley Laboratory, University of California, Livermore, CA, from 1978 to 1979. Upon
returning to Japan, he led MRI development projects. He now leads projects onenvironmental measurement and analysis, and noninvasive higher-order brain
function analysis. He is a Chief Research Scientist at the Central ResearchLaboratory, Hitachi, Ltd. and Professor of the Research Institute for Electronic
Science, Hokkaido University, Sapporo, Japan.
Dr. Koizumi received many awards such as the Ohkochi Memorial Prize,
the Prize of Science and Technology Minister, and the IR-100. He is a memberof the Board of Directors of the Chemical Society of Japan.

Yasushi Miyashita received the Ph.D. degree inphysiology from the University of Tokyo, Tokyo,
Japan, in 1979.
He is currently a Professor and Chairman of the
Physiology Department, the University of Tokyo
School of Medicine, Tokyo, Japan. His research interests include neural basis of cognition in primates,
functional imaging, and image processing with MRI,
MEG, and optical devices.