1534

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 49, NO. 12, DECEMBER 2002

Performance of an MEG Adaptive-Beamformer
Technique in the Presence of Correlated Neural
Activities: Effects on Signal Intensity andTime-Course Estimates
Kensuke Sekihara3 , Srikantan S. Nagarajan, David Poeppel, and Alec Marantz

Abstract—The influence of temporarily correlated source
activities on neuromagnetic reconstruction by adaptive beamformer techniques was investigated. It is known that the spatial
filter weight of an adaptive beamformer cannot perfectly block
correlated signals. This causes two major influences on thereconstruction results: time course distortions and reductions
in reconstructed signal intensities. Our theoretical analysis andnumerical experiments both showed that the reduction in signal
intensity for sources with a medium degree of correlation is small.
The time-course distortion for such sources, however, may be discernible. Our analysis also showed that the magnitude correlation
coefficient between two correlated sources can be accurately estimated by using the beamformer outputs. A method of retrieving
the original time courses using estimated correlation coefficientswas developed. Our numerical experiments demonstrated thatreasonably accurate time courses can be retrieved from considerably distorted time courses even when the signal-to-noise ratio islow.
Index Terms—Adaptive beamformer, biomagnetism, functional
neuroimaging, magnetoencephalographic (MEG) inverse problems, magnetoencephalography, neuromagnetic signal processing.

I. INTRODUCTIONT

HE goal of neuromagnetic imaging is to visualize neuralactivities at a fine time resolution in the order of milliseconds, and to provide functional information about brain
dynamics [1], [2]. Toward this goal, a number of algorithms
for reconstructing spatio-temporal source activities frommagnetoencephalographic (MEG) measurements have been
Manuscript received January 8, 2002; revised July 15, 2002. This work was
supported by Grants-in-Aid from the Kayamori Foundation of Informational
Science Advancement, Grants-in-Aid from the Suzuki Foundation, andGrants-in-Aid from the Ministry of Education, Science, Culture and Sports inJapan (C13680948). This work was also supported in part by the Whitaker
Foundation and by National Institute of Health (R01-DC004855-01). Asterisk
indicates corresponding author.
K. Sekihara is with the Department of Electronic Systems and Engineering,Tokyo Metropolitan Institute of Technology, Tokyo 191-0065, Japan (e-mail:ksekiha@cc.tmit.ac.jp).
S. S. Nagarajan was with the Department of Bioengineering, University of
Utah, Salt Lake City, UT 84112-9202 USA. He is now with the Department ofRadiology, the University of California, San Francisco, CA 94143 USA.
D. Poeppel is with the Department of Linguistics and Biology, University of
Maryland, College Park, MD 20742 USA.
A. Marantz is with the Department of Linguistics and Philosophy, Massachusetts Institute of Technology, Cambridge, MA 02139 USA.
Digital Object Identifier 10.1109/TBME.2002.805485

developed. Among such algorithms, a class of techniques called
an adaptive beamformer have attracted great interest recently.
Adaptive-beamformer techniques were originally developedin the fields of array signal processing, including radar, sonar,and seismic exploration [3]. Since first introduced by Robinson
et al. [4], they have been successfully applied to neuromagnetic source reconstruction problems [5]–[9]. The details ofadaptive beamformer techniques, as well as their comparison
to least-squares-based source-reconstruction techniques suchas the minimum-norm method [10], are given in [11].
Adaptive-beamformer techniques, however, require an implicit assumption that the time courses of the source activities
are orthogonal to each other; i.e., all source activities are completely uncorrelated [12]. This assumption may not necessarily
hold in regard to neuromagnetic measurements. In general, cortical and subcortical activities are supposed to be more or lesspartially correlated. Indeed, this correlation itself has become an
active research target [13]–[15] because the correlation between
the neural activities is considered to represent some sort of functional connectivity between different regions in the brain [16].
Therefore, when applying adaptive-beamformer techniques, weshould be aware of how robust they are to the source correlation,
and what kind of influences arise when they are used to reconstruct correlated sources.
This paper attempts to answer these questions. It is known
that the weight vector of adaptive beamformers cannot perfectly
block correlated signals. This causes two major influences on
the source-reconstruction results: distortion in time-course estimates and reduction in reconstructed source intensities. We
performed theoretical analysis under the condition that two correlated sources exist and derived theoretical relationships between the degree of these influences and the degree of source
correlations. We also develop a method of retrieving the originaltime courses when two sources are correlated. In this paper, Section II briefly reviews existing adaptive-beamformer techniques
for neuromagnetic reconstruction. Section III presents our theoretical analysis regarding the influence of the source correlationon the adaptive beamformer performance. Section IV presents
the results of several numerical experiments that validate the arguments in Section III. Throughout this paper, plain italics indicate scalars, lower case boldface italics indicate vectors, andupper case boldface italics indicate matrices.

0018-9294/02$17.00 © 2002 IEEE

SEKIHARA et al.: PERFORMANCE OF AN MEG ADAPTIVE-BEAMFORMER TECHNIQUE

