䉬

Human Brain Mapping 15:199 –215(2002) 䉬
DOI 10.1002/hbm.10019

Application of an MEG Eigenspace Beamformer
to Reconstructing Spatio-Temporal Activities of
Neural Sources
Kensuke Sekihara,1* Srikantan S. Nagarajan,2 David Poeppel,3Alec Marantz,4 and Yasushi Miyashita51Department of Electronic Systems and Engineering, Tokyo Metropolitan Institute of Technology,
Tokyo, Japan
2
Bioengineering Department, University of Utah, Salt Lake City, Utah
3
Department of Linguistics and Biology, University of Maryland, College Park, Maryland
4
Department of Linguistics and Philosophy, Massachusetts Institute of Technology,
Cambridge, Massachusetts
5
Department of Physiology, The University of Tokyo, School of Medicine, Hongo, Bunkyo-ku, Tokyo, Japan䉬䉬

Abstract: We have applied the eigenspace-based beamformer to reconstruct spatio-temporal activities ofneural sources from MEG data. The weight vector of the eigenspace-based beamformer is obtained byprojecting the weight vector of the minimum-variance beamformer onto the signal subspace of a measurement covariance matrix. This projection removes the residual noise-subspace component that considerably degrades the signal-to-noise ratio (SNR) of the beamformer output when errors in estimating thesensor lead field exist. Therefore, the eigenspace-based beamformer produces a SNR considerably higher
than that of the minimum-variance beamformer in practical situations. The effectiveness of the eigenspace-based beamformer was validated in our numerical experiments and experiments using auditory
responses. We further extended the eigenspace-based beamformer so that it incorporates the information
regarding the noise covariance matrix. Such a prewhitened eigenspace beamformer was experimentallydemonstrated to be useful when large background activity exists. Hum. Brain Mapping 15:199 –215, 2002.
© 2002 Wiley-Liss, Inc.

Key words: magnetoencephalography; biomagnetism; MEG inverse problems; beamformer; functional
neuroimaging; neuromagnetic signal processing䉬䉬

INTRODUCTION

Grant sponsor: Japan Science and Technology Corporation; Grant
sponsor: Kayamori Foundation of Informational Science Advancement; Grant sponsor: the Whitaker Foundation.
*Correspondence to: Kensuke Sekihara, PhD, Tokyo Metropolitan
Institute of Technology, Asahigaoka 6-6, Hino, Tokyo 191-0065,
Japan. E-mail: ksekiha@cc.tmit.ac.jp
Received for publication 30 March 2001; accepted 17 October 2001

© 2002 Wiley-Liss, Inc.

Among the various kinds of functional neuroimaging methodologies, magnetoencephalography (MEG)
has the the major advantage that it can provide finetime resolution of the millisecond order [Hämäläinen
et al., 1993]. Neuromagnetic imaging can thus be used
to visualize neural activities with such a fine time
resolution and to provide functional information
about brain dynamics [Roberts et al., 1998]. Toward䉬 Sekihara et al. 䉬this goal, a number of algorithms for reconstructing
spatio-temporal source activities have been investigated.
In our study, we explore the possibility of applying
a class of techniques referred to as the adaptive beamformer to this spatio-temporal reconstruction of theneural-source activities. The adaptive beamformer
provides a versatile form of spatial filtering, and it has
been originally developed in the fields of array signal
processing, including radar, sonar, and seismic exploration [van Veen and Buckley, 1988]. One well-known
technique of this kind, the minimum-variance beamformer, has already been successfully applied to solve
the MEG/EEG source-localization problem [Grossand Ioannides, 1999; Robinson and Vrba, 1999; Spencer et al., 1992; van Veen et al., 1997]. We found,however, that the minimum-variance beamformer
generally is very sensitive to errors in the forward
modeling or errors in estimating the data covariancematrix when it is applied to the reconstruction ofsource activities at each instant in time. Because such
errors are almost inevitable in neuromagnetic measurements, the minimum-variance beamformer generally provides noisy spatio-temporal reconstruction results, as is demonstrated later in this study.
One technique has been developed to overcome thepoor performance of the minimum-variance-based
technique caused by the forward modeling errors[Cox et al., 1987] and covariance-matrix estimation
errors [Carlson, 1988]. The technique, referred to as
the diagonal loading, uses the regularized inverse of
the measurement covariance matrix, instead of its direct matrix inverse, when calculating the weight vector. Although this technique has been applied to the
MEG source localization problem [Gross and Ioannides, 1999; Robinson and Vrba, 1999], it has been
known that the regularization leads to a trade-off between the spatial resolution and the SNR of the beamformer output. The results from our experiments demonstrate this trade-off relationship.
We propose to apply the eigenspace-based beamformer [Feldman and Griffiths, 1991; van Veen, 1988]
to reconstructing spatio-temporal activities of neuralsources. The weight vector of the eigenspace-based
beamformer is obtained by projecting the weight vector of the minimum-variance beamformer onto thesignal subspace of the measurement covariance matrix. This projection removes a residual noise-subspacecomponent that considerably degrades the signal-tonoise ratio (SNR) of the beamformer output when theabove-mentioned errors exist. Therefore, the eigenspace-based beamformer produces a SNR significantly
higher than that from the minimum-variance beam-䉬

