1358

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 57, NO. 6, JUNE 2010

Array-Gain Constraint Minimum-Norm Spatial Filter
With Recursively Updated Gram Matrix For
Biomagnetic Source Imaging
Isamu Kumihashi and Kensuke Sekihara∗ , Fellow, IEEE

Abstract—This paper proposes a novel spatial filter for biomagnetic source imaging. The proposed spatial filter is derived basedon a modified version of the minimum-norm spatial filter andis designed to have a performance close to that of the adaptive
minimum-variance spatial filter through the use of an estimatedcovariance matrix. In this method, the theoretical form of the measurement covariance matrix is estimated as an updated gram matrix in a recursive procedure. Since the proposed method does not
use the sample covariance matrix, it is free of the well-known weaknesses of the minimum-variance spatial filter, namely, the proposed
spatial filter does not require a large number of time samples, andit can even be applied to single-time-sample data. It is also robust
to source correlation. We have validated the method’s effectivenessby our computer simulations as well as through experiments using
auditory-evoked magnetoencephalographic data.
Index Terms—Biomagnetic source reconstruction, magnetoencephalography (MEG), minimum-norm method, minimumvariance spatial filter, source reconstruction, spatial filter.

I. INTRODUCTION
HE SPATIAL filter is a popular method used to reconstruct
a source distribution from bioelectromagnetic data [1]. It isa linear operator that applies a set of linear weights to measured
data to obtain a source reconstruction. When the weight of a
spatial filter only depends on the geometry of the measurements,
such a spatial filter is referred to as a nonadaptive spatial filter,
which includes the minimum-norm-based filters discussed in
Section II-B. In contrast to nonadaptive spatial filters, the weightof an adaptive spatial filter depends not only on the measurement
geometry but also on the measurement covariance matrix. A
representative adaptive spatial filter is the minimum-variancespatial filter described in Section II-C.
It is generally true that a nonadaptive spatial filter has a large
source location bias, particularly when applied to 3-D volume
reconstruction [1], [2]. In contrast, the source location bias of an
adaptive spatial filter is generally small and can provide reasonable 3-D source reconstruction [1], [2]. The spatial resolution

T

Manuscript received July 31, 2009; revised November 15, 2009; accepted
December 16, 2009. Date of publication February 17, 2010; date of current
version May 14, 2010. This work was supported by Grants-in-Aid from the
Ministry of Education, Science, Culture and Sports in Japan (C20500394). Asterisk indicates corresponding author.
I. Kumihashi is with the Department of Systems Design and Engineering,Tokyo Metropolitan University, Tokyo 191-0065, Japan (e-mail: kumihashi136_
ac@yahoo.co.jp).
∗ K. Sekihara is with the Department of Systems Design and Engineering,Tokyo Metropolitan University, Tokyo 191-0065, Japan (e-mail: ksekiha@
cc.tmit.ac.jp).
Color versions of one or more of the figures in this paper are available online
at http://ieeexplore.ieee.org.
Digital Object Identifier 10.1109/TBME.2010.2040735

of an adaptive spatial filter is significantly higher than that of
nonadaptive spatial filters [1], [2].
However, adaptive spatial filters have well-known weaknesses. First, since the computation of their weights requiresa sample covariance matrix of the measured data, a large number of time samples are needed to obtain an accurate sample
covariance matrix. In many cases of bioelectromagnetic measurements, however, it is rather difficult to obtain a large numberof time samples. Second, the adaptive spatial filter is known tobe sensitive to the source correlation, and it generally fails to reconstruct source activities when they are highly correlated [3],[4]. Typical cases include the bilateral activation of auditory
cortices [5].
In this paper, we propose a novel spatial filter based on amodified version of the minimum-norm filter, called the arraygain constrained minimum-norm filter. The proposed method isdesigned to estimate the theoretical form of the measurement
covariance matrix through recursively updating the gram matrix. Thus, the proposed method is termed “array-gain constraint
minimum-norm filter with recursively updated gram matrix”
abbreviated as the AGMN-RUG spatial filter in this paper. The
proposed AGMN-RUG method is able to provide 3-D volume
reconstruction with a spatial resolution considerably higher thanthat of existing minimum-norm-based methods. Nonetheless,the method is free of the aforementioned weaknesses occurring with adaptive spatial filters, because it does not dependon the sample covariance matrix. The method is also robust
to the source correlation, and this robustness can be seen in
our computer simulation as well as in our experiments using
auditory-evoked magnetoencephalographic (MEG) data.
After presenting reviews of the nonadaptive minimum-norm
filters and the adaptive minimum-variance filters in Section II,
this paper proposes a novel AGMN-RUG spatial filter in
Section III. The proposed spatial filter is validated by our computer simulation in Section IV and by applications to auditoryevoked MEG data presented in Section V. Throughout thispaper, plain italics indicate scalars, lower case boldface italics indicate vectors, and upper case boldface italics indicate
matrices.
II. FORMULATION OF REPRESENTATIVE SPATIAL FILTERS
A. Spatial Filter Source Reconstruction
Let us define the biomagnetic field measured by the mth sensor at time t as bm (t), and the measured data as the column vectorb(t) = [b1 (t), b2 (t), . . . , bM (t)]T . Here, M is the total number

