IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 45, NO. 1, JANUARY 1998

87

MEG Covariance Difference Analysis: A Method
to Extract Target Source Activities by Using
Task and Control Measurements
Kensuke Sekihara,* Member, IEEE, David Poeppel, Alec Marantz, Colin Phillips,
Hideaki Koizumi, and Yasushi Miyashita

Abstract— A method is proposed for extracting target dipolesource activities from two sets of evoked magnetoencephalographic (MEG) data, one measured using task stimuli and theother using control stimuli. The difference matrix between thetwo covariance matrices obtained from these two measurements
is calculated, and a procedure similar to the MEG-multiple signal
classification (MUSIC) algorithm is applied to this difference
matrix to extract the target dipole-source configuration. Thisconfiguration corresponds to the source-configuration difference
between the two measurements. Computer simulation verified
the validity of the proposed method. The method was applied
to actual evoked-field data obtained from simulated task-andcontrol experiments. In these measurements, a combination ofauditory and somatosensory stimuli was used as the task stimulus
and the somatosensory stimulus alone was used as the control
stimulus. The proposed covariance difference analysis successfully extracted the target auditory source and eliminated the
disturbance from the somatosensory sources.
Index Terms—Array signal processing, biomagnetics, biomedical electromagnetic imaging, biomedical signal processing, functional brain imaging, inverse problems.

I. INTRODUCTIONNONINVASIVE measurement of human brain functionsis attracting much interest. Such measurements often use
various kinds of sensory stimuli, especially when measuring
higher-order cognitive functions. One general difficulty here
is extracting information on the target brain activity, which isof primary interest in the experiments. This is partly because
spontaneous brain activity overlaps the target activity. Moresubstantially, sensory stimuli generally elicit not only thetarget cortical activities but also other activities that are not

Manuscript received January 16, 1997; revised July 2, 1997. Asterisk
indicates corresponding author.
*K. Sekihara, is with the Mind Articulation Project, Japan Science andTechnology Corporation (JST), Yushima 4-9-2, Bunkyo, Tokyo 113 Japan(e-mail: ksekiha@po.iijnet.or.jp).
D. Poeppel is with the Biomagnetic Imaging Laboratory, University ofCalifornia, San Francisco, CA 94143-0628 USA.
A. Marantz is with the Department of Linguistics, Massachusetts Instituteof Technology, Cambridge, MA 02139 USA.
C. Phillips is with Linguistic and Cognitive Science, University of
Delaware, Newark, DE 19716 USA.
H. Koizumi is with the Mind Articulation Project, Japan Science andTechnology Corporation (JST), Tokyo 113 Japan.
Y. Miyashita is with the Department of Physiology, The University ofTokyo, School of Medicine, Hongo, Tokyo 113 Japan. He is also with theMind Articulation Project, Japan Science and Technology Corporation (JST),Tokyo 113 Japan.
Publisher Item Identifier S 0018-9294(98)00248-1.

the main interest in the experiment. The existence of theseother activities can make interpreting the experimental resultsdifficult.
Therefore, in most studies using positron emission tomography (PET) or functional magnetic resonance imaging(fMRI), neuropsychologists carefully design their experiments
to extract only the target activities. A common example forsuch experimental designs contains two kinds of stimuli: a task
stimulus and a control stimulus. The task stimulus generally
elicits the target cortical activities as well as other activities
associated with them. The control stimulus is designed toelicit only these associated activities. Then, by calculating thedifference between the images measured with each kind of
stimulus, the target activities can be extracted.
Among the various modalities that provide functional brain
images, magnetoencephalography (MEG) [1] has a distinct
advantage—it can provide a temporal resolution of less than
1 ms. MEG measurements have been used to localize various
kinds of cortical activities including the primary sensory areas[2]–[4], spontaneous activities [5], and higher-order corticalactivities [6], [7]. However, experiments using task and control
stimuli are relatively rare in MEG, probably because appropriate data-analysis methods for such measurements do notexist.
In this paper we propose a novel method for extracting thetarget source information from evoked MEG data measured
with task and control stimuli. This method is based on thecovariance difference algorithm originally proposed to eliminate the influence of nonwhite noise with unknown statisticalproperties in sonar-signal processing [8]. Two sets of data aremeasured under different stimulus conditions; the difference
between the covariance matrices obtained from these two
measurements is then calculated. A procedure similar to the
MEG-multiple signal classification (MUSIC) algorithm [9] isapplied to this difference matrix to extract the target source
configuration, which corresponds to the source-configurationdifference between the two measurements.
After describing the proposed method, we present resultsfrom computer simulation and from application to evokedresponse measurements simulating task-and-control-type experiments. Both sets of results strongly suggest that the proposed algorithm is effective. In this paper, we use plain italics
to indicate scalars, lower-case boldface italics to indicate vectors, and upper-case boldface italics to indicate matrices. The

0018–9294/98$10.00  1998 IEEE

88

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 45, NO. 1, JANUARY 1998

