92

IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 18, NO. 1, JANUARY 1999

Rapid Communications
Time–Frequency MEG-MUSIC Algorithm
Kensuke Sekihara,* Srikantan Nagarajan,
David Poeppel, and Yasushi Miyashita

Abstract—We propose a method that incorporates the time–frequency
characteristics of neural sources into magnetoencephalographic (MEG)
source estimation. The method is based on the multiple-signalclassification (MUSIC) algorithm and it calculates a time–frequency
matrix in which diagonal and off-diagonal terms are the auto andcrosstime–frequency distributions of multichannel MEG recordings,
respectively. The method averages this time–frequency matrix over
the time–frequency region of interest. The locations of neural sources
are then estimated by checking the orthogonality between the noise
subspace of this averaged matrix and the sensor lead field. Accordingly,
the method allows us to estimate the locations of neural sources fromeach time–frequency component. A computer simulation was performed
to test the validity of the proposed method, and the results demonstrate
its effectiveness.
Index Terms— Biomagnetism, biomedical signal processing, inverse
problems, time–frequency analysis.

Fig. 1. Schematic view of the sensor-source configuration used in thecomputer simulation.

indicates the matrix transpose and I indicates the unit matrix. The
eigenvalues are numbered in decreasing order.

I. INTRODUCTIONT

HE noninvasive measurement of magnetic fields generated
from human cortical neural activities, referred to as magnetoencephalography (MEG) [1], has been found to be a powerfultool in studies of human neurophysiology and neural information
processing. One major problem here is the MEG inverse problem [1],a problem of estimating neural current distributions from the magnetic
field measured outside a human head. Because the neural current
distribution is inherently three-dimensional (3-D), the estimation
problem is generally ill posed. To reduce this ill posedness, theestimation needs to incorporate some prior knowledge regarding
the source characteristics. Such characteristics can include possible
source locations, the source spatial extent, the total number of sources,or the source-frequency characteristics. In this paper, we propose anovel method that incorporates source time–frequency characteristics
into the source estimation. The method combines time–frequency
analysis [2] with the MEG multiple-signal-classification (MUSIC)
algorithm [3] and allows estimation of neural sources from eachtime–frequency component. Throughout this paper, plain italics
indicate scalars, lower case boldface italics indicate vectors, andupper case boldface italics indicate matrices. The superscript T
Manuscript received July 15, 1998; revised December 7, 1998. This work
was supported in part by the Centre National de la Recherche Scientifique.
The Associate Editor responsible for coordinating the review of this paper
and recommending its publication was M. Viergever. Asterisk indicatescorresponding author.
*K. Sekihara is with the Mind Articulation Project, Japan Science andTechnology Corporation (JST), Yushima 4-9-2, Bunkyo, Tokyo 113-0034
Japan (e-mail: sekihara@ma.jst.go.jp).
S. Nagarajan is with the Keck Center for Integrative Neuroscience, University of California, San Francisco, CA 94143-0732 USA.
D. Poeppel is with the Department of Linguistics and Biology, Universityof Maryland, College Park, MD 20742 USA.
Y. Miyashita is with the Mind Articulation Project, Japan Science andTechnology Corporation (JST), Tokyo 113-0034 Japan and is also with the
Department of Physiology, The University of Tokyo, School of Medicine,Tokyo 113-0033 Japan.
Publisher Item Identifier S 0278-0062(99)02029-7.

II. METHOD
Let us define the magnetic field measured by the mth detector coilat time t as bm (t) and a 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. A total of P current-dipole sources is assumed to generate thebiomagnetic field, and the locations of these sources are denoted as(x 1 ; x2 ; . . . ; x P ): The magnitude of the pth dipole-source moment isdefined as sp (t): The source magnitude vector is defined as s(t) =T(s1 (t); s2 (t); . . . ; sP (t)) : A spherical homogeneous conductor [4]is assumed and two tangential components, the  and  components,of the dipole-source moment are considered. The dipole orientation
is defined as its normal vector  p (t) = (p (t); p (t))T wherek p (t)k = 1: We also define a 2P 2 P matrix that expresses theorientations of all P dipole sources as 9 (t) such that

9 (t) =

 1 (t)

0

111

0

0

 2 (t)

