642

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 47, NO. 5, MAY 2000

Estimating Neural Sources from
Each Time-Frequency Component ofMagnetoencephalographic Data
Kensuke Sekihara*, Member, IEEE, Srikantan S. Nagarajan, David Poeppel, Satoru Miyauchi, Norio Fujimaki,
Hideaki Koizumi, and Yasushi Miyashita

Abstract—We have developed a method that incorporates
the time-frequency characteristics of neural sources into magnetoencephalographic (MEG) source estimation. This method,referred to as the time-frequency multiple-signal-classification
algorithm, allows the locations of neural sources to be estimatedfrom any time-frequency region of interest. In this paper, weformulate the method based on the most general form of thequadratic time-frequency representations. We then apply it to
two kinds of nonstationary MEG data: gamma-band (frequency
range between 30–100 Hz) auditory activity data and spontaneous
MEG data. Our method successfully detected the gamma-bandsource slightly medial to the N1m source location. The method
was able to selectively localize sources for alpha-rhythm bursts
at different locations. It also detected the mu-rhythm sourcefrom the alpha-rhythm-dominant MEG data that was measured
with the subject’s eyes closed. The results of these applications
validate the effectiveness of the time-frequency MUSIC algorithm
for selectively localizing sources having different time-frequency
signatures.
Index Terms—Biomagnetism, biomedical signal processing,
functional brain imaging, inverse problems, time-frequency
analysis.

I. INTRODUCTIONTHE NONINVASIVE measurement of magnetic fields generated from human cortical neural activities, referred to as
magnetoencephalography (MEG) [1], has been found to be a
powerful tool in studies of human neurophysiology and neural
information processing. One major problem with MEG measurements is the inverse problem [1], a problem of estimating

Manuscript received December 10, 1998; revised November 29, 1999. This
work was carried out as part of the MIT-JST International Cooperative Research
Project “Mind Articulation.” Asterisk indicates corresponding author.
*K. Sekihara is with the Mind Articulation Project, Japan Science and Technology Corporation (JST) Yushima 4-9-2, Bunkyo-ku, Tokyo, 113-0034, Japan(e-mail: sekihara@ma.jst.go.jp).
S. S. Nagarajan is with the Keck Center for Integrative Neuroscience, University of California, San Francisco, CA 94143-0732 USA. He is also with the
Scientific Learning Corporation, Berkeley, CA 94704 USA.
D. Poeppel is with the Department of Linguistics and Biology, University of
Maryland, College Park, MD 20742 USA.
S. Miyauchi and N. Fujimaki are with the Kansai Advanced Research Center,
Communications Research Laboratory, Hyogo, 651-2401, Japan.
H. Koizumi is with the Mind Articulation Project, Japan Science and Technology Corporation (JST) Bunkyo-ku, Tokyo, 113-0034, Japan.
Y. Miyashita is with the Mind Articulation Project, Japan Science and Technology Corporation (JST) Bunkyo-ku, Tokyo, 113-0034, Japan. He is also with
the Department of Physiology, The University of Tokyo, School of Medicine,
Bunkyo-ku, Tokyo, 113-0033, Japan.
Publisher Item Identifier S 0018-9294(00)03282-1.

neural current distributions from the magnetic field measured
outside a human head. Because this estimation problem is inherently ill-posed, the estimation needs to incorporate some priorknowledge regarding the source characteristics. Such characteristics can include possible source locations, the source spatial
extent, the total number of sources, or the frequency characteristics of source activities.
Neuronal activities in a human brain are basically nonstationary. Time-frequency analysis has, therefore, played an important role in analyzing and characterizing electrophysiological data such as that from MEG or encephalography (EEG)[2]–[4]. We have developed a method that incorporates the timefrequency characteristics of the source activity into the source
estimation. The method, referred to as the time-frequency multiple-signal classification (MUSIC) algorithm, allows us to estimate the locations of neural sources from any time-frequency
region of interest (ROI). We have reported, in a short communication, the formulation of a method based on the Cohen-classquadratic time-frequency representations [5]. The present paper
presents a formulation based on the most general form of thequadratic time-frequency representations. It also presents the results of applying the developed method to two types of nonstationary MEG data: stimulus-locked auditory gamma-band dataand spontaneous MEG (alpha and mu rhythm) data.
A recent study investigated the possibility of incorporating
source estimation with time-frequency analysis in the form of atwo step procedure [6]. They first synthesize a signal that only
contains the target time-frequency components of interest, andthen apply a single-dipole localization procedure to this synthesized signal. This method relies on a technique referred to as
time-frequency synthesis. This technique, however, is tediousand may be computationally intensive when an ROI with anarbitrary shape is used. In contrast, our method, which estimates neural sources directly from the time-frequency domain
data, is more straightforward. Note that a method similar to our
method was independently developed and reported recently inthe radar-signal processing field [7].
The purpose of this paper is to describe the formulationof the time-frequency MUSIC algorithm, and to present thedetails of our investigation that validate this formulation.
Following the method’s formulation, which is based on general
quadratic time-frequency representations, we present theresults from numerical experiments that not only validated themethod’s formulation, but also showed the basic superiorityof the method over existing source-localization methods. We