superscript T indicates the matrix transpose. The eigenvalues
are numbered in decreasing order. For simplicity, we omit
the explicit time notation unless any confusion due to thisomission arises.
II. METHOD
A. Definitions
We define the magnetic field measured by the mth
detector coil at time tk as bm (tk ); and vector b(tk ) =(b1 (tk ); b2 (tk ); 1 1 1 ; bM (tk ))T as a set of measured data at
tk where k = 1; 2; 1 1 1 ; K: Here, K is the total numberof time points, and M is the total number of detector
coils. The spatio-temporal data matrix B is defined as B =[b (t1 ; b (t2 ); 1 1 1 ; b (tK )]: We assume that a total of P currentdipole sources generate a biomagnetic field. The sphericalhomogeneous conductor model [10] is assumed, and twotangential components, the  and  components, of the source
moment are considered. The magnitude of the pth dipolesource moment is defined as Sp (tk ): Its orientation is definedas the dipole’s normal vector  p (tk ) = (p (tk ); p (tk )), wherek (tk )k = 1: The source magnitude vector at tk is defined ass (tk ) = (S1 (tk ); S2 (tk ); 1 1 1 ; SP (tk ))T : The source temporal
behavior S is defined as S = [s(t1); s(t2 ); 1 1 1 ; s(tK )]:The lead field vectors for the  and  components of theT
pth source are defined as l p = (lp;
1 ; lp;2; 1 1 1 ; lp;M ) and


T


l p = (lp;1 ; lp;2; 1 1 1 ; lp;M ) : We define the lead field vector
for the pth source at tk as l p (tk ): This l p (tk ) is obtained using

B. MUSIC Algorithm
The MUSIC algorithm approach has recently been introduced [9], [11], [12] to avoid this highly multidimensional
search. A distinct advantage of this algorithm is that regardlessof the number of dipole sources, it gives a suboptimal estimateof the source locations by using only a three-dimensional
search in the solution space. We define the measured-datacovariance matrix as R, the covariance matrix of the dipolesource activities as Q: Using (3), we get

R  B B T = L(S S T )LT + N N T  L (Q)LT + 2I(5)where it is assumed that the noise in the measured data iswhite Gaussian noise with variance 2 , and that the noise
is uncorrelated with the signal. Unless some of the source
activities are perfectly correlated with each other, the rankof R is equal to the number of sources P: Therefore, R
has P eigenvalues greater than 2 and M 0 P eigenvalues
equal to 2 : We denote the eigenvectors of R as fe j g,
where j = 1; 2; 11 1 ; M; and define matrices E S and E N as
E S = [e1 ; 1 1 1 ; eP ] and E N = [eP +1 ; 1 1 1 ; eM ]: The span of
the columns in E S is called the signal subspace and that in
E N called the noise subspace.
To estimate the locations of the dipole sources(x 1 ; x2 ; 1 1 1 ; xP ); the MUSIC algorithm takes advantageof the fact that the lead field vector at each xp is orthogonal
to the noise subspace. Namely

L T ej = 0 for j = P + 1; 1 1 1 ; M:

(6)

The lead field matrix for the entire set of P dipole sources
is defined as

Thus, the source locations can be obtained by checking theorthogonality between the lead field vector and the noise
subspace projector, E N E TN : The measure that evaluates thisorthogonality is called the MUSIC localizer, and is proposed
to be [9]L (tk ) = [l1 (tk ); l 2 (tk ); 1 1 1 ; l P (tk )]:J (x) = 1=min (L (x)T E N E TN L(x); L(x)T L(x))

l p (tk ) = p (tk )l p + p (tk )l p :

(1)

(2)

For simplicity, we assume in this paper that all dipole
sources have fixed orientations during measurement. Because
L (tk ) becomes time-independent under this assumption, therelationship between B and S can be expressed, using thesimpler notation L

B = LS + N :(3)

N is the noise matrix defined by N = [n(t1 ); n(t2 ); 1 1 1 ;n(tK )], where n(tk ) is the additive noise at time tk : Theconventional way of estimating the locations of the dipole
sources is based on minimizing the least squares cost function

F = kB 0 L^ S^ k2 = k I 0 L^ L^ T L^ 01L^ T B )k2:
0

0

1

1

(4)

Here, I is the unit matrix, and the estimates of L and S are^ and S^ , respectively. This minimization requires
denoted as L
a 3P -dimensional search, where P is again the number ofsources. Generally, for such a highly multidimensional optimization search, there is no guarantee of obtaining a correct
solution unless we can set the initial estimate very close to
the true solution.

(7)