1

..
.

..
.
0

1
111

..

.

0

:

(1)

0
 P (t)

The lead field vectors for the  and  components of the sourceT
at x are defined as l  (x) = (l1 (x); l2 (x); . . . ; lM(x)) and l (x) =



T


(l1 (x); l2 (x); . . . ; lM (x)) : Here, lm (x) and lm (x) express the mth
sensor output induced by the unit-magnitude source moment directedin the  and  directions, respectively. We define the lead field matrix
for the source at x as L(x) = [l  (x); l  (x)], which represents thesensitivity of the sensor array at location x: The composite leadfield matrix for the entire set of P dipole sources is defined as
Lc = [L(x1 ); L(x2 ); . . . ; L(xP )]: Then, the relationship betweenb(t) and s(t) is expressed asb(t) = [Lc9 (t)]s(t) + n(t)where n(t) is the additive noise.

0278–0062/99$10.00  1999 IEEE

(2)IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 18, NO. 1, JANUARY 1999

93

(a)

(b)

(c)
Fig. 2. Generated waveforms of  components of: (a) first source; (b) second source; and (c) third source.

To formulate the time–frequency MUSIC algorithm, we assume
that all dipole sources do not change their orientations during the time
period of interest and we use time-independent 9 to express the dipole
source orientations. The time–frequency distribution matrices for b(t)and s(t) are defined as C b (t; f ) and C s (t; f ), which can be obtainedby using any of the Cohen-class time–frequency distributions [2],[5], i.e.,

C b (t; f ) =

1

01

1
and

C s (t; f ) =

1
01

1

8( 0 t;  0 f )
101b( + =2)bT ( 0 =2)e02i d d d
(3)

TABLE I
SOURCE PARAMETER VALUES ASSUMED FOR THE COMPUTER SIMULATION

its off-diagonal elements are the crosstime–frequency distributions
between different channel recordings. The matrix C s (t; f ) is a
P 2 P matrix. Its diagonal elements are the autotime–frequency
distributions of the source activities and its off-diagonal elements
are the crosstime–frequency distributions between different source
activities.
Using (2)–(4), we can derive

C b (t; f ) = (Lc9 )C s (t; f )(9 T LcT )+ C sn (t; f ) + C ns (t; f ) + C n (t; f )8( 0 t;  0 f )
101s( + =2)s ( 0 =2)e02i d d d:T(4)Here, 8(t; f ) is the convolution kernel, which determines the characteristics of the resultant Cohen-class time–frequency distributions.
The matrix C b (t; f ) is an M 2 M matrix. Its diagonal elements are
the autotime–frequency distributions of the channel recordings andwhere

C sn (t; f ) =

1
01

1

(5)

8( 0 t;  0 f )1(Lc9 )s( + =2)011 nT ( 0 =2)e02i d d d

(6)

94

IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 18, NO. 1, JANUARY 1999

(a)

(b)
Fig. 3. (a) Waveform of the generated magnetic field from one representative channel located above the left hemisphere. (b) Power spectrum of thegenerated magnetic field obtained by averaging the spectra from all channels.

C ns (t; f ) =

1
01

1

8( 0 t;  0 f )
101n( +  =2)

1 sT ( 0  =2)(9 T LcT )e02i d d d (7)andC n (t; f ) =

1
01

1

8( 0 t;  0 f )
101n( +  =2)nT ( 0  =2)e2i d d d:(8)

Let us define the region of interest in the time–frequency domain as: This region contains the signal of interest and we localize sources
by using time–frequency components in this region. The matrices 7 b
and 7 s are obtained by averaging C b (t; f ) and C s (t; f ) over thistarget time–frequency region

7b =, i.e.,C b (t; f ) dt df and 7 s =C s (t; f ) dt df (9)where s s indicates the integral over this target time–frequency
region. We denote the number of sources whose activities havetime–frequency components in the target region as P : Note thatthe matrix 7 s is a P 2 P matrix.

Fig. 4. Smoothed pseudo-Wigner–Ville time–frequency distribution of computer generated neuromagnetic data. The time–frequency maps were calculated from all channel recordings, and the results obtained by averaging allthese maps are shown. The frequency is normalized such that the Nyquist
frequency is equal to 0.5. The three quadrangular areas denoted by 1 ; 2 ;and 3 are the target regions used for the source-estimation experiments.