0018–9294/00$10.00 © 2000 IEEE

SEKIHARA et al.: ESTIMATING NEURAL SOURCES FROM EACH TIME-FREQUENCY COMPONENT OF MEG DATA

then apply the time-frequency MUSIC algorithm to auditory
gamma-band evoked-field data and spontaneous MEG (alphaand mu rhythm) data. These data were used because the their
origins have been relatively well investigated, and we canassess the validity of the time-frequency MUSIC results bycomparing them with results from previous investigations[8]–[12]. The results of these applications validated the basic
effectiveness of the time-frequency MUSIC algorithm whenanalyzing such nonstationary MEG signals as mixtures ofvarious time-frequency components.
Throughout this paper, a lower-case boldface font indicates avector, and an upper-case boldface font indicates a matrix. The
superscript indicates the matrix transpose and the superscript
indicates the Hermitian conjugate transpose of a matrix.
II. METHOD
A. Definitions
Let us define the magnetic field measured by theth detector coil at timeasand a vector
as a set of measured
is the total number of detector coils. A total
data where
current-dipole sources are assumed to generate the bioof
magnetic field, and the locations of these sources are denoted. The magnitude of the th dipole-sourceas. The source magnitude vector is
moment is defined as. The sphericaldefined ashomogeneous conductor model [13] is used, so we consider
two tangential components, the and components, of thesource moment while its radial component is assumed to be
zero. The dipole orientation is defined as its normal vector, where. We also define amatrix that expresses the orientations of all dipole
such thatsources as
..
.
..
.

..

(1)

.

In this paper, for simplicity, we assume that no dipole sourceschange their orientation during the time period of interest, andwe use time-independent to express the dipole-source orientations.
andcompoThe lead-field vectors for theare defined asnents of the source moment at
and. Here,andexpress
the th sensor output induced by the unit-magnitude source
directed in the and directions, respectively. We define thelead-field matrix for the source at aswhich represents the sensitivity of the sensor array at location. The composite lead-field matrix for the entire set of
dipole sources is defined as.
andis expressed as
Then, the relationship between
(2)
where

is the additive noise.

643

B. Conventional Time-Domain MUSIC Algorithm
The MEG inverse problem is the problem of estimating
the source locations
and the source activity
from the measurement. The conventional way to solve
this inverse problem is to minimize the least-squares error
whereis the estimated measurement.
Estimating the source locations, in this minimization, requires adimensional nonlinear search. Generally, for such a highly
multidimensional search, there is no guarantee of obtaining
a correct solution unless we can set the initial estimate very
close to the true solution. The MUSIC algorithm approachhas been introduced [14], [15] in order to avoid this highly
multidimensional search. A distinct advantage of this algorithm
is that regardless of the number of dipole sources, it gives a
suboptimal estimate of the source locations by using only athree-dimensional search in the solution space.
We define the measured-data covariance matrix as , and the. Using
covariance matrix of the dipole-source activities as(2) and assuming that the noise and signal are uncorrelated, we
have

(3)
indicates the ensemble average, and it is assumed thatwhere
the noise is white Gaussian noise with variance . This equation shows the relationship between the covariance matrices
of measurement and source activities; this relationship is thebasis of deriving the MUSIC algorithm. We denote the eigenas, where(the eigenvectors ofvalues are numbered in decreasing order), and define matrixas. The eigenvectors
are those corresponding to the noise-level eigenvalues, and theis called the noise subspan of the columns in the matrix
space.
To estimate the locations of the sources, the MUSIC algorithm takes advantage of the fact that the sensor lead-field with, in which
represents thean optimum orientation
normal vector in the optimum orientation, is orthogonal to theat the true source locations. This orthognoise subspace of
onality can be evaluated by the following function [15], whichis called the MUSIC localizer(4)
indicates the generalized minimum eigenvalue
whereof the matrix pair given in parenthesis. Equation (4) is calculated
in a volume where sources can exist, and each location wherereaches a peak is chosen as the location of one source.
C. Time-Frequency Domain MUSIC Algorithm
1) Covariance-Matrix Relationship in the Time-Frequency
Domain: Here, we describe the time-frequency domain
MUSIC algorithm, which makes it possible to estimate neuralsources from each time-frequency component of MEG data.
Time-frequency signal representations are generally classified

644IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 47, NO. 5, MAY 2000