where min (1; 1) indicates the generalized minimum eigenvalueof the matrix pair given in parenthesis. In this equation, L(x)is expressed as L(x) = [l  (x); l (x)] where l  (x) and l  (x)are the lead field vectors for the  and  components of asource at x:
The MUSIC localizer is calculated in a volume where
sources can exist, and each location where the localizer reaches
a peak is chosen as the location of one dipole source. Note thatthe localizer shown in (7) is derived under the assumption thatthe source orientations are fixed during measurement. It has
been proven, however, that this localizer is also effective for
dipole sources whose orientations vary during measurement[9]. It is also worth noting that the eigenvector corresponding
to the minimum eigenvalue in (7) gives the estimate of thedipole-source orientation.
C. Covariance Difference Analysis
The essence of covariance difference analysis is calculating
the difference between covariance matrices obtained from two
measurements with different stimulus conditions. This analysis
provides the difference in source distributions between the two

SEKIHARA et al.: MEG COVARIANCE DIFFERENCE ANALYSIS

89

measurements, and thus is an appropriate tool for analyzing
the data taken in MEG measurements with a pair of task-andcontrol stimulus conditions.
We assume that T sources are elicited by the task stimulus
and that, among these T sources, C sources are also elicitedby the control stimulus. In this paper, the sources elicited
only by the task stimulus and not by the control stimulus
are called the target sources and the sources elicited by thecontrol stimulus are called the control sources. We define thenumber of target sources as S , such that S = T 0 C
The covariance matrix for the task measurement, T , can be
expressed as

P

P

P

P

P

P
R

holds true. This relationship corresponds to the orthogonality
relationship shown in (6) for the standard MUSIC algorithm.
Thus, corresponding to the localizer in the standard MUSIC
algorithm shown in (7), the locations of the target sources canbe obtained by scanning the localizer calculated from

J (x) = 1=min(L(x)T E Z E TZ L(x);LL(x)T L(x))
where E Z = [ePS +1 ; 1 1 1 ;eeM ]:
0

(15)

0

the above discussion, we made the assumption that all
P : theIncontrol
sources are activated by the task stimulus. This
may not always be true. That is, some control sources may
not appear in the task stimulus condition. The covariancedifference analysis, however, is still effective in such cases
if the condition described below is satisfied.
Let us assume that, among the C control sources, D
sources are only activated in the control stimulus condition
and ~D sources are activated both in the task and the control
conditions. A number from S +1 to S + D is given to eachsource in the first group of control sources. That is, j where= S + 1 111
S + D represent the moment magnitudesof the sources in this group. A number from S + D + 1
to S + C is given to each source in the second group,and their moment magnitudes are represented by j , where= S + D + 1 111
S + C Note that the relationship
C = D + ~D holds.
The covariance matrix for the control sources activated onlyby the control stimulus, D , is thus given by

RT = [l 1;ll2; 1 1 1 ;llPT ]QT [l1;ll2; 1 1 1 ;llPT ]T + W T (8)PP
where QT and W T , respectively, are the source and noise
covariance matrices for the task condition. The covariancePmatrix for the control measurement, RC , is expressed as
P
P P
RC = [l1;ll2; 1 1 1 ;llPC ]QC [l1;ll2; 1 1 1 ;llPC ]T + W C (9)
S
j
P
;
;
P
P
where QC and W C are the source and noise covarianceP P
matrices for the control condition.
PP
To make the covariance difference algorithm valid, weSneed two assumptions: first, the correlation between the target j
P
P
;
;
P
P
:
sources and the control sources are negligibly small, and sec- P
P P
ond, the noise statistical property remains unchanged between
the two measurements. Under the first assumption, we have
Q

PC 
2
Q
O
2
1 1 1 hSPS +1 SPS +PD i 3
hSPS +1 i
QT = OPSS QPS :
(10)
7
..
..
C
PC
..
QD = 64
5: (16).
.
.

2Here, O  is a  2  matrix with all of its elements equal tohSPS +PD SPS +1 i 1 1 1 hSPS +PD izero. The covariance matrix for the target sources, QS , is a
The covariance matrix for the control sources activated by
PS 2 PS matrix explicitly given by~ , is given byboth the task and control stimulus, QD
2 hS 2 i
1 1 1 hS1 SPS i 3
2
3
1
2
i
1
1
1
h
S
h
S
..
..
PS +PD +1 SPS +PC i
+
P
+1
P
5
S
D
:
QS = 4 ...
(11)
.
.
7
..
..
..
Q~ D = 64
5:.
.
.
hSPS S1 i 1 1 1 hSP2 S i
2hSPS +PC SPS +PD +1 i 1 1 1
hSPS +PC i
The covariance matrix for the control sources, QC , is a(17)
PC 2 PC matrix given by~ ; the source-covariance matrix for the task stimulus
2
Using QDhSP2 S +1 i 1 1 1 hSPS +1 SPT i 3
condition
can be expressed as
6
7
.
.
..
..
QC = 4 ..
(12)
5:
"
#
.
P~D
Q
O
S
P
hSPT SPS +1 i 1 1 1 hSP2 T i
S
QT = OPS Q~ :
(18)
D
P~D
The second assumption leads to W T = W C : Therefore
1R, which is the difference between RT and RC , can be At this point, we make further assumption that the correlation
expressed as

1R = RT 0 RC = [l1;ll2; 1 1 1 ;llPS ]QS [l1;ll2; 1 1 1 ;llPS ]T :
(13)