former; thus, it is a suitable method for spatio-temporal reconstruction of the neural source activities.
After a brief introduction to the minimum-variancebeamformer, we describe the principle of the eigenspace-based beamformer. We also present a theoreticalcomparison between these beamformer techniques
and numerical experiments that validate our arguments. We also describe the application of the beamformer techniques to auditory-evoked responses. The
results demonstrate that the eigenspace-based beamformer techniques are highly effective for the spatiotemporal reconstruction of source activities.
MATERIALS AND METHODSDefinitions and Fundamental Relationships
Let us define the magnetic field measured by themth detector coil at time t as bm(t), and a column
vector b(t) ⫽ [b1(t), b2(t),…, bM(t)]T as a set of measured
data where M is the total number of detector coils andthe superscript T indicates the matrix transpose. A
spatial location is represented by a 3D vector r:r⫽ (x,y,z). A total of Q current sources are assumed togenerate the neuromagnetic field, and the locations ofthese sources are denoted as r1, r2,…, rQ. The moment
magnitude of the qth source at time t is denoted as s(rq,t), and the source magnitude vector is defined as s(t)⫽ [s(r1, t), s(r2, t)…, s(rQ, t)]T. The orientation of the qthsource is defined as a 3D column vector ␩(rq, t) ⫽ [␩x(rq, t), ␩y (rq, t), ␩z (rq, t)]T whose ␨ component (where␨ equals x, y, or z in this study) is equal to the cosineof the angle between the direction of the source moment and the ␨ direction. We assume that the orientation of each source is time independent. Omitting thetime notation t, we define a 3Q ⫻ Q matrix that
expresses the orientations of all Q sources as ⌿ such
that

⌿⫽

冤

␩共r1 兲
0···0·
␩共r2 兲 •
0
·
·
·
·
··
0
•
·
·
0
···
0 ␩共rQ兲冥.

Let us define l␨m (r) as the mth sensor output inducedby the unit-magnitude source located at r and directedin the ␨ direction. The column vector l␨(r) is defined asl␨(r) ⫽ [l ␨1 (r), l ␨2 (r),. . .,l ␨M (r)] T . We define the lead field
matrix, which represents the sensitivity of the wholesensor array at r, as L(r) ⫽ [lx(r), ly(r), lz(r)]. We define,for later use, the lead-field vector in the source-moment direction as l(r); it is obtained by using l(r)200 䉬䉬 MEG Eigenspace Beamformer 䉬⫽ L(r)␩(r). The composite lead field matrix for theentire set of Q sources is defined as
L c ⫽ 关L共r1 兲, L共r2 兲, . . . , L共rQ兲兴.

(1)

The relationship between b(t) and s(t) is then expressed asb共t兲 ⫽ 共Lc⌿兲s共t兲 ⫹ n共t兲,(2)where n(t) is the additive noise.
Let us define the measurement covariance matrix as
Rb; i.e., Rb ⫽ 具b(t)bT(t)典, where 具  典 indicates the ensemble average (this ensemble average is usually replaced with the time average over a certain time window). Let us also define the covariance matrix of thesource-moment activity as Rs; i.e., Rs ⫽ 具s(t)sT(t)典.
Then, using equation (2), we get the relationship between the measurement covariance matrix and thesource-activity covariance matrix such thatR b ⫽ 共Lc⌿兲Rs共⌿TLcT兲 ⫹ ␴ 2 I,

(3)where the noise in the measured data is assumed to be
the white Gaussian noise with the variance of ␴2 andI is the unit matrix.
Let us define the jth eigenvalue and eigenvector of
Rb as ␭j and ej, respectively. Unless some source activities are perfectly correlated with each other, therank of Rs is equal to the number of sources Q. Therefore, according to equation (3), Rb has Q eigenvalues
greater than ␴2 and M–Q eigenvalues that are equal to␴2. Let us define the matrices ES and EN as ES ⫽[e1,…,eQ] and EN ⫽ [eQ ⫹ 1, …, eM]. The column span of
ES is the maximum-likelihood estimate of the signal
subspace of Rb and the span of EN is that of the noise
subspace [Scharf, 1991]. Using equation (3), it can beshown that, at source locations, the lead field matrix
with the correct source orientation is orthogonal to thenoise subspace of the measurement covariance matrix[Schmidt, 1981], i.e.,关L共r兲 ␩共r兲兴 TEN ⫽ lT共r兲EN ⫽ 0.

ŝ共r, t兲 ⫽ wT共r兲b共t兲,

where ŝ(r,t) is the estimated magnitude of the sourcemoment at r and time t. In this equation, w(r) is acolumn vector characterizing the filter weight. Notethat because this weight vector is calculated for any
spatial location r, the source-moment distribution canbe reconstructed by scanning the output of the spatial
filter over a region of interest in a perfectly postprocessing manner. One well-known spatial filter ofthis kind is the minimum-variance distortion-lessbeamformer originally developed for seismic-arraysignal processing [Capon, 1969]. In this technique, thefilter weight vector wm(r) is obtained by minimizing
wTm (r)Rb wm (r) under the constraint of lT(r)wm(r) ⫽ 1.
The explicit form of the weight vector for the minimum-variance beamformer is known to beRb⫺ 1 l共r兲w m 共r兲 ⫽ T.
l 共r兲Rb⫺ 1 l共r兲To estimate the source moment, we have focused onusing the class of techniques referred to as the spatial
filter. The spatial filter techniques use the following
simple linear operation for estimating the source moment,䉬(6)This minimum variance beamformer has been one of
the most popular spatial-filter techniques in various
signal-processing fields. It has also been applied toneuromagnetic source localization [Robinson andVrba, 1999; Spencer et al., 1992; van Drongelen et al.,
1996; van Veen et al., 1997].
Eigenspace-Based Beamformer
The eigenspace-based beamformer [Feldman and
Griffiths, 1991; van Veen, 1988; Yu and Yeh, 1995] provides an output SNR much higher than that of the minimum-variance beamformer in practical situations. Let
us decompose the measurement covariance matrix Rb
into its signal and noise subspace components; i.e.,R b ⫽ ES⌳SEST ⫹ EN⌳NENT.

(7)

Here, we define the matrices ⌳S and ⌳N as⌳ S ⫽ diag关␭ 1 , . . . , ␭ Q兴and⌳N ⫽ diag关␭ Q ⫹ 1 , . . . , ␭ M兴,(4)