into two groups: linear and quadratic representations [16].
Although the linear representations include the well-known
short-time Fourier transform (STFT) and the wavelet transform,
quadratic time-frequency representations are known to provide
better performance, and thus many modern time-frequency
representations have quadratic forms [17]. We thus use ageneral quadratic time-frequency representation to formulate
the method, but it can also be easily formulated using linear
time-frequency representations.
Let us denote the time-frequency representation matrix of theasand that ofas. In thevector signal
most general quadratic form, they are expressed as

and

(10)
Let us introduce the target time-frequency region which contains the signal of interest whose source configuration is to be investigated. We assume that the noise and signal are uncorrelated
in the target region, and that the noise is white Gaussian and isuncorrelated between channels. Then, we obtain
and
where
again indiis the noise-power density incates the ensemble average andthe time-frequency domain. We also assume that the noise has
the ergodic property, so the ensemble average can be approximated with the average over the target region, i.e.,

(5)
and

and
(6)
is an
matrix. Its diagonal eleThe matrix
ments are the auto-time-frequency distributions of the channel
recordings and its off-diagonal elements are the cross-time-frequency distributions between different channel recordings. The
is amatrix. Its diagonal elements are thematrix
auto-time-frequency distributions of the source activities andits off-diagonal elements are the cross-time-frequency distributions between different source activities.
is a kernel, which deIn (5) and (6),termines the characteristics of the resultant quadratic
time-frequency distributions [17]. When this kernel has, the resultant
the formtime-frequency representations are called the Cohen class.
When the kernel has a form of the affine smoothing, i.e.,, whereisa constant that converts scale to frequency, the quadratic
time-frequency representations are called the affine class. The
relationship between several specific choices of the kernel andthe properties of the resultant time-frequency distributions has
been well studied in the field of time-frequency analysis. A
detailed discussions on such studies are found in [16]–[18].
Using (2), (5), and (6), we can derive

We also define the matrices
andobtained by averaging
andover the target time-frequency region

and
(11)
Then, we finally derive the relationship(12)This equation shows the relationship between the measurement
and source-activity covariance matrices in the time-frequency
domain [its conventional time-domain counterpart is shown in(3)], and is the basis of deriving the time-frequency MUSIC
algorithm.
2) Deriving the Time-Frequency MUSIC Localizer: We
asdenote the noise-level eigenvectors of, whereis the number of sources whose
activity has time-frequency components in the target region.
Then, since the noise-level eigenvectors satisfy the relationship,, (12) leads to

(7)
(13)

where
and
The matrices
fore, when the matrix
(8)

(9)

are full-column-rank matrices. Thereis a full-rank matrix, we get
for

(14)

The above equation indicates that source locations can be obtained by checking the orthogonality between the optimally-ori-

SEKIHARA et al.: ESTIMATING NEURAL SOURCES FROM EACH TIME-FREQUENCY COMPONENT OF MEG DATA

Fig. 1. (a) Waveform of the generated magnetic field from one representative
channel located approximately above the sources. (b) Power spectrum of
the waveform shown in (a). (c) Power spectrum of the generated magnetic
field obtained by averaging the spectra from all channels. The frequency isnormalized such that the Nyquist frequency is equal to 0.5.

645

Fig. 3. Results of calculating the MUSIC localizer on the plane x = 1 cm. The
time-frequency MUSIC algorithm was used with the target region set as: (a)and (b). (c) The conventional MUSIC algorithm was used. The contours
show the relative values of the MUSIC localizer.

Fig. 4. The x, y , and z coordinates used to express the estimated results in
Figs. 6, 8, 12, and 14. The midpoint between the left and right preauricular points
is defined as the coordinate origin. The axis directed away from the origin toward
the left preauricular point is defined as the +y axis and that from the origin to
the nasion is the +x axis. The +z axis is defined as the axis perpendicular toboth these axes and is directed from the origin to the vertex.
Fig. 2. Spectrogram obtained by averaging the spectrograms from all channel
recordings. The two quadrangular areas denoted byand, are the target
regions used for the source-estimation experiments.

ented sensor lead field
and the noise-level eigenvectors.
Therefore, the localizing function can be calculated in thefollowing manner:(15)

where the matrix
is defined as.
The only difference between the proposed and the conventionalMUSIC algorithms is that the proposed algorithm uses thisinstead of theused in (4), and that since the matrixis. Equation
complex-valued, we use the Hermitian conjugate(15) was obtained under the assumption that dipole sources do
not change their orientations during the time period of interest.
It can easily be shown that this localizer in (15) is effective evenwhen this assumption does not hold and the source orientations
change during the measurement [5].

646IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 47, NO. 5, MAY 2000

(a)Fig. 5. MEG waveforms from two kinds of auditory measurements. (a)Regular auditory N1m measurement. (b) Auditory gamma-band measurement.
The recordings from 60 sensors covering the subject’s right hemisphere are
overlapped.