1R

P

has S nonzero eigenvalues
It was shown in [8] thatand0 S zero-level eigenvalues, and that eigenvectors
corresponding to zero-level eigenvalues are orthogonal to thelead field vector at the target source locations. That is, defining
such eigenvectors as 0j , the relationshipT 0= S +1 111[ 1 2 1 1 1 PS ] j = 0 for
(14)

M P

l ;ll ; ;ll

e

e

j P

; ;M

between the two group of the control sources is negligibly
small. Then, the source-covariance matrix for the control
stimulus condition is given by

"

#

P~D
QC = OQPDD OQ~PD :D
P~DThus, in this case, the difference matrixpressed as(19)

1R can be ex-

1R = RT 0 RC = [l1; 1 1 1 ;llPS ;llPS +1; 1 1 1 ;llPS +PD ](20)
1 1Q[l 1 ; 1 1 1 ;llPS ;llPS +1 ; 1 1 1 ;llPS +PD ]T90IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 45, NO. 1, JANUARY 1998

where


PD 
1Q = OQPSS 0OQPS :
D
PD

(21)

1R

It can be shown that the zero-level eigenvectors ofareorthogonal to [ 1 1 1 1 PS +PD ] [8]. That is, a localizer searchusing (15) not only detects the target sources but also thecontrol sources that are activated only by the control stimulus.
These two groups of sources can, however, be discriminatedby checking the signs of the diagonal elements in matrix,
which can be estimated usingl ; ;ll

1Q

i

h

h

1dQ = 0LT L101LT 1R 0LT L101LT

iT

(22)

where

L

L = [l1; 1 1 1 ;llPS +PD ]:To calculate , the orientations of the dipole sources need
to be estimated using the eigenvector corresponding to thesmallest eigenvalue in (15). An example of applying (22) todiscriminate the two groups of sources will be presented in
Section III-D.
D. Summary of Algorithm
The covariance difference algorithm is summarized as follows.
Step 1) The covariance matrices of the measured spatiotemporal data are calculated for the task andcontrol measurements. Then, the difference matrix between these two covariance matrices iscalculated.
Step 2) Eigen-decomposition of this difference matrix is
performed, and the subspace spanned by zerolevel eigenvectors is defined by separating thezero-level eigenvalues from the nonzero eigenvalues.
Step 3) The localizer defined by (15) is calculated
throughout the field of view, and each location
where ( ) reaches a peak is determined to be
the location of one dipole source.
Step 4) When it is possible that the results obtained in
Step 3 contain control source locations, the signsof the diagonal elements ofare checked using(22) to separate the target sources from control
sources.

Jx

1Q

Fig. 1. Source and detector configuration assumed in computer simulation.

are expressed in centimeters. Three signal dipole sources wereassumed to exist on a plane defined as= 1.0: the first
dipole source was at (01.1, 1.0, 05.6), the second was at(2.8, 1.0, 05.4), and the third was at (4.5, 1.0, 07.3). The
source and detector configuration for this simulation is shown
schematically in Fig. 1.
The spherical homogeneous conductor model [10] withthe origin set at (1.0, 1.0, 010) was used in this computer
simulation. The simulated magnetic field was calculated at 1ms intervals from 0 to 400 ms. To generate the simulated
biomagnetic field, thecomponents of these sources weremodeled using exponentially damped sinusoidal functions, andthe components were set to zero. The waveforms assigned
to each component of the three sources are shown in Fig. 2.
Here, the correlation coefficients between any pair of the three
waveforms were less than 0.3.
Uncorrelated Gaussian noise was added to make the final
signal-to-noise ratio (SNR) equal to 0.3 for the single-epoch
data in the task condition. The SNR was defined by the ratioof the Frobenius norm of the signal-magnetic-field data matrix
to that of the noise matrix. One hundred sets of the magneticfield data were generated and averaged to create the magnetic
field data used for the source estimation. The final SNR ofthis averaged field data was approximately three for the task
condition.

y







B. Comparison with the Waveform-Subtraction Method
III. COMPUTER SIMULATION
A. General Data-Generation Conditions
We used computer simulation to test the validity of the
proposed method. A 37-channel magnetometer whose coilalignment was the same as that of the Magnes biomagnetic
measurement system (Biomagnetic Technologies Inc., San
Diego) was assumed. The
direction was defined as the
direction perpendicular to the plane of the detector coil locatedat the center of the coil alignment, and was equal to zero at
this coil plane. The values of the spatial coordinates ()

zzx;y; zWe first assumed that the task condition activates all threesources and that the control condition activates only the first
and third sources. In this case, the second source is thetarget source. In actual task-and-control-type measurements,
the onset of the control-source response may differ between
the task and control measurements. To take this possibility into
account, the onset of control-source activation was set to bedelayed by 30 ms in the control measurement. The activation
waveforms shown by the solid lines in Fig. 2 are used in thetask condition and those shown by the broken lines are usedin the control condition. The averaged magnetic-fields for theSEKIHARA et al.: MEG COVARIANCE DIFFERENCE ANALYSIS