II. ADAPTIVE-BEAMFORMER TECHNIQUES FOR
NEUROMEGNETIC RECONSTRUCTION
A. Definitionsth
We define the magnetic field measured by theas, and a column vector
detector coil at time
as a set of measured datais the total number of detector coils and superscriptwhere
indicates the matrix transpose. The spatial location is repre.
sented by a three-dimensional (3-D) vector :discrete sources are assumed to generate the
A total ofneuromagnetic field. These discrete sources are modeled using
the equivalent current dipoles [17]. The number of sources
is assumed to be less than the number of sensors ; i.e., themeasured neuromagnetic signal is assumed to be a low-rank
signal [18], [19]. The locations of these sources are denoted as. The moment magnitude of the th source at, and the source magnitude vector
time is denoted as.
is defined as
The orientation of the th source is defined as a 3-D column
whosevectorcomponent (where equals , , or in this paper) is equal
to the cosine of the angle between the direction of the source
moment and the direction. We assume that the orientation ofeach source is time independent. Omitting time notation , wematrix that expresses the orientations of all
define asources as such that
..
.
..
.

..

.

We defineas the output of the th sensor; the output
is induced by the unit-magnitude source located at and diis defined asrected in the direction. The column vector. Next, the lead field matrix,
which represents the sensitivity of the whole sensor array at ,. The lead-field vector
is defined asin the source-moment direction is defined, for later use, as
where. The composite lead field matrix for theentire set of sources is defined as

The relationship between

andThe measurement covariance matrix is denoted as , which is, whereindicates the enobtained fromsemble average, which is replaced with the time average over acertain time window in practice. It should be noted that whenis not equal to zero,is not equal to the covariance matrix and should rather be referred to as the second-order momentismatrix in such a case. However, according to convention,referred to as a covariance matrix in this paper. The covariancematrix of the source-moment activity is defined here as ; i.e.,.
and the source acUsing (3) and assuming that the noise
are uncorrelated, we get the relationship between thetivitymeasurement covariance matrix and the source-activity covariance matrix such that
(5)
where the noise in the measured data is assumed to be the white
Gaussian noise with a variance of and is the identity matrix.
are defined as
The th eigenvalue and the eigenvector ofand , respectively. Unless some source activities are perfectly
is equal to the number
correlated with each other, the rank ofhas eigenvaluesof sources . Therefore, according to (5),eigenvalues equal to . For later use,
greater than andis defined as. The column
the matrix
is the maximum-likelihood estimate of the signal
span of[20].
subspace of
It should be pointed out that, in (5), interference magnetic
fields originated from external noise sources such as brain background activities can be considered as a part of signals, and theirinfluence can be accounted for in the first term of the right-handside of this equation. The noise here indicates only noises withinternal origins, such as those originated from the sensor coils,
the superconducting quamtum interference devices (SQUIDs),and their associated electronics. Therefore, the assumption of
spatially uncorrelated white Gaussian noise is acceptable for deriving (5).
B. Scalar Minimum-Variance Beamformer
To estimate the source moment, we have focused on the classof techniques referred to as the adaptive beamformer [21]. This
type of technique uses the following spatial-filter operation for
estimating the source moment:

(1)

(6)

is then expressed asis the source-moment time course obtained as
where
the beamformer output. One well-known adaptive beamformer
is the minimum-variance distortion-less beamformer originallydeveloped in the field of array signal processing [22]. Thisby minimizing
technique derives the weight vector
under the constraint of. The
explicit form of the weight vector for the minimum-variancebeamformer is known to be(2)is the additive noise. Defining
whererelationship is rewritten as

1535

as

, this
(3)

where
(4)
is the composite lead field matrix whose column
The matrix
is equal to the lead field vector in the source-moment direction.

(7)
isNote that to use this formula, the source orientation. However,is generally unknown,needed to calculate

1536

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 49, NO. 12, DECEMBER 2002

although several techniques have been developed to estimate it[12], [23].

III. BEAMFORMER PERFORMANCE WHEN
SOURCES ARE CORRELATED
A. Signal Leakage and Errors in Time-Course Estimate

C. Vector Minimum-Variance Beamformer Formulation
A vector-type beamformer uses a set of three weight vec,
, and, each of which estimates the ,tors,, and components of the source moment, thus enabling thesource orientation as well as the source-moment magnitude tobe estimated. A set of weight vectors for a vector-extended minimum-variance beamformer is given by [24], [25]The weight vector for the minimum-variance beamformer is obtained by minimizing the output power ;under the constraint.
When the beamformer pointing location is equal to the source
location , neglecting the noise term, the output power isrewritten as(8)Using these weight vectors, the , , andsource moment are given bycomponents of the(9)is the estimated source moment in the direcwheretion. The estimate of the source-moment orientation is denoted, which can be extracted from (9) such thatas

(10)
The magnitude of the source moment is estimated from(14)
where sources are assumed to be uncorrelated with each other,is used. Thus, the weight
and the relationship
vector obtained by minimizing should have the property that
where
. That is, the resultant weightshould only pass the signal from , and block the signal fromother sources, and the blocking capability of the beamformer
weight is expressed as

(11)
(15)
. It can bewhereis equivalent toseen that the composite weight vectorin (7). In thethe weight vector for the scalar beamformer
for simplicity, but the results offollowing analysis, we use.
the analysis are also valid foris the Kronecker delta defined as
(
)
where
(
).
and
In general cases where partially correlated sources exist, however, the blocking capability of the weight vector is expressedin the following equation:

D. Eigenspace-Projected Beamformer
The minimum-variance beamformer is known to be very
sensitive to errors in the forward modeling or errors in estimating the measurement covariance matrix. In the case of
neuromagnetic measurements, since such errors are almost
inevitable, the minimum-variance beamformer generally
provides noisy results. The eigenspace-projection beamformer
provides an output signal-to-noise ratio (SNR) much higher
than that of the minimum-variance beamformer when such
errors mentioned above exist [26]. The extension of theminimum-variance beamformer to the eigenspace-projection
beamformer is attained by projecting the weight vector of theminimum-variance beamformer onto the signal subspace of themeasurement covariance matrix. That is, redefining the weight, the eigenspace-projection
vector obtained from (7) as
beamformer is given by [26](12)or redefining the weight vectors obtained from (8) as, the vector-extended eigenspace
beamformer is given by [7]