We assume that the noise and signal are uncorrelated in thetarget region. This assumption leads to s s C ns (t; f ) dt df =s s C sn (t; f ) dt df = 0: We also assume that the noise is white

IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 18, NO. 1, JANUARY 1999

95

(a)

(b)

(c)

(d)

Fig. 5. MUSIC-localizing function on the plane x = 1 cm. The time–frequency MUSIC algorithm [(13)] was used with the target region set at (a)(b) 2 ; and (c) 3 : (d) The conventional MUSIC algorithm [(17)] was used.

Gaussian, so the expected value of C n (t; f ); hC n (t; f )i is equal to2 I where h1i indicates the ensemble average and 2 is the noise
DDpower density. Accordingly, if the target region is sufficiently large,
replacing the ensemble average by the integration over the target
region we have s s C n (t; f ) dt df =  2 I where  2 is the totalnoise power in the target region. Therefore, we finally derive therelationship7 b = (Lc 9 )7 s (9 T LcT ) + 2I :(10)

Denoting the noise-level eigenvectors of 7 b as uj ; where j+ 1; . . . ; M , (10) becomesP(7 b 0 2I )uj = (Lc9 )7 s (9 T LcT )uj = 0;
j = P + 1; . . . ; M:

=

(11)

The matrix Lc is a full-column-rank matrix. Thus, when the matrix
7 s is a full-rank matrix, we get(9T) = 0;

for j = P

LcT uj+ 1; . . . ; M:

(12)

The above equation indicates that source locations can be obtainedby checking the orthogonality between the sensor lead field with
an optimum orientation L(x) opt where  opt represents the normal
vector in the optimum orientation, and the noise-level eigenvectors

uj (j = P

+ 1; . . . ; M ):This orthogonality can be evaluated by calculating the following
function:

J (x) = 1=min (L ( )T()

x Z N Z TN L x ; LT(x)L(x))(13)1;

where min (1; 1) indicates the minimum generalized eigenvalue ofthe matrix pair in parenthesis and the matrix Z N is defined as
Z N = [uP +1 ; . . . ; uM ]: This J (x), referred to as the localizing
function, is calculated point-by-point in a volume where sources can
exist. Each location where J (x) forms a peak is chosen as the locationof one dipole source.
Equation (13) was obtained under the assumption that dipole
sources do not change their orientations during the time period ofinterest. It is worth noting that when this assumption does not holdand the source orientations change, this localizing function can still
be used to detect source locations. In such cases, instead of (10) wehave the relationship7 b = Lc7~ sLcT + 2I

where 7~ s = s s

C~ s (t; f ) =

1

01

1

(14)

C~ s (t; f ) dt df; and C~ s (t; f ) is obtained from8( 0 t;  0 f )101

9 ( +  =2)s( +  =2)

1 sT ( 0  =2)9 T ( 0  =2)e02i d d d:(15)
Therefore, assuming that the matrix 7~ s is a full-rank matrix, we also
obtain the relationship

LcT uj = 0;

for j = P

+ 1; . . . ; M:

(16)

96

IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 18, NO. 1, JANUARY 1999

(a)

(b)