91

(a)

(b)

(c)
Fig. 2. The waveforms assigned to the  components of (a) the first, (b) the second, and (c) the third source in the computer simulation. In Section III-B,
the waveforms indicated by the solid lines are used in the task condition and those by the broken lines in the control condition.

task and control conditions are shown in Fig. 3(a) and (b).
One of the single-epoch data is shown in Fig. 3(c).
The standard MUSIC localizer in (7) was applied to theaveraged field data; the results of calculating the localizer on
the y = 1:0 plane are shown in Fig. 4(a) and (b). Fig. 4(a)shows the results for the task condition and (b) shows thosefor the control condition. The contours in these figures show
the relative value of the localizer, and each area where thelocalizer reaches a peak is considered to be the location of
one dipole source. In Fig. 4(a) the localizer clearly detects
the three signal sources, and in Fig. 4(b) it clearly detects
the first and third sources. The results of applying covariancedifference analysis are shown in Fig. 4(c). In this case, only
the target source (the second source) remains—the influenceof the control sources (the first and third sources) has been
completely eliminated.
A naive way to extract information about the target activity
from the task and control measurements is to directly subtract
the field waveform elicited under the control condition fromthat under the task condition. The results of this simple waveform subtraction are shown in Fig. 4(d). Here, the control field
waveform was subtracted from the task waveform, and thestandard MUSIC localizer [(7)] was applied to this subtracted
waveform. These results show that the simple subtraction

method cannot separate the target source activity from thecontrol activities due to the 30-ms time offset existing in
the control source activations between the two conditions. A
comparison between Fig. 4(c) and (d) clearly demonstrates theeffectiveness of the proposed covariance difference analysis.
C. Effects of Strong Correlation Between
Task and Control Sources
We conducted the computer simulation when the target
source had a strong correlation with one of the control sources.
The waveform indicated in Fig. 2(a) was assigned to both thefirst and second sources. The same data generation procedure
as in Section III-B was repeated, and the covariance difference
analysis was applied to the generated field data. The results areshown in Fig. 5. Here, the elimination of the control source
influence obviously failed. These results demonstrate that the
proposed covariance difference algorithm is not effective whenstrong correlation exists between the task and control sources.
D. Influence of a Control Source not
Elicited in the Task Condition
We next assumed that the task condition activates the first
and third sources and that the control condition activates the92IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 45, NO. 1, JANUARY 1998

(a)

(b)

(c)
Fig. 3. The 37-channel overlapped display of the field waveform generated in the computer simulation in Section III-B: (a) the task condition, (b) the control
condition, and (c) one of the 100 single-epoch data generated in the task condition. This data is used to obtain the results shown in Fig. 7(b).

(a)

(b)

(c)

(d)

Fig. 4. Results of the MUSIC source localization when the task condition activates all three sources and the control condition activates the first and thirdsources: (a) Results of applying the standard MUSIC localizer [(7)] to the data obtained under the task condition, (b) results of applying (7) to the dataobtained under the control condition, (c) results of the covariance difference analysis using (15), and (d) results of the simple waveform-subtraction method.
The contours show the relative values of the localizer on the y = 1:0 plane, and each area where the localizer reaches a peak is considered to be the locationof one dipole source. The circle is the boundary of the sphere used for the forward calculation, depicting the approximate size of a human brain on this plane.

second and third sources. In this case, the first source is thetarget source and the second source is activated only in thecontrol condition. The source waveforms shown by the solidlines in Fig. 2(a), (b), and (c) were assigned to the first, second,
and third sources both for the task and control conditions. The
results obtained by applying the standard MUSIC localizer
to the task and control data are shown in Fig. 6(a) and (b),respectively. According to the discussion in Section II-C, thecovariance difference analysis should detect the control (thesecond) source as well as the target (the first) source. The
results of the covariance difference analysis are shown in
Fig. 6(c). In these results, the first and second sources areby using (22);detected. We then estimated the matrix
the diagonal elements were found to be 33.8 (the first source)

1Q

SEKIHARA et al.: MEG COVARIANCE DIFFERENCE ANALYSIS

Fig. 5. The results of the covariance difference analysis when strong correlation exists between the task and control sources. The waveform shown in
Fig. 2(a) was assigned to both the first and second sources in this computer
simulation.