3) Summary of the Algorithm: The procedure of implementing the time-frequency MUSIC algorithm is summarized
as follows.
is calculated
Step 1) The time-frequency matrix
with an appropriate choice of the time-frequency
representations. The th diagonal term of this matrix is given by the auto-time-frequency distributionthof the th-channel recording, and its
off-diagonal term is given by the cross-time-frequency
th andth channels.
distribution between theis obtained by averaging
Step 2) The matrix
over the predetermined target region , and the noiselevel eigenvectors are obtained by applying the eigen.
decomposition of
Step 3) The localizer defined by (15) is calculated throughout
reachesthe source space, and each location wherea peak is determined to be the location of one dipole
source.
III. NUMERICAL EXPERIMENTS
We performed numerical experiments to test the validity of
the arguments in Section II. We used the coil configuration of a148-channel whole-head Magnes 2500 WH biomagnetic measurement system (Biomagnetic Technologies Inc., San Diego,
CA). The coordinate origin was defined at the center of the coilarray. The direction was defined as the direction perpendicular to the plane of the detector coil located at this center. The
direction was defined as that from the posterior to the anterior,and the direction was defined as that from the left to the righthemisphere. Two signal sources were assumed to exist on thecm) with theircoordinates equal tosame plane (
and
.

(b)

(c)
Fig. 6. Results of applying the conventional time-domain MUSIC algorithm tothe data shown in Fig. 5(a). The time window ranging from 0–130 ms was used
to calculate R . The contours show the relative values of the localizer in (4);the projections of the localizer values onto the transverse, coronal, and sagittal
planes are shown. The circles depicting a human head represent the projectionsof the sphere used to calculate the forward solution.

The simulated magnetic field was calculated at
351 time points. To generate the simulated neurocomponents of the two sourcesmagnetic field, thewere amplitude- and frequency-modulated; i.e.,. The values forwere set at,
for the first,forsource andthe second source where the unit of is timepoint. Their
components were set to zero.
Gaussian noise uncorrelated between different sensor recordings was added to make the final signal-to-noise ratio (SNR)
equal to 0.85. The SNR was defined as the ratio of the Frobenius norm of the signal-magnetic-field data matrix to that of
the noise matrix. The waveform of the generated magnetic fieldfrom one representative channel located approximately above

SEKIHARA et al.: ESTIMATING NEURAL SOURCES FROM EACH TIME-FREQUENCY COMPONENT OF MEG DATA

647

(a)

Fig. 7. Spectrogram for the frequency range above 15 Hz obtained from
the MEG data in Fig. 5(b). The spectrogram was calculated using the first
singular temporal vector, which was obtained by applying the singular value
decomposition to the 60-channel spatio–temporal matrix of the averaged data.
The windowing function was set equal to a 255-point Hanning window.

the two sources is shown in Fig. 1(a). Its power spectrum isshown in Fig. 1(b), and the power spectrum averaged across all
sensor channels is shown in Fig. 1(c). The spectrograms were
calculated from all channel recordings, and the results obtainedby averaging all these spectrograms are shown in Fig. 2. Here, aHanning window with a size of 77 data points was used [19].
Figs. 1 and 2 indicate that the signals from the two sources
cannot be separated in either the time domain or in the frequency
domain, but they are clearly separated in the time-frequency domain.
The time-frequency MUSIC algorithm was then applied by
setting the target time-frequency region as the regions indicatedandin Fig. 2. The results are shown in Fig. 3. The conby
tours in Fig. 3(a) and (b) shows the relative value of the MUSICcm, and each
localizing function in (15) on the plane
area where the localizing function reaches a peak is considered) coordinates
to be the location of one dipole source. The (andof the peak locations in Fig. 3(a) and (b) are, respectively. These results show that the time-frequency MUSIC algorithm accurately localized these two signal
sources, and these results verify the validity of our argument
in Section II. The conventional time-domain MUSIC algorithm
was also applied to the same computer-generated data, but theconventional algorithm was not able to resolve the two sources[Fig. 3(c)].

IV. EXPERIMENTS USING MEG DATAThe MEG data presented in this section was recorded byusing a 148-channel whole-head biomagnetometer (Magnes
2500WH) that was installed at the Communications ResearchLaboratory, Tokyo, Japan. We use the head coordinates explained in Fig. 4 to express the estimated source locations.

(b)

(c)
Fig. 8. The results of calculating time-frequency MUSIC localizer [(15)] withthe target region set at the gamma-band region marked as the square in Fig. 7.

A. Source Localization for Auditory Gamma-Band Activity
Two kinds of auditory measurements were conducted with ahealthy male volunteer (age 39, right-handed). First, we measured the regular auditory magnetic response by using a 1-kHz
pure-tone stimulus applied to the subject’s left ear. The sampling interval in the data acquisition was 1 kHz, the averageinter-stimulus interval was 2 s, and an on-line bandpass filter
ranging from 1–400 Hz was applied. A total of 100 epochs wereaveraged, and this averaged data is shown in Fig. 5(a). Here weselected recordings from 60 sensors covering the subject’s righthemisphere and the waveforms for those recordings are overlapped in this figure. The peak for the N100m (the peak locatedaround the latency of 100 ms) is clearly visible. Second, we
measured the gamma-band auditory activity using a click soundwith a duration of 1 ms. The click sound was applied to the subject’s left ear; the inter-stimulus interval was 1 s and a total of
400 epochs were averaged. In this measurement, to maintain thesubject’s attention to the stimulus, the subject was asked to movehis right index finger when he heard the tone. The waveforms
from the same 60 sensors are shown in Fig. 5(b).

648IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 47, NO. 5, MAY 2000

(a)

(b)

Fig. 9. MRI overlay of the gamma-band source and the N1m source both obtained from three separate measurement sessions. The white and the black triangles
indicate the locations of the gamma-band and the N1m sources, respectively. The coronal (left) and the sagittal (right) cross sections are shown.

Fig. 10. Spectrograms calculated from the temporal singular vectors of the 148-channel recordings of the spontaneous MEG. (a) The first singular vector is used.
(b) The second singular vector is used. The windowing function was set equal to a 127-point Hanning window. The regions denoted as
,
, and
were
chosen as the target regions in the time-frequency MUSIC analysis.

We applied the conventional time-domain MUSIC algorithm
to the auditory N100m shown in Fig. 5(a). The covariance mawas calculated using this averaged response with a timetrix
window ranging from 0–130 ms. The MUSIC localizer plot isshown in Fig. 6. In this figure, the MUSIC localizing functionin (4) was calculated with an interval of 0.5 cm over a volume
in the subject’s right hemisphere, and the projections onto thetransverse, the coronal, and the sagittal planes are displayed.
The figure shows that the algorithm detected the source for the
N100m response, which is believed to be located in the primary
auditory cortex area in the right hemisphere.
To extract a representative time-course of the 60-channel
recordings shown in Fig. 5(b), we applied the singular value
decomposition to the spatio–temporal data matrix of these
averaged recordings. The spectrogram was then calculatedusing the first temporal singular vector, which corresponds
to the largest singular value. The windowing function wasset equal to a 255-point Hanning window. The results areshown in Fig. 7. The stimulus-locked gamma-band activity

is clearly observed in the latency range from 0–0.1 s. The
MUSIC localizer plot obtained using (15) is shown in Fig. 8.
Here, the procedure described in Section II-C-3 was applied
to the averaged data in Fig. 5(b). That is, the averaged datain (5). The square in Fig. 7 indicates thewas used astarget region; i.e., only the gamma-band portion within thissquare was used to obtain the plot in Fig. 8. In this figure,the time-frequency MUSIC method clearly detected a single
localized source for the auditory gamma-band activity in theright temporal area.
With this subject, we repeated the same measurements threetimes; each measurement session was separated by approximately one week. The three sets of localized results for the
gamma-band activity were overlayed onto the subject’s magnetic resonance image (MRI). The results are shown in Fig. 9.
The three locations of the gamma-band source were all medial
to the N100m source but differed somewhat: One was locatednearly 1 cm medial to the N100m source, but the other two wereclose to the N100m source.

SEKIHARA et al.: ESTIMATING NEURAL SOURCES FROM EACH TIME-FREQUENCY COMPONENT OF MEG DATA

649

B. Source Localization for Spontaneous MEG
The subject was a healthy male volunteer (age 27,
right-handed) who was asked to lie down quietly with hiseyes closed during the measurements. The spontaneous MEG
was recorded continuously for a total of six seconds with asampling frequency of 677~Hz. An on-line bandpass filter
ranging from 1~Hz to 200~Hz was applied. Singular-value
decomposition of the 148-channel spatio-temporal data matrix
was performed, and spectrograms were obtained from the first
and the second temporal singular vectors. The windowing
function was set equal to a 127-point Hanning window. The
resultant spectrograms are shown in Fig. 10. A part of the
148-channel recordings of the spontaneous MEG is shown in
Fig. 11. A couple of the alpha spindles can be observed in thisfigure.
Two distinct alpha bursts are clearly observed in Fig. 10(a),which was obtained from the first singular vector. The regionsandas shown
including these components are denoted asin Fig. 10(a). Then, the time-frequency MUSIC algorithm wasand . The MUSIC
applied with the target region set at
localizer plots are shown in Fig. 12(a) and (b). Fairly localized
sources were detected and the locations of the sources fromandwere found in the occipital region. The region marked
in Fig. 10(b) was also analyzed by the time-frequency
byMUSIC algorithm. The resultant MUSIC localizer plot is shownin Fig. 12(c) A single source is clearly detected in the left temporal area.
The peak locations in Fig. 12(a)– (c) were overlayed onto thesubject’s MRI; the results for sagittal cross sections are shown
andin Fig. 13. The two sources for the alpha bursts inare, respectively, localized near the parieto-occipital sulcus andin the occipital visual area [Fig. 13(a) and (b)]. Several previous
investigations reported that the source for the alpha rhythm was
localized near the parieto-occipital sulcus and near the calcarine
sulcus in the visual area [10]–[12]. It was also reported that eachalpha burst (spindle) has its own location [20], [21]. Our results
here are in good agreement with the results from these previous
investigations. Fig. 13(c) shows that the source for a component
was localized in the motor area, and this indicates thatin
this component is a part of the mu-rhythm [10], [12]. These
results show that the time-frequency MUSIC algorithm is so
sensitive that it was able to detect the mu-rhythm source fromthe alpha-rhythm-dominant spontaneous MEG.
V. DISCUSSION
In the results presented in Section IV-A, the three locationsof the gamma-band source were all medial to the N100m
source but differed somewhat. A previous investigation [8],[9] reported that the source of the gamma-band activity was
localized slightly medial and anterior to the N100m source.
Our results are in fairly good agreement with these previous
results. Several factors could have caused the fluctuation in
the estimated source locations of the auditory gamma-band
activity. Such factors include errors in registering the head
coordinates to the MRI pixel coordinates, the subject’s headmotion during the measurement, and the noise contained in theFig. 11. The measured spontaneous MEG. Twenty-six representative
channel recordings were selected. (The selection was made with no particularpreference, and these 26 channel sensors are fairly uniformly located on a headsurface.) This spontaneous MEG was measured for 6 s with a subject lying
down quietly with his eyes closed.

target time-frequency region. Even fluctuations in the subject’s
internal conditions, such as his attention to the stimulus [9] or
his general vigilance level, may have caused fluctuation of theestimated source locations. Further investigation using multiple
subjects is needed to draw clearer conclusions about the locationof the source of the gamma-band activity.
We also applied the conventional time-domain MUSIC algorithm to the spontaneous MEG data in Section IV-B for comparison. The results are shown in Fig. 14. Here, data from thewhole 6 s was used to obtain . The results indicate that theconventional MUSIC algorithm can resolve none of the sources
detected by the time-frequency MUSIC algorithm, and the conventional MUSIC algorithm only detected a large activated areain the occipital region. This comparison demonstrates the effectiveness of the proposed time-frequency MUSIC algorithm forselectively localizing sources having different time-frequency
signatures.
In Section II, the time-frequency MUSIC algorithm was derived from the matrix equation [(12)] with the assumption that
there are no perfectly correlated sources. The MUSIC algorithm has a well-known problem in that it tends to mis-localize
sources when they are highly correlated. It is, therefore, worth
noting that the developed method can be modified to incorporate non-MUSIC-type multidipole estimation. It can be shownthat the optimum estimates of the source locations can be obtained by directly minimizing the least-squares-error cost function given by
(16)
where

is the projection operator defined as(17)

and are the estimated values ofand . This minand
imization can be done using conventional nonlinear optimiza-

650

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 47, NO. 5, MAY 2000

(a)

(b)

(c)

Fig. 12.

The results of applying the time-frequency MUSIC algorithm to the spontaneous MEG data. The target region was set at (a)tion methods [19]. Thus, if it is suspected that correlated multiple sources exist in the time-frequency ROI, we can apply thisdirect minimization procedure although we need information
about the number of sources and a good initial guess concerning
the source locations.
The most computationally intensive part of implementing
the time-frequency MUSIC algorithm is the calculation of, which contains
times calculations of the cross-, (b), and (c).

or auto-time-frequency distribution of channel recordings. In
Section IV, we used the STFT because of its computational
efficiency. For the case presented in Section IV-A, it tookapproximately one hour1 to calculate, but we estimated

1This calculation was done by using Pentium III (450 MHz)-based PC with
MATLAB (Mathworks, Natick, MA).

SEKIHARA et al.: ESTIMATING NEURAL SOURCES FROM EACH TIME-FREQUENCY COMPONENT OF MEG DATA

Fig. 13.

651

MRI overlay of the peak locations in Fig. 12. The results in (a), (b), and (c), respectively, show the peak locations obtained in Fig. 12(a)–(c).

that it would have taken 3–4 h to calculate this matrix if we had
used the smoothed pseudo Wigner distribution, which is a typical Cohen-class quadratic representation. The current computerroutine, however, is not optimized to obtain the highest performance in the computational speed, and it is possible to attain aconsiderable speedup, for example, by programing it using theC code. For analyzing MEG data such as those presented in Section IV, any time-frequency representations could be used withthe time-frequency MUSIC algorithm because high time-frequency resolution is not needed to analyze such MEG data.
However, when analyzing other kinds of MEG data, such as
some types of sleep spindles or epileptic spikes, it may be necessary to use more sophisticated time-frequency representations,since analyzing such data may require higher time-frequency
with a resresolution. Computer routine that calculates
onable computational speed should be developed in such cases.
A time-frequency representation of multichannel MEG
recordings provides space, time, and frequency information,and many time-frequency maps whose number equals thenumber of sensor channels are obtained. Such information
may be excessive for extracting the time-frequency characteristics of the signal of interest, and we need to obtain
fewer time-frequency maps that appropriately represent thewhole time-frequency properties of the multichannel data. To
extract such representative time-frequency maps, we applied
singular-value decomposition to the measured spatio–temporal

data matrix in Section IV, and calculated the time-frequency
maps of the temporal singular vectors corresponding to significantly large singular values. Other method can also be used,
for example, one that calculates the time-frequency distribution
averaged over the channels [this is equivalent to calculating thetrace of the matrix]. A simpler method such as thismay suffice in some cases such as our numerical experiment in
Section III.
In summary, we have developed a method that incorporates
the time-frequency characteristics of neural sources into magnetoencephalographic (MEG) source estimation. The method,which we refer to as the time-frequency MUSIC algorithm, allows the locations of neural sources to be estimated from any
time-frequency ROI. The algorithm was applied to two typesof nonstationary MEG data: auditory gamma-band activity dataand spontaneous MEG data. The effectiveness of the time-frequency MUSIC algorithm was assessed by comparing the results of these applications with the results from previous investigations on the same types of nonstationary MEG data. As a next
step, we plan to apply the proposed time-frequency MUSIC algorithm to MEG data whose origin is difficult to analyze withexisting methods. Possible candidates include some types of
sleep spindles, epileptic spikes, and spontaneous waves withpathologic origins. The results of these investigations will bereported in the near future.

652

IEEE TRANSACTIONS ON BIOMEDICAL ENGINEERING, VOL. 47, NO. 5, MAY 2000

(a)

(b)

(c)
Fig. 14. The results of applying the conventional time-domain MUSIC
algorithm (4) to the spontaneous MEG data. Data from the whole 6 s was used
to calculate the covariance matrix R .

[4] E. Rodriguez, N. George, J. P. Lachaux, J. Martinerie, B. Renault, andF. J. Varela, “Perception’s shadow: Long-distance synchronization ofhuman brain activity,” Nature, vol. 397, pp. 430–433, 1999.
[5] K. Sekihara, S. Nagarajan, D. Poeppel, and Y. Miyashita, “Time-frequency MEG-MUSIC algorithm,” IEEE Trans. Med. Imag., vol. 18, pp.
92–97, 1999.
[6] M. Sun, S. Qian, X. Yan, S. B. Baumann, X.-G. Xia, R. E. Dahl, N.
D. Ryan, and R. J. Sclabassi, “Localizing functional activity in the brain
through time-frequency analysis and synthesis of the EEG,” Proc. IEEE,
vol. 84, pp. 1302–1311, 1996.
[7] M. Amin and A. Belouchrai, “Time-frequency MUSIC,” in Proc. SPIE
Symp. Radar Processing, Technology, and Applications, San Diego, CA,
July 1998.
[8] C. Pantev, T. Elbert, S. Makeig, S. Hampson, C. Eulitz, and M. Hoke,“Relationship of transient and steady-state auditory evoked fields,” Electroenceph. Clin. Neurophysiol., vol. 88, pp. 389–396, 1993.
[9] C. Pantev and T. Elbert et al., “The transient auditory evokedgamma-band field,” in Oscillatory Event-Related Brain Dynamics, C.
Pntev et al., Eds. New York: Plenum, 1994, pp. 219–230.
[10] R. Hari and R. Salmelin, “Human cortical oscillations: A neuromagnetic
view through the skull,” Trends Neurosci., vol. 20, pp. 44–49, 1997.
[11] C. Tesche and M. Kajola, “A comparison of the localization of spontaneous neuromagnetic activity in the frequency and time domains,” Electroenceph. Clin. Neurophysiol., vol. 87, pp. 408–416, 1993.
[12] R. Salmelin and R. Hari, “Caracterization of spontaneous MEG rhythms
in healthy adults,” Electroenceph. Clin. Neurophysiol., vol. 91, pp.
237–248, 1994.
[13] J. Sarvas, “Basic mathematical and electromagnetic concepts of the biomagnetic inverse problem,” Phys. Med. Biol., vol. 32, pp. 11–22, 1987.
[14] R. O. Schmidt, “A signal subspace approach to multiple emitter locationand spectral estimation,” Ph.D. dissertation, Stanford University, Stanford, CA, 1981.
[15] J. C. Mosher, P. S. Lewis, and R. M. Leahy, “Multiple dipole modelingand localization from spatio–temporal MEG data,” IEEE Trans. Biomed.
Eng., vol. 39, pp. 541–557, 1992.
[16] F. Hlawatsch and G. F. Boudreaux-bartels, “Linear and quadratic timefrequency signal representations,” IEEE Signal Processing Magazine,
vol. 9, pp. 21–67, 1992.
[17] L. Cohen, Time-Frequency Analysis: Theory and Applications. Englewood Cliffs, NJ: Prentice-Hall, 1994.
[18] L. Cohen, “Time-frequency distributions—A review,” Proc. IEEE, vol.
77, pp. 941–981, 1989.
[19] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,
Numerical Recipes in Fotran: The Art of Scientific Computing, 2nded. New York: Cambridge Univ. Press, 1992.
[20] S. J. Williamson, J.-Z. Wang, and R. J. Ilmoniemi et al., “Method for
localizing sources of human alpha activity,” in Advances in Biomagnetism, S. J. Williamson et al., Eds. New York: Plenum Press, 1989,
pp. 257–260.
[21] R. M. Chapman et al., “EEG/MEG and alpha activity: Current statusof biomagnetic research,” in Biomagnetism Clinical Aspects, M. Hokeet al., Eds. Ansterdam, The Netherlands: Elsevier Science, 1992, pp.
25–31.

ACKNOWLEDGMENT
The authors would like to thank Dr. R. Takino for his help
in performing the MEG measurements. A part of the computer
simulation reported in this paper used the Time-Frequency
Toolbox provided by France’s CNRS (Centre National de laRecherche Scientifique).
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] C. Tallon-Baudry, O. Bertr, C. Delpuech, and J. Pernier, “Stimulus specificity of phase-locked and nonphase-locked 40-Hz visual responses inhuman,” J. Neurosci., vol. 16, pp. 4240–4249, 1996.
[3] C. Tallon-Baudry, O. Bertrand, C. Delpuech, and J. Pernier, “Oscillatory-band (30–70 Hz) activity induced by a visual search task in humans,”
J. Neurosci., vol. 17, pp. 722–733, 1997.

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