Minimum-Variance Beamformer

(5)

(8)

where diag[…] indicates a diagonal matrix whose diagonal elements are equal to the entries in the parenthesis. Using equations (6) and (7), and defining ␮⫽ 1/[lT (r)Rb⫺ 1 l(r)], we express the weight vector forthe minimum-variance beamformer as

201 䉬

w m 共r兲 ⫽ ␮Rb⫺ 1 l共r兲 ⫽ ␮⌫Sl共r兲 ⫹ ␮⌫Nl共r兲,(9)䉬 Sekihara et al. 䉬

where

SNR 共MV兲 ⫽⌫ S ⫽ ES⌳S⫺ 1 EST,and⌫N ⫽ EN⌳N⫺ 1 ENT.

In equation (9), the second term on the right side,␮⌫Nl(r), should ideally be equal to zero because thelead-field vector l(r) is orthogonal to EN at the source
locations as indicated by equation (4). Various factors,
however, prevent this term from being zero, and a
non-zero ␮⌫Nl(r) seriously degrades SNR as explainedin the next section. Therefore, the eigenspace-based
beamformer uses only the first term of equation (9) to
calculate its filter weight vector we(r); i.e.,⌫Sl共r兲
w e 共r兲 ⫽ ␮⌫Sl共r兲 ⫽ T.
l 共r兲Rbl共r兲(10)Note that we(r) is equal to the projection of wm(r) onto
the signal subspace of Rb. The following relationshipholds [Feldman and Griffiths, 1991]:
w e 共r兲 ⫽ ESESTwm共r兲.

(11)Comparison Between Minimum-Variance andEigenspace Beamformers
Although the minimum-variance beamformer ideally has exactly the same SNR as that of the eigenspace-based beamformer, the SNR of the eigenspace
beamformer is significantly higher in practical applications, as demonstrated in the following sections. The
reason of this high SNR can be understood as follows.
Let us assume that a signal source with a moment
magnitude equal to s(t) exists at r. Then, the estimateof s(t), ŝ(t), is derived by ŝ(t) ⫽ wTb(t) ⫽ wTl(r)s(t) andthe average power of ŝ(t), P̂s, is expressed asP̂ s ⫽ 具ŝ共t兲 2 典 ⫽ wTl共r兲具s共t兲 2 典lT共r兲w ⫽ wTl共r兲PslT共r兲w⫽ Ps关lT共r兲w兴 T关lT共r兲w兴,P̂ n ⫽ wT具n共t兲nT共t兲典w ⫽ ␴ 2 wTw,

(13)

where we again assume that 具n(t)nT(t)典 ⫽ ␴2I. Thus,
the output SNR of the minimum-variance beamformer, SNR(MV), is expressed as [Chang and Yeh,
1992, 1993],

䉬

(14)

The SNR for the eigenspace-based beamformer,
SNR(ES), is obtained as
SNR 共ES兲 ⫽P̂ s Ps 关lT共r兲⌫Sl共r兲兴 2.
⫽
P̂ n ␴ 2 关lT共r兲⌫S2 l共r兲兴

(15)

The only difference between equations, namely (14)and (15) is the existence of the second term lT (r)⌫ 2N l(r)in the denominator of the right-hand side of equation(14). It is readily apparent that SNR(MV) and SNR(ES)are equal if we can use an accurate noise subspaceestimate and an accurate lead-field vector, because theterm lT (r)⌫ 2N l(r) is exactly equal to zero in this case. Itis, however, generally difficult to attain the relationship, lT (r)⌫ 2N l(r) ⫽ 0. One obvious reason for thisdifficulty is that when calculating ⌫ 2N in practice, instead of using Rb, the sample covariance matrix R̂bmust be used; R̂b is calculated from R̂b ⫽⌺ kK ⫽ 1 b(t k )bT (t k ) where K is the number of time points.
Use of the sample covariance matrix inevitably makeslT (r)⌫ 2N l(r) have a non-zero value and causes SNR of
the minimum-variance beamformer to degrade [Richmond, 1998].
Another factor that is specific to MEG and causeslT (r)⌫ 2N l(r) to have a non-zero value is that it is almost
impossible to use a perfectly accurate lead-field vector.
This is because a conductivity distribution in the brain
usually be approximated by using some kind of conductor model, such as the spherically homogeneous
conductor model [Sarvas, 1987], to calculate the leadfield matrix. Although this error may be reduced to acertain extent by using a realistic head model, theerror cannot be perfectly avoided. Let us define theoverall error in estimating l(r) as ⑀. Assuming that
储l(r)储 2 ⬎⬎储ε储 2 , we can get(12)where Ps is the average power of s(t) defined by Ps ⫽具s(t)2典. Similarly, the average noise power contained in
the estimated results, P̂n is given by关lT共r兲⌫Sl共r兲兴 2
P̂ s Ps.
⫽ 2 T
P̂ n ␴ 关l 共r兲⌫S2 l共r兲 ⫹ lT共r兲⌫N2 l共r兲兴

SNR 共MV兲 ⫽关lT共r兲⌫Sl共r兲兴 2
P̂ s
Ps
.
ⱌ 2 T
␴ 关l 共r兲⌫S2 l共r兲 ⫹ ε T⌫N2 ε兴
P̂ n

(16)

Note that, in the denominator of the right-hand side ofthis equation, the norm of the matrix ε T ⌫ 2N ε has anorder of magnitude proportional to 储ε储 2 /␭ 2 , where ␭N
represents one of the noise-level eigenvalues of Rb.
The eigenvalue ␭N is usually significantly smaller than
the signal-level eigenvalues. Therefore, equation (16)indicates that even when the error 储ε储 is very small,the term ε T ⌫ 2N ε may not be negligibly small comparedto the first term in the denominator.