0018-9294/$26.00 © 2010 IEEE
Authorized licensed use limited to: Univ of Calif San Francico. Downloaded on June 24,2010 at 04:31:39 UTC from IEEE Xplore. Restrictions apply.

KUMIHASHI AND SEKIHARA: ARRAY-GAIN CONSTRAINT MINIMUM-NORM SPATIAL FILTERof sensors and the superscript T indicates the matrix transpose. A
spatial location is represented by a 3-D vector r, as r = (x, y, z).
A source vector at r and time t is defined as a 3-D column vector s(r, t), where s(r, t) = [sx (r, t), sy (r, t), sz (r, t)]T , andsx (r, t), sy (r, t), and sz (r, t) are the x, y, and z componentsof the source vector s(r, t).
xy(r), lm(r), andWe denote the outputs of the mth sensor as lmzlm (r) when the unit-magnitude source at r is directed in thex, y, and z directions, respectively. Then, the column vectorslx (r), ly (r), and lz (r) are defined asili (r) = [l1i (r), l2i (r), . . . , lM
(r)]T

(1)

where i = x, y or z. The vectors, lx (r), ly (r), and lz (r), express
the sensor array sensitivity for a source at r directed in the x,y, and z directions. The lead-field matrix is expressed as the
M × 3 matrix:
L(r) = [lx (r), ly (r), lz (r)].

(2)

This matrix L(r) represents the sensitivity of the whole sensor array for a source at r. The relationship between the measurement vector b(t) and the 3-D source vector s(r, t) is thenexpressed asL(r)s(r, t) d3 r
(3)b(t) =Ω

where d3 r indicates the volume element, and the integral is
performed over the entire space within which sources may exist.
This space is referred to as the source space and denoted Ω inthis paper.
Using the measurement vector b(t), the spatial filter reconstructs the source vector s(r, t) usings(r, t) = [sx (r, t), sy (r, t), sz (r, t)]T = W T (r)b(t).

beam responses to have an ideal delta-function-like shape [1],i.e, the weight matrix of the minimum-norm spatial filter isderived using
W (r) = arg min F(W )(5)where M × 1 weight vectors wx (r), wy (r), and wz (r),respectively detect sx (r, t), sy (r, t), and sz (r, t). This weightmatrix characterizes the properties of the spatial filter. Combining (3) and (4), we derive (omitting explicit time notation t)
W T (r)L(r  )s(r  ) d3 r .
(6)s(r) =Ω

Here, W T (r)L(r  ) is called the resolution kernel, andexpresses the relationship between the true source distributions(r) and the estimated source distribution s(r). One way ofinterpreting the resolution kernel is to consider it as a function ofr  using a fixed r. In such a case, W T (r)L(r  ), called the beam
response, expresses the sensitivity of the spatial filter pointing
at r to a source located at r  . In other words (when r = r  ), thebeam response represents the gain on unwanted leakage signals
from sources located elsewhere from the filter pointing location.
B. Minimum-Norm-Based Spatial Filters
Next, we show that the minimum-norm-based spatial filters
can be derived using optimization formulations that force their

(7)

W (r)

where





tr{ W T (r)L(r  ) − δ(r − r  )I

F(W ) =Ω



W T (r)L(r  ) − δ(r − r  )IT

} d3 r . (8)Here, I is the identity matrix, tr{·} indicates the matrix trace,and δ(r) is the delta function. The cost function F(W ) in (8)indicates the total power of the leakage signals. Therefore, using
the previous optimization, we wish to derive a spatial filter thatpasses the signal from the pointing location r but minimizes thetotal power of the leakages from all other locations.
To find the W (r) that minimizes this F(W ), we calculate
the derivative ∂F(W )/∂W (r) and set it to zero, i.e.,
∂F(W )L(r  )LT (r  )W (r) − L(r  )δ(r − r  ) d3 r
=2∂W (r)
Ω


T


3 
=2
L(r )L (r ) d r W (r) − 2L(r)Ω= 2 [GW (r) − L(r)] = 0(9)where G, called the gram matrix, is an M × M matrix definedas
L(r)LT (r) d3 r.
(10)
G=
Ω

From (9), we can derive the weight matrix given by(4)

where s(r, t) is the estimated source vector at r and time t. In(4), the weight matrix W (r) is an M × 3 matrix defined as
W (r) = [wx (r), wy (r), wz (r)]1359

W (r) = G−1 L(r).

(11)