SEKIHARA et al.: ESTIMATING NEURAL SOURCES FROM EACH TIME-FREQUENCY COMPONENT OF MEG DATA

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–1998, he was a Postdoctoral Fellow
in the Keck Center for Integrative Neuroscience at
the University of California, San Francisco (UCSF).
Since 1999, he has been a full-time Research Scientist at Scientific Learning Corporation, Berkeley, CA,
and is an Adjunct Assistant Professor in the Department of Otolaryngology at
UCSF. In January 2000, he will start as an Assistant Professor in the Department of Bioengineering at the University of Utah, Salt Lake City. His researchinterests in neural engineering include bio-electromagnetism, systems and computational neuroscience and statistical signal processing.

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

Satoru Miyauchi received the B.S. degree in psychology from Waseda University, Tokyo, Japan, in
1977, and the Ph.D. degree in medical sciences at
Toho University, Tokyo, Japan, in 1990.
He is currently a Section Chief of Auditory
and Visual Informatics Section, Kansai Advanced
Research Center, Communications Research Laboratory, Hyogo, Japan. His research uses functional
magnetic resonance imaging (fMRI) and magnetoencephalography (MEG) to investigate the neuralbasis of visual attention and visual imagery.

653

Norio Fujimaki received the B.S., M.S., and Ph.D.
degrees in electronics engineering from the University of Tokyo, Tokyo, Japan, in 1975, 1977, and 1980,
respectively.
He worked in Fujitsu Laboratories Ltd., Atsugi,
Japan, from 1980 to 1999. In 1999, he joined Communications Research Laboratory. His research interests are superconducting devices and brain-function measurements using electroencephalograpy, magentoencephalography, and functional magnetic resonance imaging.

Hideaki Koizumi was born in Tokyo, Japan, in1946. He graduated from the Department of Pure andApplied Sciences at the University of Tokyo, Tokyo,
Japan in 1971. He received the doctoral degree inphysics from the University of Tokyo in 1976 for his
dissertation creating polarized Zeeman-effect atomic
absorption spectrometry.
He joined Hitachi, Ltd., Tokyo, Japan. After
receiving the Ph.D. degree, he remained as Visiting
Faculty at the Lawrence Berkeley Laboratory,
University of California, Livermore, from 1978 to
1979. Upon returning to Japan, he led MRI development projects. He now
leads projects on environmental measurement and analysis and noninvasive
higher-order brain function analysis. He is a Chief Research Scientist at theCentral Research Laboratory, Hitachi, Ltd. and a Professor of the ResearchInstitute for Electronic Science, Hokkaido University, Sapporo, Japan.
Dr. Koizumi received many awards such as the Ohkochi Memorial Prize, thePrize of Science and Technology Minister, and IR-100. He is also a member ofthe board of directors of the Chemical Society of Japan.

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