202 䉬䉬 MEG Eigenspace Beamformer 䉬

Extension to Prewhitened Eigenspace
Beamformer
Real-life MEG data often contain interference arising from background brain activities [Sekihara et al.,1997]. Such interference does not affect the final reconstruction results if the sources of interference are spatially well separated from the signal source of interest.
The influence, however, may not be negligible if theinterference source is located close to the signal source(within the range of the spatial resolution). Such interference is known to cause spatially non-white noisein neuromagnetic measurements, and if the information regarding the noise covariance matrix is obtained,
the effect of such interferences can be significantly
reduced by using the prewhitened eigenspace beamformer technique described below.
When the noise is nonwhite, equation (3) becomesR b ⫽ 共Lc⌿兲Rs共⌿ L 兲 ⫹ Rn,
T

T
c

(17)

where Rn is the noise covariance matrix obtained fromRn ⫽ 具n(t)nT(t)典. Let us denote ẽj as an eigenvector
obtained by solving the generalized eigenvalue problem,R b ẽ j ⫽ ␭˜ jRnẽj.

(18)

In this case, it is easy to show that the following
orthogonality relationship holds,
关L共r兲 ␩兴 TẼN ⫽ 0,

(19)

where ẼN ⫽ [ẽQ⬘ ⫹ 1,…,ẽM] and Q⬘ is the number ofsources. Thus, the signal subspace projector can beformed by using ẼS ẼST , where Q̃S ⫽ [ẽ1,…,ẽQ] and theweight vector for the prewhitened eigenspace beamformer, w̃e(r), is obtained fromw̃ e 共r兲 ⫽ ẼSẼSTwm共r兲.

(20)

The effectiveness of the beamformer based on thisequation in removing background interferences isdemonstrated in the following sections.

Figure 1.
The source and detector configuration used in the numericalexperiments. The circle shows the cross section of the sphereused for the forward calculation. The square shows the reconstruction region used for our numerical experiments.

nel Magnes™ biomagnetic measurement system (Biomagnetic Technologies Inc., San Diego) was used inthese experiments. The coordinate origin was definedas the center of the detector coil located at the center ofthe coil alignment. Three signal sources were assumed
to exist on a plane defined as x ⫽ 1.0 cm; their locations were (1, ⫺1, ⫺6), (1, 1, ⫺6), and (1, 1.6, ⫺7.2) cm.
The source and detector configuration is shown schematically in Figure 1. The spherically homogeneous
conductor model with the origin set at (1, 0, ⫺11) (cm)was used. The magnetic field was generated at 1-msec
intervals from ⫺300 to 400 msec. The time courses of
the three source activities assumed in our numericalexperiments are shown in Figure 2. The first and second sources were the sources of interest, and the thirdsource simulated background interference. Gaussian
noise was added to the generated magnetic field sothat the SNR was set at 8 and the SNR was defined by
the ratio of the Frobenius norm of the signal-magneticfield data matrix to that of the noise matrix. The generated 37-channel magnetic-field recordings are also
shown in Figure 2.

NUMERICAL EXPERIMENTS
Spatio-Temporal Reconstruction Experiments
Data Generation
We conducted a series of numerical experiments totest the effectiveness of the eigenspace-based beamformer techniques. The coil alignment of the 37-chan-䉬

The reconstruction region was set as an area definedby ⫺2 ⱕ y ⱕ 2 and ⫺8 ⱕ z ⱕ ⫺5 (cm) (as indicated by
the square in Fig. 1), and the interval of the reconstruction grid was 1 mm in both directions. In our experi-203 䉬䉬 Sekihara et al. 䉬ments, when applying beamformer techniques, wefirst used the method described in the Appendix for
estimating the optimum source orientation ␩ˆ opt andwe then estimated the lead field vector l(r) ⫽L(r)␩ˆ opt at each scanning grid point. To display the
results of the spatio-temporal reconstruction, threetime instants at 220, 268, and 300 msec (marked bythree vertical broken lines in Fig. 2) were selected. The
amplitude of the second source was zero at 220 msec,
all the sources had non-zero amplitudes at 268 msec,
and only the second source had a non-zero amplitudeat 300 msec (Fig. 2). The measurement covariancematrix was calculated by using a time window between 0 – 400 msec. The power, averaged over thistime window, of the second source is 1.5 times stronger than that of the first source.
Results From Minimum-Variance Beamformer
Figure 2.
The moment time courses of the three sources assumed in the
numerical experiments. Time courses from the first to the thirdsources are shown from the top to the third row, respectively.
Each waveform is normalized by its maximum value and the threevertical broken lines show the three time points 220, 268, and 300
msec. The simulated 37-channel magnetic field recordings areshown in the bottom row.

We first applied the minimum-variance beamformer to the generated data set in Figure 2. The
results of the spatio-temporal reconstruction are
shown in Figure 3a. Contour maps in the upper-left,
upper-right and lower-left positions indicate the distributions of the source-moment power ŝ(r, t)2 at thethree time instants of 220, 268, and 300 msec, respectively. The lower-right contour map shows the time-

Figure 3.
The results obtained using the minimum-variance beamformer. (a) structed region is indicated by the square in Figure 1. (b) EstimatedSnapshots of source activities at 220 msec (upper left), 268 msec time courses of the first source (top), the second source (middle),(upper right), and 300 msec (lower left). The time-averaged re- and the third source (bottom). The three vertical broken lines
construction is shown in the lower-right position. The recon- show the three time points 220, 268, and 300 msec.

䉬

204 䉬

䉬 MEG Eigenspace Beamformer 䉬