(c)
Fig. 6. Smoothed pseudo-Wigner–Ville time–frequency distribution of the neuromagnetic data used for the time jitter experiments. The time–frequency maps
were calculated from all channel recordings, and the results obtained by averaging all these maps are shown. (a) Time–frequency map obtained from the
waveform averaged across all epochs. (b) Time–frequency map obtained by averaging maps from all raw epochs. (c) Time–frequency MUSIC localizing
function on the plane x = 1 cm. The target region was set at 1 shown in Fig. 4 and the matrix 7 b was calculated using 7 b = hs s C b (t; f ) dt df iepoch :Since the source locations that satisfy (16) also satisfy (12), the localizing function in (13) is effective even when the source orientations
are not time independent and change during the period of interest.
It is informative to compare (13) with the localizing function usedin the conventional time-domain MUSIC algorithm. This conventional localizing function is expressed as [3]T L(x); LT (x)L(x))
J (x ) = 1=min (LT (x)E N E N(17)where E N = [eP +1 ; . . . ; eM ] and ej (j = P + 1; . . . ; M ) indicatesthe noise-level eigenvectors of the measured-data covariance matrix.
This measured-data covariance matrix is obtained by averaging
b(t)bT (t) during the time period of interest. The only difference
between the proposed and the conventional MUSIC algorithms is thatthe proposed algorithm uses Z N instead of E N in the conventionalalgorithm.
III. COMPUTER SIMULATION
We performed a computer simulation to test the validity of the
proposed method. We used the coil configuration of a 148-channel
whole-head Magnes 2500WHTM biomagnetic measurement system(Biomagnetic Technologies Inc., San Diego, CA). The coordinate
origin was defined at the center of the coil array. The z direction was
defined as the direction perpendicular to the plane of the detector coillocated at this center. The x direction was defined as that from theposterior to the anterior, and the y direction was defined as that fromthe left to the right hemisphere.
Three signal dipole sources were assumed to exist on the sameplane (x = 1:0 cm). The configuration of this computer simulation isshown schematically in Fig. 1. The simulated magnetic field was calculated at 351 time points. To generate the simulated neuromagnetic
field, the  components of the three sources wj (t); (j = 1; 2; 3)
were all Gaussian amplitude modulated and frequency modulated,i.e., wj (t) = exp[0(t 0 tj )2 =(2j2 )] cos[2 (aj t + bj )(t 0 tj )]: The
 components of the three sources were set to zero. The values for
tj ; j ; aj ; bj and the assumed source locations are listed in Table I,
and the generated waveforms for the  components of the threesources are shown in Fig. 2.
Uncorrelated Gaussian noise was added to make the final signalto-noise ratio (SNR) equal to two. The SNR was defined as the ratio

IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. 18, NO. 1, JANUARY 1999

of the Frobenius norm of the signal-magnetic-field data matrix to thatof the noise matrix. The waveform of the generated magnetic fieldfrom one representative channel located above the left hemisphere isshown in Fig. 3(a) and the power spectrum obtained by averaging
the spectra from all channels is shown in Fig. 3(b). These figuresindicate that the signals from the three sources cannot be separated
in either the time domain or in the frequency domain.
The smoothed pseudo-Wigner–Ville transform was calculated inorder to obtain the time–frequency distribution of the computer
generated data. The convolution kernel in (3) and (4) was chosensuch that it was equal to a multiplication of a function of t andthat of f each derived by Fourier transforming the hamming window
with a window size of 39 data points. The time–frequency maps were
calculated from all channel recordings, and the results obtained byaveraging all these maps are shown in Fig. 4. Three time–frequency
components are clearly resolved.
The time–frequency MUSIC algorithm was then applied by setting
the target time–frequency region as the regions indicated by 1 ; 2 ;and 3 in Fig. 4. The results are shown in Fig. 5. The contours
in Fig. 5(a)–(c) show the relative value of the MUSIC localizing
function in (13) on the plane x = 1 cm and each area where thelocalizing function reaches a peak is considered to be the locationof one dipole source. The (y; z ) coordinates of the peak locations in
Fig. 5(a), (b), and (c) are (3.9, 06.5), (4.8, 08.1), and (5.4, 09.5),
respectively. Comparison with the original source locations listed inTable I shows that the time–frequency MUSIC algorithm accurately
localized the three sources, and these results verify the validity ofthe proposed algorithm. The conventional MUSIC algorithm [3] wasapplied to the same computer-generated data [Fig. 5(d)], but theconventional localizing function in (17) was not able to resolve thethree sources. The results in Fig. 5 demonstrate the effectiveness ofthe proposed method.
It is sometimes possible that a transient event-related signal is nottime locked to the stimulus and has a time jitter among multipleepoch recordings [6]. We next show that a minor modification of the
proposed algorithm allows source locations to be estimated from sucha nontime-locked event-related signal. The simulated magnetic fielddata was calculated assuming that the first source has a time jitter tz ,i.e., its  component has a form w1 (t 0 tz ): Here, tz was generated byusing a uniform random number ranging from -50 to 50 time points,
and a total of 64 raw-epoch data was generated. The SNR of eachraw-epoch data was set at 0.25. Two kinds of time–frequency maps
were obtained; one was obtained by using the simulated magneticfield waveforms averaged across all raw epochs and the other was
obtained by averaging time–frequency maps across all raw epochs.
The results are shown in Fig. 6.
In the time–frequency distribution in Fig. 6(a), the time–frequency
component generated from the first source has almost disappeared
because of the time jitter. Conversely, this component is clearlyobservable (although it is significantly blurred compared to thatin Fig. 4) in the time–frequency distribution in Fig. 6(b), which
is obtained by averaging a time–frequency map from each set of
raw-epoch data. This observation suggests that if we calculate thematrix 7 b such that 7 b = hs s C b (t; f ) dt df iepoch where
h1iepoch indicates the average over epochs we can localize sources
for nontime-locked activities. The localization function in (13) was
calculated with 7 b obtained in the above-mentioned manner and withthe same target region 1 shown in Fig. 4. The results are shownin Fig. 6(c). Here, the first source is clearly detected. The (x; y )coordinates of the peak in this figure are (4.1, 06.5), which are very
close to the assumed location of the first source. The results obtained97

here show that the proposed algorithm is, in principle, effective for
localizing an evoked source whose activation is not time-locked tothe stimulus.
IV. DISCUSSION
Although the quadratic Cohen-class distributions can provide highresolution time–frequency maps, they are known to contain significant
amount of spurious crossterms when the distribution contains multiple
components. The convolution kernel 8(t; f ) in (3) or (4) should
therefore be chosen so as to reduce the amount of such spurious
cross terms at the minimum sacrifice of the resolution. As a result,the kernel determines the tradeoff between the resolution and theamount of such spurious crossterms. The relationship between several
specific choices of the kernel and the properties of the resultant
time–frequency distributions has been well studied in the field oftime–frequency analysis. A detailed discussion on the choices of thekernel and properties of the resultant time–frequency distributions
are found in [2] and [5].
In our computer simulation, we used the smoothed pseudoWigner–Ville representation, which uses the kernel separable interms of time and frequency. This choice, however, is somewhat
arbitrary and any other Cohen-class representation can be used forour computer generated data. For actual spontaneous or event-related
MEG, the choice of the kernel must be made by taking several
properties of the MEG data to be analyzed into consideration. Suchproperties include the number of time–frequency components, howclose they are, and how localized or diffused they are. The optimum
kernel for analyzing various types of transient MEG data needs to
be investigated.
The proposed algorithm is suited to localizing sources for transientand nonstationary MEG signals. There are several situations in
MEG measurements in which a transient signal is observed but its
origin is not adequately established. A well-known example of suchsignals is visual or auditory gamma-band activity, which may have arelationship with feature binding in a cognitive process or even withconsciousness [6], [7]. Other examples of such signals include humanspontaneous MEG [8]. We plan to apply the proposed algorithm tosuch nonstationary signals and we will report these results in thenear future.
REFERENCES[1] M. Hämäläinen, R. Hari, R. J. IImoniemi, J. Knuutila, and O. V.
Lounasmaa, “Magnetoencephalography-theory, instrumentation, and applications to noninvasive studies of the working human brain,” Rev.
Mod. Phys., vol. 65, pp. 413–497, 1993.
[2] L. Cohen, Time–Frequency Analysis: Theory and Applications. Englewood Cliffs, NJ: Prentice-Hall, Inc., 1994.
[3] 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.
[4] J. Sarvas, “Basic mathematical and electromagnetic concepts of thebiomagnetic inverse problem,” Phys. Med. Biol., vol. 32, pp. 11–22,
1987.
[5] F. Hlawatsch and G. F. Boudreaux-bartels, “Linear and quadratic
time–frequency signal representations,” IEEE Signal Processing Mag.,
vol. 9, pp. 21–67, Apr. 1992.
[6] 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.
[7] R. Llinas and U. Ribary, “Coherent 40-Hz oscillation characterizes
dream state in humans,” Proc. Nat. Acad. Sci., vol. 90, pp. 2078–2081,
1993.
[8] R. Hari and R. Salmelin, “Human cortical oscillations: a neuromagnetic
view through the skull,” Trends Neurosci., vol. 20, pp. 44–49, 1997.