and 042.71 (the second source). The signs of these diagonal
elements indicate that the first source is the target source
and the second source is the control source. These results
demonstrate the validity of the discussion in Section II-C.
E. Influence of Time Jitter in Control Source Response
It is possible that the onset of the control source response
will have a time jitter and vary from epoch to epoch inactual cognitive experiments. In such cases, the proposed
covariance difference analysis cannot effectively eliminate thecontrol source influence if it is applied to the data obtainedby averaging multiple epochs. The same computer simulation
as described in Section III-B was again performed exceptthat the onset of the control source in the task measurement
was randomly changed within a maximum of 30 ms whengenerating data for each epoch. The covariance difference
analysis was applied to the data averaged over 100 sets ofthese single-epoch data.
The results are shown in Fig. 7(a). Obviously, the control
source influence was not eliminated. The results obtained when
the covariance difference analysis was applied to a pair of thesingle-epoch data are shown in Fig. 7(b). One of the single
epoch data used here is shown in Fig. 3(c). The results containa severe blur, due to the low signal-to-noise ratio of the single
epoch data, even though the control source influence wasremoved. This blur, however, can be reduced by averaging
the difference matrix obtained using each set of single-epochdata. The results of using this averaged-difference matrix are
shown in Fig. 7(c). To obtain these results,100 single-epoch
difference matrices were averaged, and the proposed localizer[(15)] was calculated using this averaged difference matrix.
The quality of these results is almost equal to that of the resultsin Fig. 4(c). Therefore, when one suspects that the onset of thecontrol-source response may vary, the procedure mentioned
in this subsection should be applied to eliminate the control
source influence.
IV. APPLICATION TO EVOKED-FIELD MEASUREMENTS
SIMULATING TASK- AND CONTROL-TYPE EXPERIMENTS
We applied the proposed method to source localization formeasured evoked responses; these measurements simulated
task- and control-type experiments. The evoked responses
1 These elements are expressed using an arbitrary scale.

93

were recorded using the 37-channel Magnes magnetometer
installed at the Biomagnetic Imaging Laboratory, Universityof California, San Francisco. All measurements were done ina magnetically shielded room.
In these experiments we applied an auditory stimulus, asomatosensory stimulus, or both to a male volunteer. The
auditory stimulus was a 1000-Hz pure tone with a 200ms duration; it was applied to the subject’s right ear. The
somatosensory stimulus was a 30-ms-duration tactile pulse (17psi) delivered to the distal segment of the right index finger.
In the combination stimulus, the auditory and somatosensory
stimuli started at the same time. These auditory, somatosensory, and combination stimuli were repeatedly given one after
the other without changing the head position relative to thesensor array. The data was acquired at a sampling frequencyof 1 kHz for the prestimulus interval of 300 ms and thepoststimulus interval of 800 ms and averaged for 256 epochsof each stimulus condition. An on-line bandpass filter with abandwidth from 1–400 Hz was used and no post-processing
digital filter was applied. The sensor array was placed on
the left hemisphere and positioned to best record the M100
auditory response. The mean interstimulus interval was 2 s,
randomly varied between 1.75 and 2.25 s.
The x; y; and z coordinates used to express the localization
results are depicted in Fig. 8. The results of applying thestandard MUSIC algorithm are shown in Figs. 9, 10, and 11.
The data from 0–300 ms post stimulus onset were used forthe analysis. In these figures, the localizer shown in (7) was
calculated with an interval of 0.5 cm within a volume definedas 04  x  6, 03  y  6, and 3  z  11: Thesefigures show the projections of the localizer values onto thetransverse, coronal, and sagittal planes; the relative positionsof these planes are shown in Fig. 8. The circles depicting ahuman head represent the projections of the sphere used to
calculate the forward solutions.
The results obtained using the data measured with only
the auditory stimulus are shown in Fig. 9. The sharp single
peak corresponds to a source in the auditory cortex. The
results from the somatosensory stimulus alone are shown in
Fig. 10. Two sources, probably corresponding to sources atthe primary (SI) and secondary somatosensory (SII) cortices,can be seen, although they are not well resolved. The resultsfrom the combination of auditory and somatosensory stimuli
are shown in Fig. 11. Two peaks can be seen in the left
hemisphere. The sharp peak probably corresponds to thesource at SI; the dull peak probably corresponds to the twosources at the auditory and secondary somatosensory cortices;
these two sources cannot be resolved because they are very
close together.
We calculated two covariance matrices for the covariancedifference analysis: one by using the data from the combined
stimulus and the other by the data from the somatosensory
stimulus. That is, the combined stimulus was assumed to be thetask stimulus, and the somatosensory stimulus was assumed
to be the control stimulus. The data portion from 0–300 ms
were also used for calculating both of the two covariancematrices. The results of the covariance difference analysis areshown in Fig. 12. A single sharp peak exists at almost the94IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 45, NO. 1, JANUARY 1998

(a)

(b)

(c)

Fig. 6. Results of the MUSIC source localization when the task condition activates the first and third sources and the control condition activates thesecond and third sources: (a) Results for the data obtained from task condition, (b) results for the data from control condition, and (c) results of thecovariance difference analysis using (15).

(a)

(b)

(c)

Fig. 7. Results of the same computer simulation as in Fig. 4 for the case when the onset of the control source response was randomly varied within
30 ms when generating each set of epoch data: (a) Results of the covariance difference analysis applied to the data obtained by averaging 100 sets ofsingle-epoch data, (b) results of the covariance difference analysis applied to a pair of single epoch data, and (c) results when 100 single-epoch difference
matrices were averaged, and the localizer in (15) was calculated using the averaged matrix.