Figure 4.
The results obtained using the minimum-variance beamformer lower-right position. (b) Estimated time courses of the first sourcewith the regularized inverse. (a) Snapshots of source activities at (top), the second source (middle), and the third source (bottom).
220 msec (upper left), 268 msec (upper right), and 300 msec The three vertical broken lines show the three time points 220,
(lower left). The time-averaged reconstruction is shown in the 268, and 300 msec.

averaged reconstruction, obtained by averaging ŝ(r, t)2
over the time window from 0 – 400 msec. The estimated time courses at the pixels nearest to the threesource locations are shown in Figure 3b. These results
show that the spatio-temporal reconstruction obtainedby using the minimum variance beamformer was
fairly noisy; that is, the estimated time courses contained a considerable amount of noise, and the threesnapshots of the source activity contained the influence of this noisy reconstruction. The time-averaged
results, however, did not contain this influence andclearly detected the three sources.
Results From Minimum-Variance Beamformer
With Regularized Inverse
We next tested the minimum-variance weight vector together with the use of the regularized inverse (Rb⫹ ␥I)⫺1. The regularization parameter was set at0.003␭11 where ␭1 is the largest eigenvalue of Rb. The
results of the reconstruction are shown in Figure 4a,b.
1This value was chosen so that the SNR of the estimated timecourses in Figure 4b is nearly equal to that of the time courses in
Figure 5b, which was obtained by using the eigenspace-based beamformer.

䉬The time course estimation in Figure 4b shows that the
SNR of the beamformer output was greatly improved,
but the contour maps in Figure 4a indicate that a
considerable amount of blur was introduced. Our results here confirm the trade-off relationship between
the spatial resolution and the output SNR.
Results From Eigenspace-Based Beamformer
We applied the eigenspace-based beamformer given
by equation (10) to the same computer-generated data
set. The reconstructed source distributions are shownin Figure 5a, and the estimated time courses areshown in Figure 5b. These figures show that thisbeamformer technique considerably improved theoutput SNR with almost no sacrifice of spatial resolution.
Results From Prewhitened Eigenspace
Beamformer
We tested the effectiveness of the prewhitened eigenspace beamformer proposed earlier in reducing
the background interference; namely the influence
from a third source. First, we estimated the noise
covariance matrix using the prestimulus time window

205 䉬䉬 Sekihara et al. 䉬

Figure 5.
The results obtained using the eigenspace-based beamformer. (a) time courses of the first source (top), the second source (middle),
Snapshots of source activities at 220 msec (upper left), 268 msec and the third source (bottom). The three vertical broken lines(upper right), and 300 msec (lower left). The time-averaged re- show the three time points 220, 268, and 300 msec.
construction is shown in the lower-right position. (b) EstimatedFigure 6.
The results obtained using the prewhitened eigenspace-based (b) Estimated time courses of the first source (top), the second
beamformer. (a) Snapshots of source activities at 220 msec (upper source (middle), and the third source (bottom). The three verticalleft), 268 msec (upper right), and 300 msec (lower left). The broken lines show the three time points 220, 268, and 300 msec.
time-averaged reconstruction is shown in the lower-right position.

䉬

206 䉬

䉬 MEG Eigenspace Beamformer 䉬

stimuli were presented to a subject’s right ear. The
sensor array was placed above a subject’s left hemisphere with the position adjusted to optimally record
the N1m auditory-evoked field. The inter-stimulus interval randomly varied between 1.75 sec and 2.25 sec,
and the average interval was 2 sec. The sampling
frequency was set at 1 kHz. An on-line filter with abandwidth from 1– 400 Hz was used, and no postprocessing digital filter was applied. To express theresults of reconstructing source activities in this section, we used the head coordinates illustrated in Figure 7.

Figure 7.
The x, y, and z coordinates used to express the reconstruction
results in the application to auditory-evoked responses. The midpoint between the left and right pre-auricular points is defined as
the coordinate origin. The axis directed away from the origin
toward the left pre-auricular point is defined as the ⫹y axis andthat from the origin to the nasion is the ⫹x axis. The ⫹z axis isdefined as the axis perpendicular to both these axes and is directed from the origin to the vertex.

Results From Minimum-Variance andEigenspace Beamformer
The auditory stimulus was a 1,000-Hz pure tone,
and we measured 256 epochs. The auditory-evoked
response averaged across these epochs is shown between 0 –300 msec in Figure 8a. We applied the eigenspace-based beamformer to reconstruct the spatio-temporal source activities from this auditory data set. The

from ⫺300 – 0 msec. We then applied the prewhitened
eigenspace beamformer (equation (20)) with this estimated noise covariance matrix to reconstruct thesource-moment distribution (Fig. 6a). The estimatedtime courses are shown in Figure 6b. These results
show that the influence of the third source was nearlycompletely removed. Also, the results in Figure 6b
show that the time courses of the first and the secondsources were not affected by the prewhitening procedure.
APPLICATION TO AUDITORY-EVOKEDRESPONSES
Data acquisition Condition
To further demonstrate the superiority of the eigenspace-based beamformer over the minimum-variance beamformer, we applied these beamformer techniques to auditory-evoked responses. Auditoryevoked fields were measured by using the 37-channel
Magnes™ magnetometer installed at the Biomagnetic
Imaging Laboratory, University of California, San
Francisco. Healthy male volunteers participated in the
MEG measurements. Their written informed consents
were obtained and the MEG measurements were approved by the Committee on Human Research, University of California, San Francisco. All measurements
were done in a magnetically-shielded room. Auditory䉬

Figure 8.
(a) Thirty-seven channel recordings of auditory evoked magnetic
field measured using a 1-kHz pure tone. A total of 256 epochs
were averaged. Three vertical broken lines indicate time instants
at 87, 90, 101 msec, at which the snapshots of the source activity
are shown in Figures 9a, 10a, and 11a. (b) Eigenvalue spectrum of
Rb obtained from the time window between 0 and 300 msec of thedata set shown in (a).