The spatial filter in (11) is equal to the well-known minimumnorm spatial filter [6].
Several variants of the minimum-norm filter can be derived
by adding some constraints to the optimization formulation in(7). An example of such constraint is one on the filter gain at
the pointing location. We first consider the unit-gain constraint
expressed as W T (r)L(r) = I. The constrained optimization
is explicitly formulated, in this case, as
W (r) = arg min F(W ),W (r)subject toW T (r)L(r) = I(12)where F(W ) is defined in (8). The weight matrix satisfying(12) is known to have the following form:
W (r) = G−1 L(r)[LT (r)G−1 L(r)]−1 .

(13)This spatial filter is called the unit-gain constraint minimumnorm filter reported in [7].
The constrained optimization,
W (r) = arg min F(W ), subject to


W (r)

 T

T
W (r)L(r  ) W T (r)L(r  ) d3 r = IΩ

Authorized licensed use limited to: Univ of Calif San Francico. Downloaded on June 24,2010 at 04:31:39 UTC from IEEE Xplore. Restrictions apply.

(14)

1360

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 57, NO. 6, JUNE 2010

leads to the weight matrix expressed as
W (r) = G−1 L(r)[LT (r)G−1 L(r)]−1/2 .

(15)This spatial filter is known as standardized low-resolution electromagnetic tomography (sLORETA) [8].
C. Minimum-Variance Spatial Filter
The minimum-variance spatial filter is a representative adaptive spatial filter, whose weight matrix is formulated using themeasurement covariance matrix R. Using the unit-gain constraint, W T (r)L(r) = I, the weight matrix is derived using
the optimization,
W (r) = arg min tr{W T RW },W (r)subject to W T (r)L(r) = I.

(17)

This minimum-variance spatial filter is known to have significantly higher spatial resolution and much smaller source location bias, compared to the minimum-norm-based spatial filters
mentioned in the preceding subsection [2].
For biomagnetic imaging, the norm of the lead field L(r)
is spatially nonuniform and the minimum-variance spatial filter
is affected by this nonuniformity of the lead-field norm, as canbe seen in the weight expression in (17). Particularly when thespherical homogeneous conductor model [9] is used for computing the lead field, the source reconstruction results containa false intensity increase around the center of the sphere [1],because L(r) becomes zero and the weight becomes infinity
at the center.
When L(r) is spatially nonuniform, instead of using theunit-gain constraint W T (r)L(r) = I, it is more reasonable touse the constraint W T (r)L(r) = Λ, where00lx (r)
(18)
0 .
0
ly (r)
Λ=
0
lz (r)
In this case, the weight matrix is expressed asT

W (r) = R−1 L(r)[L (r)R−1 L(r)]−1Ω

where · indicates averaging. This expression for R is used inderiving the spatial filter proposed in the next section.

(16)

The resultant weight matrix is expressed as [3],W (r) = R−1 L(r)[LT (r)R−1 L(r)]−1 .

 is used for R,
these weights, a sample covariance estimate R
and a large number of time samples are generally needed toobtain an accurate sample covariance matrix. However, in manybioelectromagnetic applications, it is difficult to obtain a large
number of time samples. In the following section, we develop
a spatial filter that outperforms the minimum-norm-based filters without using the sample covariance matrix. We derive an
expression for R needed in our arguments in the next section.
Using (3), the measurement covariance matrix is expressed as
TL(r)s(r, t)sT (r, t)LT (r)dr (21)R = b(t)b (t) =III. PROPOSED AGMN-RUG SPATIAL FILTER
The proposed spatial filter is based on the array-gainconstraint minimum-norm filter, which is derived using the following optimization:
W (r) = arg min F(W ),W (r)W T (r)L(r) = Λ

subject to(22)where Λ is the gain matrix defined in (18), and the cost function
F(W ) is defined, in this case, as
tr{ W T (r)L(r  ) − δ(r − r  )Λ P (r  )F(W ) =Ω

T
 T
W (r)L(r  ) − δ(r − r  )Λ } d3 r

(23)



where P (r ) is a 3 × 3 positive semidefinite matrix. This P (r  )is a matrix weighting function that controls location-dependent
weighting when the resolution kernel is optimized. The resultant
spatial filter is derived as
−1

T

−1

W (r) = Ḡ L(r)[L (r)Ḡ L(r)]−1
where

(24)


L(r)P (r)LT (r) d3 r.

Ḡ =

(25)

Ω

(19)

where L(r) is the lead-field matrix consisting of the normalized
columns, i.e.,lx (r)ly (r)lz (r),,
L(r) =.
(20)
lx (r) ly (r) lz (r)In (19), the weight is independent of the norm of the lead field,
and we can avoid the artifacts due to the nonuniformity of thelead-field norm. Since lx (r), ly (r), and lz (r) represent
the gains of the sensor array for a source directed in the x-, y-,and z-directions, the constraint W T (r)L(r) = Λ is called thearray-gain constraint, and the spatial filter in (19) is called thearray-gain (constraint) minimum-variance filter.
The weight matrices of the minimum-variance filters in (17)and (19) contain the covariance matrix R. When computing