Fig. 8. The x; y; and z coordinates used to express the localization results
shown in Figs. 9–12. The midpoint between the left and right preauricularpoints is defined as the coordinate origin. The axis directed away from theorigin toward the left preauricular point is defined as the y axis, and thatfrom the origin to the nasion is the x axis. The z axis is defined as theaxis that is perpendicular to both these axes and directed from the origin to
the vertex. The relative positions of the transverse, sagittal and coronal planes
are also shown.

+

+

+

same location where the source from the auditory stimulus
was localized (Fig. 9). These results clearly show that thecovariance difference analysis can extract the target (auditory)
source location and eliminate the disturbance from the control(somatosensory) sources.
V. CONCLUSION
As mentioned in Section III-B, a naive method to extract the
information about the target activities is to simply subtract themagnetic-field waveform obtained by the control measurementfrom that by the task measurement. The advantages of the

proposed covariance difference algorithm over this simple
subtraction method are summarized as follows.
First, the simple subtraction method is effective only when
the onset of the control-source time response in the task
measurements is exactly equal to that in the control measurements. Otherwise, this method fails to extract the target-source
information, as was shown in Fig. 4(d). Second, the simple
subtraction method also cannot handle the case where some of
the control sources do not appear in the task measurements. In
this case, even if all sources are correctly localized, the simple
subtraction method provides no procedure for discriminating
between the target sources and the control sources. In contrast,
covariance difference analysis can discriminate between these
two types of sources by checking the signs of the diagonal
obtained using (22), as discussed inelements of matrix
Section II-C.
The assumption that makes the proposed covariance difference analysis valid is that the brain response is linear, i.e.
the brain responses to different stimuli are additive and thecorrelation between these responses is zero. It should be notedthat the validity and the limitations of the linear modelingof brain response have been studied in the functional MRI[13], [14], and the linear modeling proves to be effective foranalyzing the fMRI or PET data in ordinary measurement
conditions [15], although the brain response is known to be,in principle, highly nonlinear.
The success of the proposed covariance difference analysis
depends on the assumption that this linear modeling is also
valid for MEG measurements. In actual MEG measurements,d

1Q

SEKIHARA et al.: MEG COVARIANCE DIFFERENCE ANALYSIS

95

Fig. 9. Results of applying the standard MUSIC localizer in (7) to data obtained with the auditory stimulus. The localizer was calculated with an intervalof 0.5 cm within a volume defined as 4 x 6; 3 y6, and 2 z 10; the projections of the localizer values onto the transverse, coronal, andsagittal planes are shown. The circles depicting a human head represent the projections of the sphere used to calculate the forward solution.

0   0   Fig. 10.

Results of applying the standard MUSIC localizer in (7) to data obtained with the somatosensory stimulus.

Fig. 11.

Results of applying the standard MUSIC localizer in (7) to data measured with combination of the auditory and somatosensory stimuli.

Fig. 12. Results of the covariance difference analysis. One covariance matrix was calculated from the auditory-somatosensory combined response and theother was calculated from the somatosensory response. The localizer in (15) was used.

however, this assumption may not always be valid, so the
proposed method would not be effective in some cases.
The validity and limitations of the linear model for theneuromagnetic measurements should be investigated, and theusefulness and limitations of the proposed method must be
evaluated in such investigations.

96

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 45, NO. 1, JANUARY 1998

In summary, we have proposed a method for extracting
target source activities, which are of primary interest in
experiments, by using two sets of evoked MEG data: one
set measured by using task stimulus and the other by using
control stimulus. The proposed method is a variant of the
MEG-MUSIC algorithm, that is, a procedure similar to theMUSIC algorithm is applied to the difference between thetwo covariance matrices obtained from the above two measurements. Computer simulation showed the validity of the
proposed method. Application to evoked-field measurements
strongly suggested the method’s effectiveness.

ACKNOWLEDGMENT
The authors wish to thank S. Masterman for helping toprepare the biomagnetic data. This work was carried out aspart of the MIT–JST International Research Project “Mind
Articulation.”