(13)

(16)
indicates theelement of. This equawheretion, first derived by Zoltowski [27], is the basis of our analysis.
Its derivation is described in Appendix I and, as shown in thisderivation, (16) exactly holds when the SNR is so high that noise
can be neglected. We assume that the target source exists atand that thesources are correlated with the target source.
sources are mutually correlated). The beam(A total of, is expressed as
former output at ,(17)
where the locations of the correlated interferences are denoted). This equation shows that the beamas , (former outputs for the correlated sources contain leakage fromthe other correlated sources, and such leakage causes errors in
the time-course estimates of the source activities.
We consider a case where, among sources, the first and thesecond sources are significantly correlated and the other sources
have no significant correlation with any other sources. We define the correlation coefficient between the first and the secondsources as , and the average power of the th source as ,

SEKIHARA et al.: PERFORMANCE OF AN MEG ADAPTIVE-BEAMFORMER TECHNIQUE

1537

TABLE I
SOURCE PARAMETER VALUES USED FOR THE
NUMERICAL EXPERIMENTS IN SECTION IV

where. The source covariance matrix and its
inverse can be expressed as in (18) and (19), shown at the bottomof the page.
Using (16) and (19), we can derive(20)

Fig. 1. The coordinate system and source-detector configuration used in
the numerical experiments. The coordinate origin was set at the center of thedetector coil located at the center of the coil array. The plane at x = 1:0 cm
is shown. The large circle shows the cross section of the sphere used for theforward calculation.

and
for

TABLE IIVALUES OF THE PARAMETERS USED FOR
CALCULATING w (t), w (t), AND w (t)(21)

Therefore, when the beamformer is pointing at the first source
location , it passes the signal from the second source with, although it blocks thethe multiplicative constant ofsignals from the other sources uncorrelated with the first source.
, is given byConsequently, the beamformer output at ,(22)
In exactly the same manner, the beamformer output at, is given by,(23)
These equations explicitly show that the beamformer output for
the first source contains the leakage from the second source andthat the output for the second source contains the leakage fromthe first source. Equations (22) and (23) are valid for high SNR,
but they still hold for considerably low SNR, as is shown in our
numerical experiments.
B. Signal-Intensity Reduction and Correlation-Coefficient
Estimation
The leakage from correlated activities not only causes errors
in the time-course estimate but also causes intensity reduction

..
.
..
.

..

(18)

.

and

..
.
..
.

..

.

(19)

1538

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 49, NO. 12, DECEMBER 2002

This equation indicates that the magnitude correlation coefficient can be accurately estimated by substituting the
beamformer outputs into (27). Appendix II shows that themagnitude of the cross coherence in the frequency domain
can also be correctly estimated by using the outputs from thefrequency-domain beamformer.
If a third correlated source exists, the accuracy of the estimated correlation coefficient should be affected by this source.
This influence can be evaluated in the following manner. We assume that the first, second, and third sources are mutually correlated, and their correlation coefficients are denoted as,,. We further assume that the correlation between the firstandand second sources is the target of the measurement and the thirdsource is an interference. According to (17), the beamformer
outputs at the target source locations are

(29)
and

(30)
is the sourcewhere is the location of the third source, andcovariance submatrix related to the correlated three sources.
is presented in Appendix III.)(This
Using (27), (29), (30), and (45), we finally obtain
Fig. 2. (a) Time courses w (t), w (t), and w (t) used for the numericalexperiments. The top panel represents w (t), the middle one, w (t), and thebottom one, w (t). Each time course is normalized by its maximum value. The
broken line in the top panel indicates time course w (t), which is obtained fromw (t) = (1 )w (t) + w (t) with  = 0:5. The correlation coefficient
between w (t) and w (t), , is 0.8 in this case. (b) Examples of the simulated
magnetic-field recordings when SNR is equal to four and  = 0.

0

in the reconstructed source activities. Using (22) and (23) with, we getthe relationship
(24)
(25)

(31)
This equation shows how the values ofandare small,
Clearly, whenandaffectis close to.

.

C. Retrieval of Original Time Courses
As discussed in the preceding sections, the beamformer
contains not only the target source activities but
outputalso the activities correlated with the target sources. However,
if we know the amount of correlated interference contained in, the original target-source time courses can be retrievedby inverting (22) and (23), i.e.,

and
(26)
Equations (25) and (26) indicate that the power of the recon. This reduction
structed sources is reduced by a factor ofof source power has been known as the signal cancellation in
the field of the array signal processing [28], [29].
We define the magnitude correlation coefficient between thefirst and the second sources calculated from the beamformer
output as
(27)

(32)
is the retrieved time course of the th source
whereand . In the right-hand side of the aforemenwhereandtioned equation, can be estimated from (28), andcan be obtained fromand.
IV. NUMERICAL EXPERIMENTS
A. Data Generation

Substituting (24)–(26) into (27), we finally obtain(28)

A series of numerical experiments were conducted to verify
the validity of the arguments presented in Section III. In particular, we have shown that the arguments still hold even for

SEKIHARA et al.: PERFORMANCE OF AN MEG ADAPTIVE-BEAMFORMER TECHNIQUE

h

1539

i

Fig. 3. Magnitude-squared average reconstructions~(r ; t) obtained by using the minimum-variance beamformer in (8) and (11). SNR is equal to 16. (a)  =0:08. (b)  = 0:5. (c)  = 0:6. (d)  = 0:7. (e)  = 0:8. (f)  = 0:95.  is the cross-correlation coefficient between the first and second sources.

Fig. 4. Intensity of the first source with respect to the correlation between thefirst and second sources. The broken line indicates the theoretical relationship1
 . In these Monte Carlo simulations, the mean intensity of the firstsource was calculated from 100 generated data sets. The error bars indicate therange of 2 standard deviations. The SNR was set to four.

p 0

6

low-SNR cases, although the arguments are based on the assumption that noise can be neglected. A sensor alignment of the
37-sensor array from Magnes1 neuromagnetometer was used
in which the sensor coils are arranged in a uniform, concentric
array on a spherical surface with a radius of 12.2 cm. The sensors are first-order axial gradiometers with a baseline of 5 cm.
The coordinate origin was set at the center of the detector coillocated at the center of the coil array. The direction was defined as that from the posterior to the anterior, and the direction was defined as that from the right to the left hemispheres.
Three point sources were assumed to exist on a plane defined as
14D Neuroimaging Inc., San Diego, CA

Fig. 5. Estimated correlation coefficient between the first and secondsources. Monte Carlo-type simulations, which generated 100 sets of magnetic
recordings, were performed and the average values of 
~ for five SNRs areplotted. The error bar shows the 2 standard deviations. The experiments wererepeated for  = 0:8 and  = 0:5. The broken horizontal lines show the true
values of the correlation coefficient.

6. [The values of the spatial coordinates
are expressed in centimeters.] The locations as well as the orientationsof the sources are listed in Table I. The source-sensor configuration and the coordinate system are illustrated in Fig. 1.
,, andwere calcuFirst, three time courses
lated by using

for

(33)

and
(34)

1540

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 49, NO. 12, DECEMBER 2002

Fig. 6. Time-course outputs obtained by using the minimum-variance beamformer for the three sources when SNR is equal to four and (a)  = 0:08, (b)  = 0:4,(c)  = 0:6, and (d)  = 0:8. Time courses from the first to third sources are shown from top to bottom. The solid lines indicate the beamformer outputs and thebroken lines indicate the original time courses.

where , , , andare the numerical parameters controlling the shapes of the time courses. The values of these parameters used in these experiments are listed in Table II. Time,, and
, are shown by the solid lines incoursesFig. 2(a). The correlation coefficients of these three time courses
10 between
and,2
10 between
are 8
and
, and 4 10 between
and.A
was obtained frommodified time course, where parameter controls the degree of theand. The time course
correlation between
when
is shown by the broken line in Fig. 2(a). The
andis 0.8 in thiscorrelation coefficient between
case. The simulated magnetic recordings were calculated by as,, andto the first, second, and third
signingsources, respectively.
The magnetic-field recordings were calculated at 1-ms
intervals from zero to 400 ms by using the spherically ho.
mogeneous conductor [17] with its center set at
White Gaussian noise was then added to the simulated fieldrecordings. The simulated recordings with an SNR of four andwith set at zero are shown in Fig. 2(b). Here, the SNR is

defined as the ratio of the Frobenius norm of the signal-magnetic-field data matrix to that of the noise matrix, i.e., the ratioof. As we explained in Section V,
the SNR of four represents a very poor SNR situation.

B. Signal-Cancellation Experiments
We generated six data sets with the correlation coefficients
between the first and second sources, , of 0.08, 0.5, 0.6, 0.7,
0.8, and 0.95. The SNR was set at 16. The weight vectors for
detecting the , , and components were obtained from (8),and the 3-D reconstruction was performed by using (11). The
reconstruction region included the three sources and an area deandon plane;
fined bythis area is displayed in Fig. 3. The interval between the reconstruction grids was 0.1 cm. The results of magnitude-squared
average reconstruction for the six values of are shown in thesix contour maps in Fig. 3, where each contour line represents a. The reconrelative value of the source magnitudestruction results contain three localized sources at the locations

SEKIHARA et al.: PERFORMANCE OF AN MEG ADAPTIVE-BEAMFORMER TECHNIQUE

1541

exceeds 0.6, and 60% of the original source intensity is still
maintained when the correlation reaches 0.8. Therefore, as far
as the signal intensity reduction is concerned, no serious influences arise from sources with weak or medium degrees of correlation. Even when sources are significantly cor–
, the beamformer can still reconstruct
relatedsuch sources if their intensities are strong enough to overcome
the signal cancellation. The results in Fig. 3 demonstrated thisfact, i.e., there is no large difference between the reconstructedand.
sources when
C. Estimation of Cross-Correlation CoefficientFig. 7. (a) Predicted time courses and the minimum-variance beamformer
outputs of the first and the second sources. The solid lines indicate thebeamformer outputs obtained when  = 0:8 and SNR is equal to four. The
broken lines indicate the time courses predicted from (22) or (23). (b) Resultsof the experiment on time-course retrieval when  = 0:8 and SNR is equal
to four. Solid lines indicate retrieved time courses and broken lines indicate
original time courses assumed.

given in Table I. These results, however, are somewhat blurry asa result of noise effects [30].
These contour maps clearly show that, compared with the intensity of the third source, the intensities of the first and thesecond sources are reduced according to their degree of correlation. This intensity reduction is theoretically predicted in (25)and (26), which indicate that the intensities of the reconstructed. Monte Carlo-typesources are reduced by a factor ofexperiments were performed to check whether this relationshipholds in nonidealistic situations such as when SNR is considerably low. One hundred sets of simulated magnetic recordings
were generated with different noise realizations for ten different
values of . The mean intensity of the reconstructed first sourceis plotted in Fig. 4. The SNR was set at four. The theoreticalis plotted with the broken line, and the error
trendbars show the range of 2 standard deviations. In Fig. 4, although a small discrepancy is observed when approaches one,the broken line is generally overlapped with the plots from theMonte Carlo experiments, indicating that the theoretical relationship holds well even at the low SNR situation.
, the intensity reduction
According to the trend
is less than twenty percent, unless the correlation coefficient

The correlation coefficient between the first and secondsources was estimated from (27) using beamformer outputs
and; these outputs are obtained using (8) and(11). To investigate the influence of noise on the accuracy ofthe calculated correlation coefficient, we performed Monte
Carlo-type experiments in which 100 values of the correlation
coefficient between the first and second sources were calculatedusing the 100 sets of simulated magnetic recordings generated
with different noise realizations. Such experiments wererepeated for five values of SNR, and the results are plotted inFig. 5 for of 0.5 and 0.8. In this figure, the error bars show therange between 2 standard deviations. This plot shows that the
influence of the noise on the estimated correlation coefficientis generally very small. The estimated correlation coefficient isbiased and the average of the estimates is slightly smaller than
the true value when SNR equals three. However, this bias isvery small and is less than 5% for the whole SNR range used
for the experiments.
D. Errors in Time-Course Estimate and Original Time Course
RetrievalReconstructed time courses of the three source activities areshown in Fig. 6. The SNR was set at four in this experiment. The,
results in Fig. 6(a)–(d), respectively, correspond to,, and. The solid lines indicate thebeamformer outputs obtained from (8) and (11), and the brokenlines indicate the original time courses assumed in the numerical experiments. When is small, original time courses can be
retrieved as beamformer outputs. However, as increases, thereconstructed time courses become distorted because of the interference from the other correlated source. This distortion isand; it is small but still disevident when.
cernible when
These distorted beamformer outputs can be theoretically predicted from (22) and (23). The predicted beamformer outputs
are shown in Fig. 7(a). The
and the actual outputs for
outputs predicted from (22) and (23) are shown by broken lines,
and the beamformer outputs are shown by solid lines. The solidand broken lines are nearly perfectly overlapped, indicating that(22) and (23) are valid even for an SNR as low as four. We next
performed a time-course retrieval experiment. In this experiment, the true time courses were retrieved by using the beamandwith (32). The results offormer outputs

1542

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 49, NO. 12, DECEMBER 2002

Fig. 8. Results from the eigenspace-projection beamformer given by (9), (11), and (13). The same simulated recordings as used for obtaining the results in Fig. 6
were again used. Time course outputs for the three sources when SNR is equal to four and (a)  = 0:08, (b)  = 0:4, (c)  = 0:6, and (d)  = 0:8. Time courses
from the first to third sources are shown from top to bottom. The solid lines indicate the beamformer outputs and the broken lines indicate the original time courses
assumed.

this experiment are shown in Fig. 7(b). In this figure, the retrieved and true time courses are nearly perfectly overlapped, although a small discrepancy exists between the two time courses
due probably to noise effects.
We also applied the eigenspace-projection beamformer tothe same data sets used for obtaining the results in Fig. 6.
The weight vector was obtained from (8) and (13), and thereconstruction was performed by using (11). Reconstructed
time courses of the three source activities are shown in Fig. 8,
in which the results in (a)–(d) correspond respectively to
,
,
, and
. These time courses
are, in principle, identical to those in Fig. 6, although the noise
was considerably reduced by the eigenspace projection [26].
The predicted beamformer outputs and the actual outputs forare shown in Fig. 9(a), and the retrieved and true
time courses are shown in Fig. 9(b). In Fig. 9(a), the predicted
outputs are nearly perfectly overlapped with the actual outputs.
In Fig. 9(b), the retrieved and true time courses are overlapped.
These results confirm that the prediction using (22) and (23)and the time-course retrieval using (32) are also effective forthe eigenspace-projection beamformer.
V. DISCUSSIONS
We have shown that when two sources are correlated, thecross-correlation coefficient can be accurately estimated fromthe beamformer outputs. However, if a third correlated source
exists, the accuracy of the estimated correlation coefficient isaffected by this source. This influence can be evaluated from, this equation can be simplified as(31). When(35)(
) is the correlation of the target sources
where(the first and second sources) with the interference source (thethird source). This equation is plotted for the three values of thein Fig. 10. This figure shows that eventarget correlation
when an additional correlated source exists, the errors due toSEKIHARA et al.: PERFORMANCE OF AN MEG ADAPTIVE-BEAMFORMER TECHNIQUE

Fig. 9. (a) Predicted time courses and eigenspace-projection beamformer
outputs of the first and the second sources. Solid lines indicate the beamformer
outputs obtained when  = 0:8 and SNR is equal to four. Broken lines indicate
the time courses predicted from (22) and (23). (b) Results of the experimenton time course retrieval when  = 0:8 and SNR is equal to four. Solid lines
indicate retrieved time courses and broken lines indicate original time courses.

Fig. 10. Effects of a third correlated source on the estimates of the target
correlation coefficient. Dotted, broken, and solid lines correspond to theestimate when  = 0:8,  = 0:7, and  = 0:6, respectively.

this interference source are less than 10%, unlessexceeds. Therefore, it is generally true that if theeighty percent ofcorrelation with the interference source is not as strong as thecorrelation between the target sources, a reasonably accurate
target-source correlation coefficient can be obtained.

1543

One may notice that, in Figs. 6 and 8, the output for the firstsource has a time course that is very similar to. Thissimilarity can be explained as follows. The case ofis considered here. In this case, the time course of the firstwas derived from.
sourcein (22) is equal to 0.53, theSince the coefficient ofis expressed as
output. This indicates that the beamformer reconcompostruction coincidentally removed nearly all the
becomes very close to. On the othernents andis equal to 1.2, thehand, since the coefficient of
output of the second-source time course is equal to. Therefore, output
is considerably distorted by the existence of thecomponent.
In our numerical experiments in Section IV, we used an SNR
of four to show that the results of our analysis in Section III,
which is based on a high-SNR assumption, are still valid for
low-SNR cases. In the following, we show that an SNR of four
is a reasonable assumption for simulating low-SNR situations.
As pointed out in Section II-A, noise with external origins suchas background brain activities is considered as a part of the signals in our analysis. Thus, the noise referred to here indicates
only that with internal origins such as noise from the probe coils,
SQUIDs, and the associated electronics. The level of such noises
in typical modern neuromagnetic instruments is known to begenerally less than 5 fT Hz [1], [31]. On the other hand, thepeak amplitude of a typical neuromagnetic signal is known to
be 50–500 fT [1]. A simple case, where the signal has a sinusoidal waveform with a peak amplitude of 200 fT, is considered here. The bandwith of the data acquisition is assumed to
be 200 Hz, so the SNR defined from the Frobenius norm is cal. This SNR can be improved byculated as
signal averaging for evoked neuromagnetic measurements witha factor equal to the square-root of the number of averagings. It
can, therefore, be concluded that an SNR of four is almost equal
to the SNR of the evoked neuromagnetic measurements with avery few number of averagings, and this SNR value represents
very low-SNR situations.
One influence of correlated sources that is not discussed inthis paper is the degradation of the spatial resolution. The spatial resolution is degraded by the source correlation, but thedegradation is generally small when the sources have a weak or
medium degree of correlation, as demonstrated in Fig. 3. Thatis, an image blur can be observed, but it is not significant in). We are currently investigating the facFig. 3(a)–(e) (tors that determine the spatial resolution in the adaptive beamformer neuromagnetic reconstruction. The results of these investigations, including the influence of the source correlation,
will be published in the near future.
In summary, we analyzed the effects of correlated sources on
the reconstruction results of MEG adaptive beamformer techniques. Two major influences were found: signal-intensity reductions and distortions in the time-course estimate. Our numerical experiments showed that the theoretical relationship between signal intensity and source correlation holds in nonidealistic situations such as when SNR is considerably low. Indeed,
as far as the signal intensity is concerned, no serious influences

1544

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 49, NO. 12, DECEMBER 2002

arise from sources with weak or medium degrees of correlation). However, distortion in time-course estimates is dis(cernible even for such correlated sources. We also developed amethod of correcting the time-course distortion for two-source
correlation, and our numerical experiments demonstrate that the
distortion can be almost perfectly corrected.
APPENDIX IThis appendix shows the derivation of (16), which was origas
inally derived by Zoltowski [27]. We define the matrix. Then, using (5), we get the relationship. Therefore, when the SNR is so high that noise can beholds, and we getneglected, the relationshipas a vector, whereis the Fourier spectrum of the th-channel recording.
We also define the Fourier transform of the source-activity
whereis
vector as. Let us
the Fourier spectrum of the th source magnitude:define the cross-spectrum matrix of the measured data as, where the superscript indicates the
Hermitian transpose. The estimated Fourier spectrum of the, is obtained assource activity at ,(38)
where the weight vector

is given by
(39)

(36)
indicates the pseudoinverse ofdefined as. Regarding (36), we also assume thatis nonsingular, and we use the fact that the transpose and the pseudoinverse operations are interchangeable. We apply the formula, whereis thevector that has all elements equal to zero except for the th element that is equal toone. Using this formula and (36), we obtain

Here,

When the first and second sources have a coherence value
at the frequency , we get the relationship
(40)
(41)
is the power spectrum of the th source at . Then,
where
substituting (40) and (41) into

(37)
Note that this equation exactly holds when noise can be, (37) is also valid for theneglected. Since
eigenspace-projection beamformer mentioned in Section II-D.
One vital assumption in the argument above is that of the
low-rank signal; i.e., the assumption that the number of signal
sources is less than the number of sensors . If this assumpis not equal to and (37) cannot betion does not hold,derived. However, since the number of sensors in commerciallyavailable neuromagnetometers has increased constantly over thelast ten years and the latest neuromagnetometers are equippedwith 200–300 sensor channels, the low-rank-signal assumption
should hold for most neuromagnetic measurements.

(42)

leads to(43)
where is the estimated magnitude coherence. Equation (43)indicates that the accurate estimate of the magnitude coherence
can be obtained by substituting the outputs of the beamformer
into (42). It should be noted that the analysis in this appendix
provides a theoretical validation of the method recently proposed by Gross et al. [15].

APPENDIX IIAPPENDIX IIIThis appendix shows that the magnitude of the cross coherence in the frequency domain can also be estimated using theoutputs from the frequency-domain beamformer implementation. We define the Fourier transform of the measurement vector

This appendix presents the inverse of the source covariancematrix in the case of a three-source correlation. We assume thatthe first, second, and third sources are mutually correlated. The,,
correlation coefficients of these sources are denoted as(45)

SEKIHARA et al.: PERFORMANCE OF AN MEG ADAPTIVE-BEAMFORMER TECHNIQUE

and. In the source-activity covariance matrix, the submatrix related to the three correlated sources is denoted as .
This matrix is expressed as(44)

Its inverse is expressed as shown in (45) at the bottom of theindicates the determinant ofandprevious page, whereis given by(46)

REFERENCES[1] M. Hämäläinen, R. Hari, R. J. IImoniemi, J. Knuutila, and O. V.
Lounasmaa, “Magnetoencephalography—Theory, instrumentation, andapplications to noninvasive studies of the working human brain,” Rev.
Mod. Phys., vol. 65, pp. 413–497, 1993.
[2] T. P. L. Roberts, D. Poeppel, and H. A. Rowley, “Magnetoencephalography and magnetic source imaging,” Neuropsych., Neuropsych., Behav.
Neurology, vol. 11, pp. 49–64, 1998.
[3] B. D. van Veen and K. M. Buckley, “Beamforming: A versatile approachto spatial filtering,” IEEE Acoust. Speech Signal Processing Mag., vol.
5, pp. 4–24, Apr. 1988.
[4] S. E. Robinson and D. F. Rose, “Current source image estimation by
spatially filtered MEG,” in Biomagnetism Clinical Aspects, M. Hoke etal., Eds. New York: Elsevier, 1992, pp. 761–765.
[5] S. E. Robinson and J. Vrba, “Functional neuroimaging by synthetic aperture magnetometry (SAM),” in Recent Advances in Biomagnetism, T.
Yoshimoto et al., Eds. Sendai, Japan: Tohoku Univ. Press, 1999, pp.
302–305.
[6] J. Gross and A. A. Ioannides, “Linear transformations of data space in
MEG,” Phys. Med. Biol., vol. 44, pp. 2081–2097, 1999.
[7] K. Sekihara, S. S. Nagarajan, D. Poeppel, A. Marantz, and Y. Miyashita,“Reconstructing spatio–temporal activities of neural sources using an
MEG vector beamformer technique,” IEEE Trans. Biomed. Eng., vol.
48, pp. 760–771, July 2001.
[8] I. Hashimoto, T. Kimura, Y. Iguchi, R. Takino, and K. Sekihara,“Dynamic activation of distinct cytoarchitectonic areas of the human
SI cortex after median nerve stimulation,” NeuroReport, vol. 12, pp.
1891–1897, 2001.
[9] I. Hashimoto, K. Sakuma, T. Kimura, Y. Iguchi, and K. Sekihara,“Serial activation of distinct cytoarchitectonic areas of the human SI
cortex after posterior tibial nerve stimulation,” NeuroReport, vol. 12,
pp. 1857–1862, 2001.
[10] M. S. Hämäläinen and R. J. Ilmoniemi, “Interpreting measured magnetic
fields of the brain: Estimates of current distributions,” Helsinki Univ.
Technol., Helsinki, Finland, Tech. Rep. TKK-F-A559, 1984.
[11] K. Sekihara and S. S. Nagarajan, “Neuromagnetic source reconstruction and inverse modeling,” in Modeling and Imaging of Bioelectric Activity—Principles and Applications, B. He, Ed. Norwell, MA: Kluwer,
to be published.
[12] K. Sekihara and B. Scholz, “Generalized Wiener estimation of threedimensional current distribution from biomagnetic measurements,” in
Biomag 96: Proc. 10th Int. Conf. Biomagnetism, C. J. Aine et al., Eds.,
1996, pp. 338–341.
[13] A. S. Gevins, S. L. Bressler, N. H. Morgan, B. A. Cutillo, R. M. White, D.
S. Greer, and J. Illes, “Event-related covariances during a bimanual visuomotor task. I. Methods and analysis of stimulus- and response-locked
data,” Electroenceph. Clin. Neurophysiol., vol. 74, pp. 58–75, 1989.
[14] J. Sarnthein, H. Petsche, P. Rappelsberger, G. L. Shaw, and A. von Stein,“Synchronization between prefrontal and posterior association cortex
during human working memory,” Proc. Nat. Acad. Sci., vol. 95, pp.
7092–7096, 1998.

1545[15] J. Gross, J. Kujara, M. Hämäläinen, L. Timmermann, A. Schnitzler, and
R. Salmelin, “Dynamic imaging of coherent sources: Studying neuralinteractions in the human brain,” Proc. Nat. Acad. Sci., vol. 98, pp.
694–699, 2001.
[16] W. Singer, “Neuronal synchrony: A versatile code for the definition ofrelations,” Neuron, vol. 24, pp. 49–65, 1999.
[17] J. Sarvas, “Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem,” Phys. Med. Biol., vol. 32, pp. 11–22, 1987.
[18] A. Paulraj, B. Ottersten, R. Roy, A. Swindlehurst, G. Xu, and T. Kailath,“Subspace methods for directions-of-arrival estimation,” in Handbookof Statistics, N. K. Bose and C. R. Rao, Eds. Amsterdam, The Netherlands: Elsevier, 1993, pp. 693–739.
[19] R. O. Schmidt, “A signal subspace approach to multiple emitter locationand spectral estimation,” Ph.D. dissertation, Dept. Elect. Eng., Stanford
Univ., Stanford, CA, 1981.
[20] L. L. Scharf, Statistical Signal Processing: Detection, Estimation, andTime Series analysis. Reading, MA: Addison-Wesley, 1991.
[21] D. H. Johnson and D. E. Dudgeon, Array Signal Processing: Concepts
and Techniques. Englewood Cliffs, NJ: Prentice-Hall, 1993.
[22] J. Capon, “High-resolution frequency wavenumber spectrum analysis,”
Proc. IEEE, vol. 57, pp. 1408–1419, Aug. 1969.
[23] 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, June 1992.
[24] M. E. Spencer, R. M. Leahy, J. C. Mosher, and P. S. Lewis, “Adaptive
filters for monitoring localized brain activity from surface potential timeseries,” in Conf. Rec. 26th Ann. Asilomer Conf. Signals, Systems, Computers, Nov. 1992, pp. 156–161.
[25] B. D. van Veen, W. van Drongelen, M. Yuchtman, and A. Suzuki, “Localization of brain electrical activity via linearly constrained minimum variance spatial filtering,” IEEE Trans. Biomed. Eng., vol. 44, pp. 867–880,
Sept. 1997.
[26] K. Sekihara, S. S. Nagarajan, D. Poeppel, A. Marantz, and Y. Miyashita,“Application of an MEG eigenspace beamformer to reconstructing
spatio–temporal activities of neural sources,” Human Brain Mapping,vol. 15, pp. 199–215, 2002.
[27] M. D. Zoltowski, “On the performance analysis of the MVDR beamformer in the presence of correlated interference,” IEEE Trans. Signal
Processing, vol. 36, pp. 945–947, June 1988.
[28] B. Widrow, K. M. Duval, R. P. Gooch, and W. C. Newman, “Signal
cancellation phenomena in adaptive antennas: Causes and cures,” IEEE
Trans. Antennas Propagat., vol. AP-30, pp. 469–478, May 1982.
[29] A. Paulraj, V. U. Reddy, and T. Kailath, “Analysis of signal cancellation
due to multipath in optimum beamformers for moving arrays,” IEEE J.
Ocean. Eng., vol. OE-12, pp. 163–172, Jan. 1987.
[30] H. Cox, “Resolving power and sensitivity to mismatch of optimum arrayprocessors,” J. Acoust. Soc. Amer., vol. 54, pp. 771–785, 1973.
[31] J. Vrba, “SQUID gradiometers in real environments,” in SQUID
Sensors: Fundamentals, Fabrication and Applications, H. Weinstock,
Ed. Norwell, MA: Kluwer, 1996, pp. 117–178.

Kensuke Sekihara received the M.S. and Ph.D. degrees from the Tokyo Institute of Technology, Tokyo,
Japan, in 1976 and 1987, respectively.
From 1976 to 2000, he was with Central ResearchLaboratory, Hitachi, Ltd., Tokyo, Japan. He was avisiting Research Scientist at Stanford University,
Stanford, CA, from 1985 to 1986, and at Basic Development, Siemens Medical Engineering, Erlangen,Germany, from 1991 to 1992. From 1996 to 2000, he
was with the “Mind Articulation” Research Projectsponsored by the Japan Science and Technology
Corporation. He is currently a Professor at Tokyo Metropolitan Instituteof Technology, Tokyo, Japan. His research interests include biomagnetic
inverse problems, and statistical estimation theory, especially its application tofunctional neuroimaging.
Dr. Sekihara is a Member of the IEEE Medicine and Biology Society and theIEEE Signal Processing Society.

1546

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 49, NO. 12, DECEMBER 2002

Srikantan S. Nagarajan received the B.S. degree inelectrical engineering from the University of Madras,
Madras, India, and the M.S. and Ph.D. degrees fromthe Department of Biomedical Engineering, Case
Western Reserve University, Cleveland, OH, in 1993
and 1995, respectively.
From 1995 to 1998, he was a Postdoctoral Fellow
at the Keck Center for Integrative Neuroscience, theUniversity of California, San Francisco (UCSF). In
1999, he was a full-time Research Scientist at Scientific Learning Corporation, Berkeley, CA, and an Adjunct Assistant Professor in the Department of Otolaryngology at UCSF. Subsequently, he was an Assistant Professor in the Department of Bioengineering
at the University of Utah, Salt Lake City. Currently, he is the Director of the
Biomagnetic Imaging Laboratory and an Assistant Professor in the Departmentof Radiology at UCSF. His research interests in neural engineering include bioelectromagnetism, systems and computational neuroscience, rehabilitation, andstatistical signal processing.

David Poeppel received the B.S. and Ph.D. degrees
in cognitive neuroscience from the Massachusetts
Institute of Technology, Cambridge, in 1990 and
1995, respectively.
Currently, he is an Assistant Professor in the
Departments of Linguistics and Biology at the
University of Maryland at College Park. His researchuses the functional neuroimaging methods magnetoencephalography (MEG), electroencephalography(EEG), positron emission tomography (PET), andfunctional 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, and the
Ph.D. degree in linguistics from the Massachusetts
Institute of Technology (MIT), Cambridge, in 1978
and 1981, respectively.
He joined the Faculty of MIT in 1990, where heis currently Professor of Linguistics and Head of the
Department of Linguistics and Philosophy. His research interests include the syntax and morphology ofnatural languages, linguistic universals, and the neurobiology of language. He is currently involved inrevising morphological theory within linguistics and in exploring MEG techniques to uncover how the brain processes language.