If we compare (24) with (19), it is the weight equation for thearray-gain constraint minimum-variance filter. The only difference between (24) and (19) is the use of Ḡ, instead of R. We
may then compare Ḡ [see (25)] and R [see (21)] and find thatthese two matrices are very similar and differ only in the manner of weighting in the integral, i.e., the covariance matrix R isexpressed as the integral of L(r)LT (r), with the weighting of
the source power matrix s(r, t)sT (r, t). On the other hand,
the gram matrix Ḡ is given by the integral of L(r)LT (r), withthe weighting of an arbitrary matrix function P (r).
Therefore, if we could use a P (r) that is similar to thesource power matrix s(r, t)sT (r, t), the gram matrix Ḡ
could approximate the covariance matrix R, and consequently,
the performance of the array-gain minimum-norm filter would
approach the performance of the minimum-variance filter.
However, since the true source distribution s(r) is unknown,

Authorized licensed use limited to: Univ of Calif San Francico. Downloaded on June 24,2010 at 04:31:39 UTC from IEEE Xplore. Restrictions apply.

KUMIHASHI AND SEKIHARA: ARRAY-GAIN CONSTRAINT MINIMUM-NORM SPATIAL FILTER

1361

we propose using the estimated source distribution s(r) forcomputing P (r).
The proposed recursive algorithm is described in the following manner.
1) By setting an initial P (r) to be the identity matrix I, theweight W (r) is derived using (24) and (25).
2) The estimated source distribution s(r, t) is obtained using(4).
3) This s(r, t) is then used to compute Ḡ at each time pointsuch thatL(r) 
(26)
s(r, t)
sT (r, t) LT (r) d3 rḠ =Ωwhereupon the weight matrix is updated using (24).
4) Finally, the updated source distribution is obtained using
the updated weight matrix with (4).
This procedure is repeated until a stopping criterion issatisfied.
One important problem in practical applications is the influence of measurement noise. In biomagnetic measurements,
the condition number of the gram matrix obtained in (26) isgenerally significantly large and the spatial filter in (24) is notrobust to the measurement noise. Since the norm of the filter
weight W (r) is the noise power gain of the spatial filter [1],one way to derive a spatial filter that is relatively tolerant to themeasurement noise is to add an additional constraint that makesW (r) not greater than a certain bound, i.e., we derive such aspatial filter using
W (r) = arg min F(W ), subject to W T (r)L(r) = Λ,
W (r)

and

W (r) ≤ WB

(27)