207 䉬䉬 Sekihara et al. 䉬

Figure 9.
Results of applying the eigenspace-based beamformer to the data directions are shown. The upper-case letters L and R indicates theshown in Figure 8. (a) Reconstructed source magnitude distribu- left and right hemispheres. (b) Estimated time course obtained attions at 87 msec (top row), 90 msec (middle row), and 101 msec a location where the intensity in the results in (a) has the maxi(bottom row). The maximum-intensity projections onto the axial mum value.
(left column), coronal (middle column), and sagittal (right column)
䉬

208 䉬

䉬 MEG Eigenspace Beamformer 䉬

Figure 10.
Results of applying the the minimum-variance beamformer to the msec (bottom row). (b) Estimated time course obtained at adata shown in Figure 8. (a) Reconstructed source magnitude location where the intensity in the results in Figure 9a has the
distributions at 87 msec (top row), 90 msec (middle row), and 101 maximum value.

䉬

209 䉬

䉬 Sekihara et al. 䉬

Figure 11.
Results of applying the the minimum-variance beamformer with msec (middle row), and 101 msec (bottom row). (b) Estimatedthe regularized inverse to the data shown in Figure 8. (a) Recon- time course obtained at a location where the intensity in the
structed source magnitude distributions at 87 msec (top row), 90 results in (a) has the maximum value.

䉬

210 䉬

䉬 MEG Eigenspace Beamformer 䉬time window ranging from 0 –300 msec was used for
calculating the covariance matrix Rb. The eigenvalue
spectrum of Rb is shown in Figure 8b. The dimensionof the signal subspace was set at 1 because the eigenvalue spectrum contained one distinctly large eigenvalue. Three time instants, 87, 90, 101 msec, wereselected near the peak vertex of N100m. These instants
are shown by three vertical broken lines in Figure 8a.
The reconstructed source-magnitude maps at thesethree latencies are shown in Figure 9a. These threesnapshots of the source activities contain a single,
sharply-localized activity in the left temporal-lobe,
probably near the primary auditory cortex area. The
time course in the maximum point in these contour
maps is shown in Figure 9b. The time course has aclear negative peak near the latency of 50 msec and alarge positive peak near the latency of 100 msec.
We next applied the minimum-variance beamformer to the same data set, and the results are shownin Figure 10a,b. These results are very noisy. The
snapshots at 87 msec and 90 msec contain false sourcesin addition to the source near the primary auditory
area. The time course of the auditory source in Figure
10b is so noisy that the N100m can hardly be identified. We then applied the minimum-variance beamformer with the regularized inverse. The results areshown in Figure 11a,b. Here, the regularization parameter ␥ was set at 0.04␭1. This value was chosen sothat the SNR of the time course in Figure 11b wasnearly equal to that of the time course in Figure 9b.
The results in Figures 9 –11 demonstrate that the regularization causes a spatial blur although it improves
the SNR in the time course estimate.
Results From Prewhitened Beamformer
We applied the prewhitened beamformer to the data
set shown in Figure 12a. The data set was one of the
100-epoch selectively averaged results obtained in a series of syllable-discrimination experiments. The auditory
stimuli used in the experiments were four kinds of syllables /dae/, /bae/, /pae/, and /tae/. The subject was
asked to discriminate the voiced syllables /dae/ and/bae/ from the voiceless syllables /pae/ and /tae/. The
subject pressed one response button when perceiving a/dae/ or /bae/ and pressed another button when perceiving a /pae/ or /tae/. Stimuli were presented to thesubject’s right ear, and the sensor array was placed
above the subject’s left hemisphere. The subject used hisleft fingers to press the response buttons. The four syllables were presented in a pseudo random order at a
variable inter-stimulus interval ranging from 1–1.5 sec.
The averaged response for the particular syllable /bae/䉬Figure 12.
Thirty-seven channel recordings of an auditory magnetic fieldevoked by the voiced syllable /bae/. A total of 100 epochs wereaveraged. The two vertical broken lines indicate time instants at
78, 101 msec, at which the snapshots of the source activity are
shown in Figures 13 and 14. (b) Prewhitened eigenvalue spectrumof Rb obtained from the time window between 0 and 200 ms ofthe data set shown in (a).

are shown in Figure 12a. The activity in the motor areaelicited by button pressing was roughly time locked tothe stimulus and thus the influence from such a motor
activity was contained to some extent in this auditory
response despite that the subject used his epsi-lateralfingers.
We applied the prewhitened eigenspace beamformer to test its effectiveness in removing the influence of the motor activities. The measurement covariance matrix was calculated from the time window
between 0 –200 msec, and the noise covariance matrix
was calculated from the time window between⫺100 – 0 msec (the whole prestimulus portion). The
dimension of the signal subspace was set at 1 because
the prewhitened eigenvalue spectrum (shown in Fig.
12b) contained one large eigenvalue. The reconstructed source-magnitude maps for two time instants, 78 and 101 msec, are shown in Figure 13a; thesetime instants are shown by two vertical lines in Figure
12a. At both time instants, the activation in the primary auditory area was clearly detected by the proposed prewhitened beamformer. The results from the211 䉬

䉬 Sekihara et al. 䉬