REFERENCES[1] M. Hämäläinen, 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.
[2] G. Romani, S. Williamson, and L. Kaufman, “Tonotopic organizationof the human auditory cortex,” Sci., vol. 216, pp. 1330–1340, 1982.
[3] R. Hari, M. Hämäläinen, and S. L. Joutsiniemi, “Neuromagnetic steadystate responses to auditory stimuli,” J. Acoust. Soc. Amer., vol. 86, pp.
1033–1039, 1989.
[4] N. Fross, R. Hari, R. Salmelin, A. Ahonen, M. Hämäläinen, M. Kajola,
J. Knuutila, and J. Simola, “Activation of the human posterior parietal
cortex by median nerve stimulation,” Exp. Brain Res., vol. 99, pp.
309–315, 1994.
[5] R. Hari and R. Salmelin, “Human cortical oscillations: A neuromagnetic
view through the skull,” Trends Neurosci., vol. 20, pp. 44–49, 1997.
[6] K. Sasaki, T. Tsujimoto, A. Nambu, R. Matsuzaki, and S. Kyuhou,“Dynamic activities of the frontal association cortex in calculating andthinking,” Electroenceph. Clin. Neurophysiol., vol. 74, pp. 58–75, 1989.
[7] M. J. Liu, P. B. C. Fenwick, J. Lumsden, C. Lever, K.-M. Stephan, and
A. A. Ioannides, “Averaged and single-trial analysis of cortical activation
sequence in movement preparation, initiation, and inhibition,” Human
Brain Mapping, vol. 4, pp. 254–264, 1996.
[8] A. Paulraj and T. Kailath, “Eigenstructure methods for direction ofarrival estimation in the presence of unknown noise fields,” IEEE Trans.
Acoust., Speech, Signal Processing, vol. ASSP-34, pp. 13–20, 1986.
[9] 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.
[10] J. Sarvas, “Basic mathematical and electromagnetic concepts of thebiomagnetic inverse problem,” Phys. Med. Biol., vol. 32, pp. 11–22,
1987.
[11] R. O. Schmidt, “A signal subspace approach to multiple emitter locationand spectral estimation,” Ph.D. dissertation, Stanford Univ., Stanford,
CA, 1981.
[12] A. Paulraj, B. Ottersten, R. Roy, A. Swindlehusrt, G. Xu, and T.
Kailath, “Subspace methods for directions-of-arrival estimation,” in
Handbook of Statistics, N. K. Bose and C. R. Rao, Eds. Amsterdam,
the Netherlands: Elsevier Science, 1993, pp. 693–739.
[13] G. M. Boynton, S. A. Engel, G. H. Glover, and D. J. Heeger, “Linear
analysis of functional magnetic resonance imaging in human V1,” J.
Neurosci., vol. 16, pp. 4207–4221, 1996.
[14] K. J. Friston, C. J. Price, P. Fletcher, C. Moore, R. S. J. Frackowiak,
and R. J. Dolan, “The trouble with cognitive subtraction,” Neuroimag.,vol. 4, pp. 97–104, 1996.
[15] K. J. Friston, A. P. Holmes, K. J. Worsley, J.-P. Polin, C. D. Frith, and R.
S. J. Frackowiak, “Statistical parametric maps in functional imaging: A
general linear approach,” Human Brain Mapping, vol. 2, pp. 189–210,
1995.

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

David Poeppel obtained the B.S. and Ph.D. degrees
in cognitive neuroscience at Massachusetts Instituteof Technology (MIT), Cambridge, in 1990 and
1995, respectively.
He is currently an Adjunct Assistant Professor
in the Department of Radiology, University of California San Francisco, and an 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, in 1981.
He joined the faculty at MIT in 1990, where he iscurrently Professor of Linguistics in the Departmentof Linguistics and Philosophy. His research interests
include the syntax and morphology of natural languages, linguistic universals, and the neurobiologyof language. He is currently involved in revisingmorphological theory within linguistics and in exploring MEG techniques to
uncover how the brain process language.

Colin Phillips received the B.A. degree in modern
languages from Oxford University, Oxford, U.K.,
in 1990 and the Ph.D. degree in linguistics from
Massachusetts Institute of Technology (MIT), Cambridge, in 1996.
He is currently an Assistant Professor in the
Department of Linguistics and the Program in Cognitive Science at the University of Delaware. His
research combines theoretical linguistics with language processing, language acquisition and neurolinguistics, with the goal of embedding explicit
linguistic theories in cognitively and neuroscientifically grounded models of
language processes.

SEKIHARA et al.: MEG COVARIANCE DIFFERENCE ANALYSIS

Hideaki Koizumi was born in Tokyo, Japan, in1946. He was graduated from the Department ofPure and Applied Sciences, the University of Tokyo,
Tokyo, Japan, in 1971. He joined Hitachi, Ltd.
Tokyo, Japan, and received the doctoral degree inphysics from the University of Tokyo in 1976 for histhesis: “Creating Polarized Zeeman-Effect Atomic
Absorption Spectrometry.”
He stayed as Visiting Faculty at the Lawrence
Berkeley Laboratory, University of California, Livermore, CA, from 1978–1979. Upon returning toJapan, he led MRI development projects. He is a Chief Research Scientist at
the Central Research Laboratory, Hitachi, Ltd. and a Professor of the ResearchInstitute for Electronic Science, Hokkaido University, Sapporo, Japan. He now
leads projects on environmental measurement and analysis, and noninvasive
higher-order brain function analysis.
Dr. Koizumi is a Member of the Board of Directors of the Chemical Societyof Japan. He received many awards such as the Ohkochi Memorial Prize, thePrize of Science and Technology Minister, and IR-100.

97

Yasushi Miyashita received the Ph.D. degree inphysiology from the University of Tokyo, Japan, in
1979.
He is currently 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.