where the noise gain of the resultant spatial filter is not greaterthan the bound WB [1]. This formulation leads to the weight
expressed as−1−1

 L(r)]−1
 L(r)[L (r)G
W (r) = G
where
 =
G

T

(28)


L(r)P (r)LT (r) d3 r + γI = Ḡ + γI.

(29)

Ω

The scalar parameter γ works as the regularization parameter
determined based on the SNR in the measurements. We use
this regularized version of the proposed spatial filter in the experiments in the following sections. In these experiments, we
determine γ to be equal to the noise floor of the eigenvalue
 i.e., we use
spectrum of the measurement covariance matrix, R,γ that is equal to the mean of the noise-level eigenvalues of R.
Actually, considering the difference of the scaling between R
and Ḡ, we derive γ such that


Ḡ
1 
λi
(30)
γ=

Kn
R


noise

where noise λi indicates the summation over the noise eigenvalues, and Kn is total number of those eigenvalues.

Fig. 1. Coordinate system and source-sensor configuration used in the computer simulation. The coordinate origin was set at the center of the sensor coillocated at the center of the array. The plane at x = 0 cm is shown. The large
circle shows the cross-section of the sphere used for the forward calculation,and the small circles show the locations of the three sources. The square shows
the region for displaying the reconstruction results in Figs. 2– 7.

IV. COMPUTER SIMULATION
A. Reconstruction From Single-Time-Sample Data
A computer simulation was conducted to check the performance of the proposed AGMN-RUG spatial filter. In this computer simulation, we used a sensor alignment of the Magnes
2500 (4-D Neuroimaging, Inc., San Diego, CA) whole-head
sensor array, in which 148 sensors were arranged on a helmetshaped surface. Three sources were assumed to exist on thevertical single plane: (x = 0 cm); their (y, z) coordinates chosen as (−3.5, −6.5) cm, (1.0, −5.0) cm, and (4.0, −8.75) cm,
respectively. The source-sensor configuration and the coordinate system are depicted in Fig. 1. The spherical homogeneous
conductor model [9], with the sphere origin set at (0, 0, −12)cm, was used for the forward calculation. First, the sensor data
at a single time sample was calculated. Here, the intensities of
the three sources were set equal in the sensor domain. Gaussian
noise was added, and the SNR of the data was set to 16 where
the SNR was defined as the ratio of the norm of the signal vector
to the norm of the noise vector.
The source reconstruction was performed using sLORETA,
and the proposed AGMN-RUG spatial filter. Here, sLORETA
was chosen for comparison, because it is a method representative of existing nonadaptive methods and is able to perform 3-Dvolume reconstruction. The source space Ω was defined as the
3-D region between −4 ≤ x ≤ 4 cm, −5 ≤ y ≤ 5 cm, and−11 ≤ z ≤ −3 cm. The 3-D reconstruction was performed onΩ with a 1-cm voxel interval. The squared source intensity|s(r)|2 was plotted on the plane x = 0 cm on which the threesources were presumed to exist. The results from sLORETA
and the proposed spatial filter are shown in Fig. 2(a) and (b).
It can be clearly seen that the proposed spatial filter reconstructs the three sources at approximately correct locations. On

Authorized licensed use limited to: Univ of Calif San Francico. Downloaded on June 24,2010 at 04:31:39 UTC from IEEE Xplore. Restrictions apply.

1362

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 57, NO. 6, JUNE 2010

Fig. 2. Results of source reconstruction for single-time-sample data. (a) Results when the sLORETA spatial filter was applied. (b) Results when the proposed AGMN-RUG spatial filter was applied. The contours show the relativevalues of the squared source intensity |
s(r)|2 on the plane x = 0 cm on which
three sources exist.

Fig. 4. (Top to bottom) Time courses s1 (t), s2 (t), and s3 (t), used for generating the spatiotemporal dataset. The dataset was generated over 200 timesamples. The abscissa indicates the time samples and the ordinate indicates therelative amplitudes.

Fig. 3. Results of the source reconstruction obtained using the proposed
AGMN-RUG spatial filter with four different iteration numbers. (a) Numberof iterations was set to 4. (b) Number of iterations set to 16. (c) Number of iterations set to 256. (d) Number of iterations set to 1024. The contours show the
relative values of the squared source intensity |
s(r)|2 on the plane x = 0 cm.

the contrary, sLORETA produces blurred results, and fails toreconstruct the first and the second sources.
In this implementation of the proposed AGMN-RUG spatial
filter, we set the number of iterations to eight. Since this number was empirically determined, we checked how the proposed
method is sensitive to the number of iterations by using the same
single-time-sample data with four different iteration numbers.
The results obtained with iteration numbers of 4, 16, 256, and1024 are shown in Fig. 3(a), (b), (c), and (d), respectively. It
can be seen that these four sets of results are nearly identical,
suggesting that the method is rather stable with regard to thenumber of iterations. We present a brief discussion regarding
this stability of the proposed method in Section VI.
B. Reconstruction From Spatiotemporal Dataset
Next, we performed experiments on the reconstruction basedon a spatiotemporal dataset. The three time courses shown in
Fig. 4 were assigned to the three sources, and the spatiotemporal
data were generated over 200 time samples. These three timecourses are denoted s1 (t), s2 (t), and s3 (t). These time courses

Fig. 5. Results of source reconstruction for the spatiotemporal data. (a) Results
when the sLORETA spatial filter was applied. (b) Results when the proposed
AGMN-RUG spatial filter was applied. (c) Results when the minimum-variancespatial filter was applied. The contours show the relative values of the squared
source intensity |s(r, t)|2  on the plane x = 0 cm, where · indicates theaverage over the 200 time samples. (d) Reconstructed time courses of the threesources obtained using the proposed AGMN-RUG spatial filter. The three timecourses are the time courses at voxels that give the three maxima in (b).

are largely uncorrelated and the correlation coefficients between
any pair of the three are less than 0.1. The SNR was set equal
to 16 (where the SNR was defined as the ratio of the Frobenius
norm of the signal matrix to that of the noise matrix).
We applied sLORETA, the proposed spatial filter, and theminimum-variance spatial filter, to this spatiotemporal dataset.
The source-power distribution |s(r, t)|2  (where the notation
· indicates the average over the 200 time samples), was reconstructed over the source space Ω with a voxel interval of 1 cm,
and the results for the plane x = 0 cm are displayed in panels (a),(b), and (c) of Fig. 5. When implementing the proposed spatial
filter, the number of iterations was set to eight. The sample covariance matrix was computed using all 200 time samples when
implementing the minimum-variance filter reconstruction.

Authorized licensed use limited to: Univ of Calif San Francico. Downloaded on June 24,2010 at 04:31:39 UTC from IEEE Xplore. Restrictions apply.

KUMIHASHI AND SEKIHARA: ARRAY-GAIN CONSTRAINT MINIMUM-NORM SPATIAL FILTER

Fig. 6. Results of source reconstruction from the spatiotemporal data in which
the first and the second sources are highly correlated. (The correlation coefficient
was set to 0.98.) (a) sLORETA spatial filter was used. (b) Proposed AGMNRUG spatial filter was used. (c) Minimum-variance spatial filter was used. The
contours show the relative values of the squared source intensity |s(r, t)|2  on
the plane x = 0 cm, where · indicates the average over the 200 time samples.

It can be seen that sLORETA again produces rather blurred
results and fails to reconstruct the second source. The proposed
spatial filter and the minimum-variance spatial filter, however,are successfully able to reconstruct the three sources. The results of the proposed method have a spatial resolution intermediate to those of sLORETA and of the minimum-variance filter.
The reconstructed time courses of the three sources appear in
Fig. 5(d). These time courses are the reconstructed voxel timecourses obtained by the proposed AGMN-RUG spatial filter; thevoxels are those that give the three maxima seen in Fig. 5(b).
The results demonstrate that the method is capable of providing
considerably accurate time-course reconstruction.
We then performed experiments in which sources were highly
correlated. Without changing the time courses of the second andthe third sources, we assigned the time course s1 (t) to the firstsource, obtained using s1 (t) = (1 − α)s1 (t) + αs2 (t), where
α controls the degree of correlation between the first and thesecond sources. This α was set to 0.8, and the correlation coefficient between s1 (t) and s2 (t) turned out to be 0.98. The
results of the source reconstruction, in this case, are shown in
Fig. 6. The results of the minimum variance filter show a severe
influence from the source correlation known as signal cancellation [4], i.e., the intensities of the first and second sources are
significantly diminished. On the other hand, the results of the
proposed AGMN-RUG spatial filter show nearly no influenceof the source correlation, thereby demonstrating the robustnessof the proposed spatial filter to highly correlated sources.
We finally checked the performance of the proposed spatial
filter when SNR is considerably low. Here, other conditions
for data generation were the same as those for Fig. 5. The
reconstruction results from the proposed spatial filter for four1363

Fig. 7. Results of source reconstruction for four cases of SNR. Other datageneration conditions were the same as those for Fig. 5. (a) Results when SNR
was set at eight. (b) Results when SNR set at four. (c) Results when SNR set
at two. (d) Results when SNR set at one. The contours show the relative valuesof the squared source intensity |s(r, t)|2  on the plane x = 0 cm, where ·
indicates the average over the 200 time samples.

cases of SNR are shown in Fig. 7. Results when SNR was set ateight, four, two, and one are, respectively, presented in Fig. 7(a),(b), (c), and (d). These results show that the spatial resolution isaffected by the measurement SNR. When SNR is equal to one,the second and the third sources, which are approximately 5-cm
apart, are not resolved.
V. EXPERIMENTS USING AUDITORY-EVOKED MEG DATA
The effectiveness of the proposed spatial filter was further
tested using auditory-evoked MEG. Measurements were performed using the 160-channel MEG-VISION whole-cortex biomagnetometer (Yokogawa Electric, Inc., Tokyo, Japan). Here, a
1-kHz pure tone of 500 ms duration was delivered to the subject’s left ear, and the auditory-evoked field was measured at asampling rate of 10 kHz. The average interstimulus interval was
2 s. An online low-pass filter with a cutoff frequency of 2 kHz
was applied. A total of 400 trials were measured and averaged.
The averaged waveforms of the auditory-evoked field are shownin Fig. 8(a).
A single-time-sample data at the latency of 120 ms was used
for the reconstruction experiments. This time point, indicatedby a vertical broken line in Fig. 8(a), is close to the peak vertexof the M100, and the auditory cortices in the left and righthemispheres are expected to be simultaneously active at thisinstant. This is confirmed by the contour plots [see Fig. 8(b)] of
the sensor data at this time point, because the dipolar patterns
are visible above the left and right hemispheres. The resultsof the reconstruction experiments are shown in Fig. 9. In thisfigure, the reconstructed source distributions are overlaid onto
the sagittal, coronal, and axial slices of the subject MRI.
The results of the minimum-variance filter are shown in
Fig. 9(a). To obtain these results, the sample covariance matrix
was computed using a time window between 0 and 300 ms, andAuthorized licensed use limited to: Univ of Calif San Francico. Downloaded on June 24,2010 at 04:31:39 UTC from IEEE Xplore. Restrictions apply.

1364

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 57, NO. 6, JUNE 2010

the source power distribution (obtained by averaging the squared
source intensity over this time window) is shown in Fig. 9(a). In
these results, the minimum-variance spatial filter fails to reconstruct the sources of the bilateral activation of auditory cortices,
probably because the activities in these cortices are highly correlated. The results from sLORETA are shown in Fig. 9(b) andthe results from the proposed spatial filter are shown in Fig. 9(c).
Those in Fig. 9(b) and (c) were obtained using the single-timesample data at 120 ms. The results of sLORETA are significantly
blurred and biased. On the contrary, the proposed AGMN-RUG
spatial filter can reconstruct the bilateral activities of the auditory cortices with significantly high spatial resolution.
VI. DISCUSSION

Fig. 8. (a) Averaged waveforms of the auditory-evoked field used for thereconstruction experiments. The vertical broken line indicates the latency at
120 ms. (b) Sensor-field map at a latency of 120 ms. The contours show the
relative intensity of the measured magnetic field. The letters, A, P, L, and R,
respectively, show the anterior, posterior, left-hemisphere, and right-hemispheredirections.

We have mentioned that Eq. (26) is used to compute Ḡ in the
update process. However, instead, we can also useL(r) diag 
(31)
Ḡ =
s(r, t)
sT (r, t) LT (r) d3 rΩwhere diag[A] indicates the diagonal matrix whose diagonal
elements are equal to those of a matrix A. In this case, Ḡapproaches to the model covariance matrix for a case in which ateach source location there are three independent sources whose
powers are equal to sx (r, t)2 , sy (r, t)2 , and sz (r, t)2 . We
found that (31) gives, in most cases, spatial resolution higher
than that from (26), although the difference is generally small.
Therefore, we have used (31) in our experiments in Sections IV
and V.
To derive the proposed spatial filter, we used the optimization
formulation in (22). A slight modification of this formulation iscapable of providing another type of recursively implemented
spatial filter, which is derived using
W (r) = arg min F(W ),subject to

W (r)





T
W T (r)L(r  ) P (r  ) W T (r)L(r  ) d3 rΩ= P (r)(32)

where the cost function F(W ) is defined as
F(W ) =
tr{ W T (r)L(r  ) − δ(r − r  )IΩ
T
P (r  ) W T (r)L(r  ) − δ(r − r  )I } d3 r .

(33)

By setting P (r) = |s(r, t)| I, the resultant spatial filter is expressed as
2

−1

−1

W (r) = |s(r, t)|Ḡ L(r) LT (r)Ḡ L(r)Fig. 9. Results of the source reconstruction for the auditory-evoked field
shown in Fig. 8. The reconstruction results are superimposed onto the subject’s MRI, and the axial, sagittal, and coronal slices are shown in the left,
middle, and right panels. These slices were chosen as the slices that contained the maximum source intensity. (a) Results obtained from the minimumvariance spatial filter. The source power distribution obtained by averaging thesquared source intensity over the time window between 0 and 300 ms is shown.
(b) Results obtained from the sLORETA spatial filter. (c) Results obtained fromthe proposed AGMN-RUG spatial filter. The two sets of results in (b) and(c) were obtained using the single-time-sample data at 120 ms.

−1/2

.

(34)

The spatial filter in Eq. (34) is very similar to the method known
as the sLORETA-FOCUSS algorithm [10], [11].1
1 Strictly speaking, the sLORETA-FOCUSS algorithm described in [10] and[11] is different from the recursive spatial filter given in (34), because thealgorithm in [10] and [11] uses the sLORETA method only for the initial
iteration, and uses the weighted minimum-norm method for other iterations.
In contrast, in (34), the sLORETA method is used for all iterations with theweighting of |s(r, t)| obtained in the preceding iteration.

Authorized licensed use limited to: Univ of Calif San Francico. Downloaded on June 24,2010 at 04:31:39 UTC from IEEE Xplore. Restrictions apply.

KUMIHASHI AND SEKIHARA: ARRAY-GAIN CONSTRAINT MINIMUM-NORM SPATIAL FILTER

The difference between the spatial filter derived above andthe proposed AGMN-RUG filter is that the previous spatial
filter explicitly contains the estimated source distribution s(r, t)
in the weight equation. Therefore, the method generally gives
a so-called sparse solution. On the contrary, in the proposed
AGMN-RUG spatial filter, the filter weight is updated in an
indirect manner, only through the updated gram matrix, andthe sparsity is not imposed on the solution. It may probably betrue that the proposed spatial filter is superior in the stability,
for example, the stability to the number of iterations as wasdemonstrated in our computer simulation. On the other hand, themethod expressed in (34) can attain the spatial resolution higher
than that of the proposed method. A performance comparison
between these two methods are currently in progress, and results
will be presented in a near future occasion.
It should be pointed out that the weights of spatial filtersderived with the array-gain constraint are unitless and the reconstructed source vector has the unit of the magnetic field, i.e.,
the intensity of the reconstructed source is equal to the dipole
moment multiplied by the norm of the lead field at the source
location. Therefore, the source dipole moment is retrieved bymultiplying the reconstructed source intensity with the inverseof the lead-field norm at each source location.
In summary, we have proposed a novel AGMN-RUG spatial
filter, which is developed based on a modified version of theminimum-norm method. The proposed spatial filter uses the estimated theoretical covariance matrix, and is designed to havea performance close to that of the adaptive minimum-variancespatial filter. Since the proposed method does not use a sample covariance matrix, a large number of time samples are notneeded. The method is also robust to the source correlation.
This robustness is demonstrated in our computer simulation andexperiments using auditory-evoked MEG data.
The proposed method, however, cannot generally replace theminimum-variance spatial filter because the proposed method
attains a spatial resolution only intermediate to those ofsLORETA and of the minimum-variance filter, as was shown in
our experiments. Therefore, in situations where the minimumvariance spatial filter can attain its full performance, there should
be no advantage of using the proposed method. We are now
in the process of conducting a series of experiments which
will evaluate the method’s performance in various conditions.
The results of those experiments will be published in a futureoccasion.
ACKNOWLEDGMENT
The authors would like to thank H. Tanaka for his help inproviding auditory-evoked MEG data. They are also grateful to
Dr. D. Palomo for his careful proofreading.
REFERENCES[1] K. Sekihara and S. S. Nagarajan, Adaptive Spatial Filters for Electromagnetic Brain Imaging. Berlin, Germany: Springer-Verlag, 2008.
[2] K. Sekihara, M. Sahani, and S. S. Nagarajan, “Location bias and spatial
resolution of adaptive and non-adaptive spatial filters for MEG sourcereconstruction,” NeuroImage, vol. 25, pp. 1056–1067, 2005.

1365[3] 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, no. 9,
pp. 867–880, Sep. 1997.
[4] K. Sekihara, S. S. Nagarajan, D. Poeppel, and A. Marantz, “Performanceof an MEG adaptive-beamformer technique in the presence of correlated
neural activities: Effects on signal intensity and time-course estimates,”IEEE Trans. Biomed. Eng., vol. 49, no. 12, pp. 1534–1546, Dec. 2002.
[5] S. S. Dalal, K. Sekihara, and S. S. Nagarajan, “Modified beamformers
for coherent ref_pubtitle region suppression,” IEEE Trans. Biomed. Eng.,
vol. 53, no. 7, pp. 1357–1363, Jul. 2006.
[6] M. S. Hämäläinen and R. J. Ilmoniemi, “Interpreting measured magnetic
fields of the brain: Estimates of current distributions,” Helsinki Univ.
Technol., Espoo, Finland, Tech. Rep. TKK-F-A559, 1984.
[7] R. E. Greenblatt, A. Ossadtchi, and M. E. Pflieger, “Local linear estimators
for the bioelectromagnetic inverse problem,” IEEE Trans. Signal Process.,
vol. 53, no. 9, pp. 3403–3412, Sep. 2005.
[8] R. D. Pascual-Marqui, “Standardized low resolution brain electromagnetic
tomography (sLORETA): Technical details,” Methods Findings Exp. Clin.
Pharmacol., vol. 24, pp. 5–12, 2002.
[9] J. Sarvas, “Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem,” Phys. Med. Biol., vol. 32, pp. 11–22, 1987.
[10] R. Khemakhem, B. Hamida, A. Ahmed-Taleb, and P. Derambure, “Newhybrid method for the 3d reconstruction of neuronal activity in the brain,”in Proc. 15th Int. Conf. Syst., Signals Image Process., 2008, pp. 405–408.
[11] K. Rafik, B. H. Ahmed, F. Imed, and T.-A. Abdelmalik, “Recursive
sLORETA-FOCUSS algorithm for EEG dipoles localization,” in Proc.
Conf. Image Process. Theory, Tools Appl., Nov. 2008, pp. 1–5.

Isamu Kumihashi received the B.E. degree from
Tokyo Metropolitan Institute of Technology, Tokyo,
Japan, in 2007, and the M.E. degree from Tokyo
Metropolitan University, Tokyo, in 2009.
He is currently with the Department of Systems
Design and Engineering, Tokyo Metropolitan University. His research interests include the neuromagnetic
source imaging and the reconstruction methods, especially their application to imaging of the cardiac
electrical activity.

Kensuke Sekihara (M’88–SM’06–F’09) received
the M.S. degree in 1976 and the Ph.D. degree in 1987
both from Tokyo Institute of Technology, Tokyo,
Japan.
From 1976 to 2000, he was with the Central Research Laboratory, Hitachi Ltd., Tokyo. From 1985
to 1986, he was a Visiting Research Scientist with
Stanford University, Stanford, CA. From 1991 to
1992, he was with the Basic Development, Siemens
Medical Engineering, Erlangen, Germany. From1996 to 2000, he was engaged in the “Mind Articulation” research project sponsored by Japan Science and Technology Corporation. He is currently a Professor with the Department of Systems Design andEngineering, Tokyo Metropolitan University, Tokyo. His research interests include the electromagnetic Source imaging, reconstruction and inverse methods,
multidimensional signal processing, and statistical signal processing, especially
their application to functional Source imaging. He is author of a book Adaptive
Spatial Filters for Electromagnetic Brain Imaging (Springer-Verlag, 2008).
Dr. Sekihara is a member of the IEEE Medicine and Biology Society and theIEEE Signal Processing Society. He is a Fellow of the International Society ofFunctional Source Imaging (ISFSI).

Authorized licensed use limited to: Univ of Calif San Francico. Downloaded on June 24,2010 at 04:31:39 UTC from IEEE Xplore. Restrictions apply.