Figure 13.
(a) Results of applying the prewhitened eigenspace-based beam- pace beamformer to the same data. Reconstructed source magformer to the data shown in Figure 12a. Reconstructed source nitude distributions at 78 msec (upper row) and 101 msec (botmagnitude distributions at 78 msec (upper row) and 101 msec tom row).
(bottom row). (b) Results of applying the non-prewhitened eigens-non-prewhitened eigenspace-based beamformer are
shown in Figure 13b for comparison. The source reconstruction at 101 msec contains a single localized
activity in the left primary auditory area. At the latency of 78 msec, however, the non-prewhitened
beamformer failed to detect the auditory source due,probably, to the existence of the background motor
activity.
DISCUSSION
In our numerical experiments the results of the spatio-temporal reconstruction obtained by using the䉬minimum variance beamformer was fairly noisy. In
these experiments, the spherically homogeneous conductor model was used for both procedures: the generation of the simulated magnetic field and the forward calculation in the reconstruction. Therefore, thecause of such noisy reconstruction was primarily because of the use of the sample covariance matrix R̂b. In
our auditory experiments, because the spherically homogeneous conductor model was used for the forward calculation, the lead field matrix was not perfectly accurate and the primary cause of the noisy
reconstruction is probably the use of such an inaccu-

212 䉬䉬 MEG Eigenspace Beamformer 䉬

Figure 14.
Results of applying the prewhitened eigenspace-based beamformer to the data shown in Figure 12awith the signal subspace dimension set at 3.

rate lead field in addition to the use of the sample
covariance matrix. Therefore, these auditory experiments had very noisy results. A comparison among
Figures 9 –11 suggests that the use of the eigenspacebased beamformer is highly effective for the spatiotemporal reconstruction of source activities.
We emphasize, though, that the effectiveness of theeigenspace-based beamformer becomes evident only
when we perform spatio-temporal reconstruction. For
time-averaged reconstruction, the results obtained using either the minimum variance or the eigenspacebased beamformer were more or less the same (asshown by the time-averaged reconstruction in Figs. 3
and 5). This is probably the reason why the poor
performance of the minimum-variance beamformer
has been somewhat overlooked in the previous investigations [Robinson and Vrba, 1999; van Drongelen etal., 1996; van Veen et al., 1997]. In these investigations,the spatio-temporal reconstruction of source activities
was not emphasized; instead, a time-averaged reconstruction of the source activities was obtained.
The signal subspace dimension Q in equation (8) isin principle determined by separating distinctly large
eigenvalues from the small eigenvalues. This separation, however, may not be easy if there is no clear
threshold in the eigenvalue spectrum. For example,
the results in Figure 13a were obtained by setting thesignal subspace dimension to one. This determination,
however, is somewhat ambiguous because, as can beseen in Figure 12b, there are two more eigenvalues
that are much smaller than the first eigenvalue but still䉬

greater than the noise-level eigenvalues. The results
obtained by setting the signal-subspace dimension at
three are shown in Figure 14. The comparison between
Figures 13a and 14 shows that when the signal subspace dimension was determined differently, the final
results were not significantly different.
In general, the overestimation of the signal subspacedimension gives intermediate results between thosefrom the minimum-variance beamformer and thosefrom the eigenspace beamformer with a correct signalsubspace dimension. Let us consider the case where
the signal subspace dimension is overestimated at Q⫹ ⌬. Then, the SNR of such an eigenspace beamformer, SNR⌬ is expressed as,
SNR ⌬ ⫽

关lT共r兲⌫Sl共r兲兴 2
P̂ s Ps,⫽ 2 T
P̂ n ␴ 关l 共r兲⌫S2 l共r兲 ⫹ ε T⌫2⌬ ε兴

(21)

where
⌫ ⌬ ⫽ 关eQ ⫹ 1 , . . . , eQ ⫹ ⌬ 兴⫻冋␭Q ⫹ 1 · · ·
0
·
·
··
·
·
·
·
·
0
· · · ␭Q ⫹ ⌬册

⫺1关eQ ⫹ 1 , . . . , eQ ⫹ ⌬ 兴 T

2Because the relationship ε T ⌫ ⌬ε ⬍ ε T ⌫ 2N ε holds, the relationship SNR(ES) ⬎ SNR⌬ ⬎ SNR(MV) always holds
2true. Moreover, if ε T ⌫ ⌬ε is small compared to ε T ⌫ S2 ε ,the eigenspace-based beamformer with an overesti-

213 䉬䉬 Sekihara et al. 䉬mated signal-subspace dimension gives nearly thesame SNR as that from the eigenspace-based beamformer with a correct signal-subspace dimension. The
results in Figure 14 exactly show this case.
In summary, we have applied the eigenspace-based
beamformer to reconstruct spatio-temporal activitiesof neural sources. This beamformer attains a SNR
significantly higher than that of the minimum-variance beamformer, particularly when errors in estimating the sensor lead field exist. We further extended theeigenspace beamformer so that it incorporates the information regarding the noise spatial correlation. The
effectiveness of these eigenspace-based beamformer
techniques was validated in our numerical experiments and experiments using auditory responses.
ACKNOWLEDGMENTS
The authors wish to thank Ms. Susanne Honma andDr. Timothy Roberts for their help in performing auditory MEG measurements. This work was carried out
as part of the MIT-JST International Cooperative Research Project “Mind Articulation”. The work was also
funded by a grant from the Whitaker Foundation to
SN.
REFERENCES
Capon J (1969): High-resolution frequency wavenumber spectrum
analysis. Proc IEEE 57:1408 –1419.
Carlson BD (1988): Covariance matrix estimation errors and diagonal loading in adaptive arrays. IEEE Trans Aerospace Electronic
Syst 24:397– 401.
Chang L, Yeh CC (1992): Performance of DMI and eigenspace-based
beamformers. IEEE Trans Antenn Propagat 40:1336 –1347.
Chang L, Yeh CC (1993): Effect of pointing errors on the performance of the projection beamformer. IEEE Trans Antenn Propagat 41:1045–1056.
Cox H, Zeskind RM, Owen MM (1987): Robust adaptive beamforming. IEEE Trans Signal Process 35:1365–1376.
Feldman DD, Griffiths LJ (1991): A constrained projection approachfor robust adaptive beamforming. Proceedings of the International Conference on Acoustics, Speech, and Signal Processing,
Toronto. p 1357–1360.
Ferrara ER, Parks TM (1983): Direction finding with an array ofantennas having diverse polarizations. IEEE Trans Antenn
Propagat 31:231–236.
Gross J, Ioannides AA (1999): Linear transformations of data spacein MEG. Phys Med Biol 44:2081–2097.
Hämäläinen M, Hari R, IImoniemi RJ, Knuutila J, Lounasmaa OV(1993): Magnetoence-phalography-theory, instrumentation, andapplications to noninvasive studies of the working human brain.
Rev Mod Phys 65:413– 497.
Li J (1994): On polarization estimation using a crossed-dipole array.
IEEE Trans Signal Process 42:977–983.
Mosher JC, Lewis PS, Leahy RM (1992): Multiple dipole modelingand localization from spatio-temporal MEG data. IEEE Trans
Biomed Eng 39:541–557.

䉬

Mosher JC, Leahy RM (1998): Recursive music: a framework for
EEG and MEG source localization. IEEE Trans Biomed Eng
45:1342–1354.
Richmond CD (1998): Response of sample covariance based MVDR
beamformer to imperfect look and inhomogeneities. IEEE Signal
Process Let 5:325–327.
Roberts TPL, Poeppel D, Rowley HA (1998): Magnetoencephalography and magnetic source imaging. Neuropsychiatry Neuropsychol Behav Neurol 11:49 – 64.
Robinson SE, Vrba J (1999): Functional neuroimaging by synthetic
aperture magnetometry (SAM). In: Yoshimoto T, et al., editors.
Recent advances in biomagnetism. Sendai: Tohoku University
Press. p 302–305.
Sarvas J (1987): Basic mathematical and electromagnetic concepts ofthe biomagnetic inverse problem. Phys Med Biol 32:11–22.
Scharf LL (1991): Statistical signal processing: detection, estimation,and time series analysis, New York: Addison-Wesley Publishing
Co. 252 p.
Schmidt RO (1981): A signal subspace approach to multiple emitter
location and spectral estimation. PhD thesis. Stanford, CA: Stanford University.
Sekihara K, Scholz B (1996); Generalized Wiener estimation of 3D
current distribution from biomagnetic measurements. In: Aine
CJ, et al., editors. Biomag 96: Proceedings of the Tenth International Conference on Biomagnetism. New York: Springer-Verlag. p 338 –341.
Sekihara K, Poeppel D, Marantz A, Koizumi H, Miyashita Y (1997):
Noise covariance incorporated MEG-MUSIC algorithm: amethod for multiple-dipole estimation tolerant of the influenceof background brain activity. IEEE Trans Biomed Eng 44:839 –847.
Spencer ME, Leahy RM, Mosher JC, Lewis PS (1992): Adaptive
filters for monitoring localized brain activity from surface
potential time series. Conference record for 26th Annual Asilomer Conference on Signals, Systems, and Computers. p
156 –161.
van Drongelen W, Yuchtman M, van Veen BD, van Huffelen AC(1996): A spatial filtering technique to detect and localize multiple sources in the brain. Brain Topogr 9:39 – 49.
van Veen BD (1988): Eigenstructure based partially adaptive array
design. IEEE Trans Antenn Propagat 36:357–362.
van Veen BD, Buckley KM (1988): Beamforming: a versatile approach to spatial filtering. IEEE ASSP Magazine 5:4 –24.
van Veen BD, van Drongelen W, Yuchtman M, Suzuki A (1997):
Localization of brain electrical activity via linearly constrained
minimum variance spatial filtering. IEEE Trans Biomed Eng
44:867– 880.
Yu JL, Yeh CC (1995): Generalized eigenspace-based beamformers.
IEEE Trans Signal Process 43:2453–2461.

APPENDIX
This appendix describes the method of obtaining areasonable estimate of the source orientation. Thismethod was used in our experiments. The source orientations can be obtained by utilizing the orthogonality relationship in equation (4). The optimum estimateof ␩ is one that minimizes the following measure ofthis orthogonality;214 䉬䉬 MEG Eigenspace Beamformer 䉬

J共 ␩, r兲 ⫽

储ENL共r兲 ␩储 2 关 ␩TLT共r兲兴ENENT关L共r兲 ␩兴⫽.
储L共r兲 ␩储 2
关 ␩TLT共r兲兴关L共r兲 ␩兴(22)This minimization problem can be solved by using thegeneralized-eigenproblem formulation. The optimum
orientation ␩ˆ opt(r) satisfies the relationship,关L T 共r兲ENENTL共r兲兴 ␩ˆ opt共r兲 ⫽ ␭ min关LT共r兲L共r兲兴 ␩ˆ opt共r兲,(23)

where ␭min indicates the minimum eigenvalue of thisgeneralized eigenvalue problem. Therefore, the optimum estimate of the source-moment orientation is
obtained by finding the eigenvector corresponding to
the minimum eigenvalues of equation (23).
When noise is non-white and has a known covariance matrix, the orthogonality relationship is ex-䉬pressed in equation (19). The optimum orientation
␩ˆ opt(r) should therefore be obtained solving
关L T 共r兲ẼNẼNTL共r兲兴 ␩ˆ opt共r兲 ⫽ ␭˜ min关LT共r兲Rn⫺ 1 L共r兲兴 ␩ˆ opt共r兲.
(24)That is, the eigenvector corresponding to the miminum eigenvalue ␭˜ min in the above equation gives the
optimum estimate of the source-moment orientation.
The method presented in this appendix was originallydeveloped for estimating the polarization of wavefronts received by an array of antennas in radar-signal
processing [Ferrara and Parks, 1983; Li, 1994; Schmidt,
1981], and its application to estimating an orientationof neuromagnetic sources has been reported [Mosher
et al., 1992; Mosher and Leahy, 1998; Sekihara and
Scholz, 1996].

215 䉬