Journal of Neural Engineering
J. Neural Eng. 14 (2017) 051001 (20pp)

https://doi.org/10.1088/1741-2552/aa7693

Topical Review

Subspace-based interference removal
methods for a multichannel biomagnetic
sensor array
Kensuke Sekihara1,2

and Srikantan S Nagarajan3

1

Signal Analysis Inc., Hachioji, Tokyo, Japan
Department of Advanced Technology in Medicine, Tokyo Medical and Dental University, 1-5-45
Yushima, Bunkyo-ku, Tokyo 113-8519, Japan3
Biomagnetic Imaging Laboratory, University of California, San Francisco, 513 Parnassus Avenue,
S362, San Francisco, CA 94143, United States of America
2E-mail: k-sekihara@nifty.com
Received 27 January 2017, revised 4 May 2017
Accepted for publication 18 May 2017
Published 18 August 2017
Abstract

Objective. In biomagnetic signal processing, the theory of the signal subspace has been
applied to removing interfering magnetic fields, and a representative algorithm is the signal
space projection algorithm, in which the signal/interference subspace is defined in the spatial
domain as the span of signal/interference-source lead field vectors. This paper extends thenotion of this conventional (spatial domain) signal subspace by introducing a new definitionof signal subspace in the time domain. Approach. It defines the time-domain signal subspaceas the span of row vectors that contain the source time course values. This definition leads
to symmetric relationships between the time-domain and the conventional (spatial-domain)
signal subspaces. As a review, this article shows that the notion of the time-domain signal
subspace provides useful insights over existing interference removal methods from a unified
perspective. Main results and significance. Using the time-domain signal subspace, it ispossible to interpret a number of interference removal methods as the time domain signal
space projection. Such methods include adaptive noise canceling, sensor noise suppression,
the common temporal subspace projection, the spatio-temporal signal space separation, andthe recently-proposed dual signal subspace projection. Our analysis using the notion of thetime domain signal space projection reveals implicit assumptions these methods rely on, andshows that the difference between these methods results only from the manner of deriving
the interference subspace. Numerical examples that illustrate the results of our arguments areprovided.
Keywords: interference removal, magnetoencephalography, sensor array processing,
signal subspace, biomagnetic imaging, biomagnetism, multi-sensor array(Some figures may appear in colour only in the online journal)

1741-2552/17/051001+20$33.00

1© 2017 IOP Publishing Ltd

Printed in the UK

Topical Review

J. Neural Eng. 14 (2017) 051001

1. Introduction

time domain can be derived by those algorithms. We provide
numerical examples that illustrate the results of our arguments
in section 6. The Appendix provides several proofs of mathematical arguments, as well as some details of the SSS algorithm needed for readers to follow our arguments.

The notion of signal and noise subspaces has been considered
useful in the signal processing of data acquired using multichannel sensor arrays [1]. In biomagnetic signal processing,
the theory of the signal subspace has been applied to removing interfering magnetic fields, and a representative algorithm
in this regard is the signal space projection (SSP) algorithm[2–4]. In this algorithm, the interference subspace is defined as
the span of interference-source lead field vectors [2, 5], and it
is estimated by utilizing so-called empty-room-noise data. The
algorithm projects the measured data onto the subspace orthogonal to the interference subspace to remove the interference.
In biomagnetic signal processing, the signal subspace has
been defined to reflect a spatial property of the multichannel
data [2, 5]. On the other hand, since biomagnetic signals haverich temporal properties, the incorporation of temporal information can be a natural and fruitful way to extend the notionof the conventional (spatial domain) signal subspace. In thispaper, we propose, for the first time (as far as we know), a
new definition of signal subspace in the time domain. The
time-domain signal subspace is defined as the span of row
vectors that contain the source time course values. Such row
vectors are referred to as the source time course vectors.
By defining the time-domain signal subspace in thismanner, we derive symmetric relationships between thetime-domain signal subspace and the conventional (spatialdomain) signal subspace. That is, while a column vector ofsignal components at a particular time point lies within thespatial-domain signal subspace, a row vector of signal timecourse from a particular sensor lies within the time-domain
signal subspace. While the column space of the signal matrix
is equal to the spatial domain signal subspace, the row spaceof this matrix is equal to the time-domain signal subspace.
As a review article, this paper does not propose a new
method, but rather a new way of looking at various existing
methods in a unified perspective. Actually, using the timedomain signal subspace, it is possible to interpret various
interference removal methods as the time domain signal space
projection. It can be shown that these methods differ only inthe manner by which they derive the interference subspace.
These methods rely on some implicit assumptions that are
generally hidden behind their formulations, but they can berevealed by our analysis using the notion of the time domain
SSP. Such methods include adaptive noise canceling [6, 7],sensor noise suppression [8], common temporal subspace projection [9], spatio-temporal signal space separation [10] andthe recently-proposed dual signal subspace projection [11].
The paper is organized as follows: After reviewing the conventional (spatial-domain) signal subspace, the time-domain
signal subspace is defined in section 2. The interferenceremoval methods utilizing the spatial-domain signal subspace
are reviewed in section 3. An important variant of the SSP
algorithm, the signal space separation (SSS) method [12–14],
is also reviewed in this section. The time-domain SSP algorithm is introduced in section 4. Section 5 presents interpretations of various existing algorithms as the time-domain
SSP. Here, we argue how the interference subspace in the

2. Signal subspaces in the spatial and time
domains
2.1. Sensor array measurements

Biomagnetic measurement is usually conducted using a sensor array, which simultaneously measures the biomagnetic
signal with multiple sensors. Let us define the measurement of the mth sensor at time t as ym(t). The measurementfrom the whole sensor array is expressed as a column vectory(t): y(t) = [y1 (t), y2 (t), . . . , yM (t)]T . Here, M is the numberof sensors, and the superscript T indicates the matrix transpose.
Throughout this paper, plain italics indicate scalars, lower-caseboldface italics indicate vectors, and upper-case boldface italics
indicate matrices. The location in the three-dimensional spaceis represented by r: r = (x, y, z). The source magnitude at r andtime t is denoted as a scalar s(r, t). The source vector is denoted
s(r, t), and the source orientation is denoted η = [ηx , ηy , ηz ]T .
We thus have the relationship: s(r, t) = s(r, t)η .
Let us assume that a unit-magnitude source exists at r. When
this unit-magnitude source is directed in the x, y, and z directions, the outputs of the mth sensor are respectively denotedxyz(r), lm(r), and lm(r). Let us define an M × 3 matrix L(r)
by lm
z
y
x
(r)]. This matrix L(r),(r), lm(r), lm
whose mth row is equal to [lmreferred to as the lead field matrix, represents the sensitivity ofthe sensor array at r. When the unit-magnitude source at r isoriented in the η direction, the outputs of the sensor array are
expressed as l(r) = L(r)η . This column vector l(r), referred toas the lead field vector, represents the sensitivity of the sensor
array in the direction of η at the location r.
The outputs of the sensor array y(t) are expressed as thesum of the signal component yS (t) and the noise ε:y(t) = yS (t) + ε.
(1)In equation (1), yS (t) is called the signal vector, which isexpressed as:yS (t) =L(r)s(r, t) dr,(2)Ω

where the integral on the right-hand side is carried out over athree-dimensional volume Ω where signal sources of interest
can exist. This Ω is called the source space. In equation (1),an M × 1 random vector ε represents additive sensor noise,
which is assumed to obey the normal distribution:p(ε) = N (ε|0, 2 I),
(3)where I is the identity matrix and 2 is the variance of thesensor noise.
We denote the time series outputs of a sensor arrayy(t1 ), . . . , y(tK ), where K is the total number of measured time
points. It is assumed that K  >  M in this paper. We define themeasured data matrix B as:
2

Topical Review

J. Neural Eng. 14 (2017) 051001

B = [y(t1 ), . . . , y(tK )] = [y1 , . . . , yK ],(4)sq = [sq (t1 ), . . . , sq (tK )],(11)where y(tj ) is denoted yj for simplicity. We also define amatrix of the signal vector such thatwhich we call the time course vector of the qth source. We
then prove that a row of the signal matrix BS is expressed as alinear sum of the time course vectors, s1 , . . . , sQ . We assume,in this paper, that the source time course vectors s1 , . . . , sQ
are linearly independent. Substituting equation (8) into equation (5), the following relationship is obtained:

BS = [yS (t1 ), . . . , yS (tK )] = [yS1 , . . . , ySK ],(5)where the jth column of BS is denoted ySj. This BS is called
the signal matrix in this paper. Then, the data model in equation (1) is expressed in matrix form as:

B = BS + Bε ,(6)BS = 

where Bε is the noise matrix whose jth column is equal to thenoise vector ε at time tj.

Let us assume that a total of Q discrete sources exist. Theirlocations are denoted by r1 , . . . , rQ , their orientations byη 1 , . . . , η Q , and their magnitudes by s1 (t), . . . , sQ (t). Then,the source distribution is expressed as:q=1sq (tK )lq 

1q=1 [sq (t1 ), . . . , sq (tK )]lq



 Q

1
q=1 lq sq




 
..
..
,
=
.
.

 

Q
Q
M
M
q=1 lq sq
q=1 [sq (t1 ), . . . , sq (tK )]lq

(12)

Q

β Sj =
lqj sq .
(13)

q=1

q=1

where δ(r) indicates the delta function. Substituting theequation above into equation (2), the signal vector yS (t) isexpressed as:

Q
Q


y
(t)
=
L(r)
s
(t)η
δ(r
−
r
)
dr
=
sq (t)lq ,
(8)
q
q
q
S
q=1

sq (t1 )lq , . . . ,Qwhere lq1 , . . . , lqM are the elements of the lead field vector lq :lq = [lq1 , . . . , lqM ]T . Denoting the jth row vector of BS by β Sj ,equation (12) shows that

Q

s(r,
t)
=
sq (t)η q δ(r − rq ),
(7)

Ω

q=1

 Q


=


2.2. Definition of signal subspace in the spatial domain

QThis equation indicates that a row vector of the signal matrix,
β Sj , is expressed as a linear sum of sq (q = 1, . . . , Q). That is,
we haveβ Sj ∈ rsp( [sT1 , . . . , sTQ ]T ),
(14)

q=1

where lq represents the lead field vector of the qth source
obtained such that lq = L(rq )η q. We assume that the numberof sources Q is smaller than the number of sensors, i.e. Q  <  M.
This assumption is referred to as the low-rank signal assumption [5, 15], and we hold this assumption throughout the paper4.
Equation (8) claims that the signal vector yS is expressed asa linear sum of the lead field vectors l1 , · · · , lQ. That is, the signal vector yS lies within a subspace spanned by l1 , · · · , lQ. The
subspace spanned by the source lead field vectors l1 , · · · , lQ is
defined as the signal subspace [5], which is denoted by ES, i.e.

where the notation rsp(·) indicates a row space of the matrix
in the parentheses.
Analogous to equations (9) and (10), it is reasonable
to define rsp( [sT1 , . . . , sTQ ]T ) as the signal subspace in time
domain KS , i.e.
KS = rsp( [sT1 , . . . , sTQ ]T ).
(15)By defining the time domain signal subspace this way, we canderive symmetric relationships between the time domain signal subspace and the spatial domain signal subspace. That is,
we have already shown the relationships:

ES = csp( [l1 , · · · , lQ ] ).
(9)Here, the notation csp(·) indicates the column space of thematrix within the parentheses. Equation (8) indicates therelationship,column of BS : ySj ∈ ES ,(16)yS (t) ∈ ES .
(10)row of BS : β Sj ∈ KS .
(17)The signal vector lies within the signal subspace, which is thesubspace formed by all possible signal vectors [1].

We can show that, with the assumption K  >  Q, the column
space of BS is equal to the spatial domain signal subspace, i.e.

2.3. Definition of signal subspace in the time domain

ES = csp(BS ).
(18)This section introduces a novel definition of signal subspacein the time domain. To do so, we define a row vector sq consisting of the time course of the qth source such that

The proof is presented in appendix A.1. With the assumption
M  >  Q, the row space of BS is equal to the time domain signal
subspace, i.e.

4KS = rsp(BS ).
(19)Since we assume that K  >  M, the assumption K  >  M  >  Q holds throughout
the paper.

The proof is presented in appendix A.2.
3

Topical Review

J. Neural Eng. 14 (2017) 051001

3. Interference removal using the spatial-domain
signal subspacey(t) = yS (t) + yI (t) + ε,(23)where yI (t) represents the interference overlapped on the signal vector yS (t). In this paper, this yI (t) represents interferences originated from outside the source space. The signal
vector yS (t) represents all signals originated from inside thesource space6. We assume that a total of P sources generate
the interference, and yI (t) is thus expressed as

3.1. Estimation of spatial-domain signal subspace

Although the signal subspace is defined in equation (9), it isdifficult to use this equation to derive the signal subspace,because the source lead field vectors are generally unknown5.
The signal subspace can be estimated using the time series
measurement of the sensor data y(t) as described in thefollowing.
According to equation (18), we can estimate the signal subspace through the estimation of the column space of BS . For
this estimation, let us express the singular value decomposition of BS (M  <  K) as:
BS =

M


P

yI (t) =σp (t)ξ p ,
(24)
p=1

where ξ p is the lead field vector of the pth interference sourcewith its amplitude of σp (t). Note that we still assume that the
row-rank signal assumption P  +  Q  <  M holds. According toequation (9), the interference subspace EI is defined as

γjS uSj (vSj )T .

EI = csp( [ξ 1 , . . . , ξ P ] ).
(25)

j=1

The interference subspace can be estimated when the datathat contain only interference are available. Such data arecalled the control data in this paper, and expressed as

Here, the singular values γjS ( j = 1, . . . , M ) are numbered indecreasing order, and uSj and vSj are the spatial and temporal singular vectors, respectively. Since BS is a matrix whose
rank is equal to Q, the singular values of BS are given as
γ1S , . . . , γQS , 0, . . . , 0. Namely, BS has only Q non-zero singular values, and the column space of BS is equal to the span of
singular vectors corresponding to the non-zero singular values, uS1 , . . . , uSQ. That is, we have the relationshipand the control data matrix is defined as
Bc = [yc (t1 ), . . . , yc (tK )]. The spatial singular vectors of Bc
that correspond to the P largest singular values are denotedby uc1 , . . . , ucP. The projector onto the interference subspace,
PI , is formulated such that

ES = csp( [uS1 , . . . , uSQ ] ).
(20)PI = [uc1 , . . . , ucP ][uc1 , . . . , ucP ]T .
(27)yc (t) = yI (t) + ε,(26)Since the signal matrix BS (and thus its singular vectors) are
unknown quantities, we cannot use equation (20) to derive thesignal subspace.
The singular vectors uS1 , . . . , uSQ are estimated using the singular vectors of B . Let us denote the spatial singular vectorsof the data matrix B that correspond to the Q largest singularvalues as u1 , . . . , uQ. Then, using the noise model in equation (3), singular vectors u1 , . . . , uQ are asymptotically equal
to uS1 , . . . , uSQ. That is, in the limit of infinite number of timesamples, we haveSince the interference vector yI (t) can be expressed as
PyI (t) = j=1 aj (t)ucj , the relationship (I − PI )yI (t) = 0
holds. Thus, by projecting the data vector onto the subspaceorthogonal to the interference subspace, we have:yS (t) = (I − PI )y(t) = yS (t) − PI yS (t) + (I − PI )ε.
(28)It is apparent from equation (28) that the projector I − PI
can remove the interference yI (t) but also affects the signal.
The influence on the signal vector is evaluated by the second
term on the right-hand side. This term is generally small when
the orthogonality between lead field vectors of signal sourceslj ( j = 1, . . . , Q) and the basis vectors of the interference subspace uci (i = 1, . . . , P) is high. The method of interference
suppression based on equation (28) is called the signal space
projection (SSP) [2, 4]. One problem we have when implementing SSP is the determination of P, the dimension of theinterference subspace. This P is usually determined by thresholding the singular values. It is obvious from the arguments
above that the underestimation of P results in a stuation that apart of yI (t) remains in yS (t) and conversely the overestimation may cause a large distortion of the signal vector yS (t) inyS (t).

ES = csp( [u1 , . . . , uQ ] ).
(21)The proof is presented in appendix A.4. In the case of a finite
number of time samples, csp( [u1 , . . . , uQ ] ) is equal to the
maximum likelihood estimate of the signal subspace ES:
(22)ES = csp( [u1 , . . . , uQ ] ).

A formal proof of equation (22) can be found in [15] and[16].
3.2. Signal space projection (SSP) algorithm

The signal space projection (SSP) is an algorithm intended toremove the interference overlapped onto the signal [2, 4]. The
algorithm is based on the theory of signal subspace and themeasurement model is

6
Therefore, yS (t) includes interferences originated from inside the source
space; such interferences include so called brain noise in MEG measurements. In other words, we consider the brain noise as a part of the signal andits removal is not argued in this paper.

5

To derive the source lead field vectors, we must know the locations andorientations of sources. These quantities are generally unknown.
4Topical Review

J. Neural Eng. 14 (2017) 051001

The equation above indicates that the SSS signal extractor ΓS
passes the signal-source lead field vector lq with no distortion.
The assumption that all the interference sources are locatedwithin the external region leads to the situation that the column space of D includes the interference subspace, i.e.

3.3. Signal space separation (SSS) algorithm

A method for deriving the signal subspace based not on the
span of source lead field vectors but on physical properties ofthe magnetic field predicted by Maxwell’s equations has been
proposed [12–14, 17]. In this method, assuming that there areno sources in a region where sensors are located (a region
called the sensor region), the sensor measurements can be
expressed by an expansion using the vector spherical harmonics, which contains natural separation between the magnetic
fields generated from the internal and external regions. Here,the internal region indicates the region closer to the origin
than the sensor region, and the external region indicates theregion farther from the origin than the sensor region.
Therefore, if the origin is properly set such that the source
space Ω is included within the internal region, and the interference sources are located within the external region, the signal
can be separated from the interference. Interference removal
based on this idea is called signal space separation (SSS). The
derivation of the SSS method is presented in the appendix B.
With the proper setting of the origin, we can derive the SSS
signal extractor ΓS:

csp(D) ⊃ EI .

Therefore, the lead field vector of an interference source ξ plies within the column space of D, and ξ p is expressed as
ND

ξp =
φj dj = Dφ,
(33)
j=1

where dj is the jth column of D, φj is the jth expansion
coefficient, and φ is a column vector containing these coefficients: φ = [φ1 , . . . , φND ]T . Again denoting the NC × 1 column vector whose elements are all zero by 0, we have therelationshipΓS ξ q = ΓS Dφ = ΓS S

C

(34)
The equation above indicates that the signal extractor ΓS completely blocks the lead field vector of an interference sourceξ q. Consequently, substituting (32) and (34) into (30), we
obtainΓS = CCT (CCT + DDT )−1 = CCT (SST )−1 ,(29)where the matrices C and D are given respectively in equations (B.9) and (B.10), and S = [C, D]. Note that C and D are
M × NC and M × ND matrices where NC and ND are defined inequation (B.14). The derivation of ΓS is presented in appendix
B.2. We use the data model in equation (23), and assume that
equation (8) holds for yS (t) and equation (24) for yI (t). Then,
multiplying the signal extractor ΓS to the data vector y(t) givesΓS y(t) = ΓS yS (t) + ΓS yI (t) =Qsq (t)ΓS lq +

q=1

P


ΓS y(t) = ΓS yS (t) + ΓS yI (t) =+sq (t)ΓS lq

σp (t)ΓS ξ p =Qsq (t)lq = yS (t).

(35)

q=1


The equation above shows that by multiplying the signal
extractor ΓS with the data vector y(t), the signal vector yS (t)is selectively extracted with no distortion. This distortionlesssignal extraction is a major advantage of the SSS algorithm
over the SSP algorithm.

σp (t)ΓS ξ p ,p=1

4. Time-domain signal space projection

Thus, the lead field vector of a signal source lq lies within thecolumn space of C , and lq is expressed as a linear sum of thecolumn vectors of C , such that

If we obtain the basis vectors of the interference subspace inthe time domain, it is possible to remove the interference byprojecting the measured data onto the subspace orthogonal to
the time-domain interference subspace. We define the interference matrix BI as

NC

l
=
αj cj = Cα,
(31)
q
j=1

BI = [yI (t1 ), . . . , yI (tK )].
(36)where cj is the jth column of C , αj is the jth expansion coefficient, and α is a column vector containing the coefficients, i.e.
α = [α1 , . . . , αNC ]T . Thus, denoting an ND × 1 column vector
whose elements are all zero by 0, and using the derivation inequation (B.20), we can derive the relationship

Then, the data model in equation (23) is expressed as
B = BS + BI + Bε .
(37)Let us define the time course vector of the pth interference
source, σ p, as

 

 
α
α
T
T −1
ΓS lq = ΓS Cα = ΓS S
= C C (SS ) S
0
0

 

 
α
α
= C ST (SST )−1 S≈ C (ST S)−1 ST S
0 [1:N ]0 [1:N ]= Cα = lq .

P

p=1

csp(C) ⊃ ES .



Q

q=1

(30)
where the sensor noise term is dropped for simplicity.
The assumption that the source space Ω is located within
the internal region leads to the relationship

C

 

 
0
0
≈ C (ST S)−1 ST S
= C0 = 0.
φ
φ [1:N ]σ p = [σp (t1 ), . . . , σp (tK )].
(38)Then, the interference subspace in the time domain, KI , isdefined as

C

KI = rsp( [σ T1 , . . . , σ TP ]T ).
(39)

(32)
5

Topical Review

J. Neural Eng. 14 (2017) 051001

We assume that time course vectors σ T1 , . . . , σ TP are linearly
independent in this paper.
If basis row vectors of the interference subspace in thetime domain are estimated as ψ 1 , . . . , ψ P, defining Υ as
Υ = [ψ T1 , . . . , ψ TP ]T , we can compute the projector onto KI
such thaty(t) = ZyR (t) + d(t),(44)where Z is an M × J coefficient matrix of this multi-variate
regression, and the residual signal d(t) represents the interference removed results. Here, the coefficient matrix Z is
obtained by solving the minimization problem:
Z = argminy(t) − ZyR (t)2 ,(45)ΠI = ΥT (ΥΥT )−1 Υ.
(40)

Z

where · indicates time average. That is, Z is determined so
as to maximize the correlation between y(t) and yR (t). Theinterference-removed results are expressed as:

We denote the jth row of BI by β Ij . Since β Ij is expressedP
as the linear sum of ψ j: β Ij = j=1 cj ψ j , we thus haveβ Ij (I − ΠI ) = 0. Therefore, projecting the data matrix B onto
S
the subspace orthogonal to KI , the estimated signal matrix B
is given by S = B(I − ΠI ) = (BS + BI + Bε )(I − ΠI )B
= BS − BS ΠI + Bε (I − ΠI ).

 R (t) = y(t) − y(t)yTR (t) yR (t)yTR (t) −1 yR (t).
d(t) = y(t) − Zy(46)
To understand the relationship between ANC and the timedomain SSP, Let us rewrite equation (46) using a matrix form.
To do so, the matrix of the reference sensor data is defined as
BR, such that

(41)

The method of removing the interference BI based onequation (41) is referred to as the time-domain signal space
projection (time-domain SSP). The influence of the time
domain SSP on the signal component is assessed by the second term BS ΠI on the right-hand side of equation (41). This
term becomes small when the correlations between the timecourses of the signal and interference sources are small. Thiscan be considered an advantage of the time-domain SSP overthe spatial-domain SSP. This is because in many real-life
applications, the time courses of the signal and interference
sources are expected to differ significantly, but the orthogonality of lead field vectors between signal and interference
sources may not be so high. In the next section, we show thata number of existing interference removal methods can beinterpreted as the time domain SSP, and that these methods
differ only in their manner of deriving the basis vectors of thetime-domain interference subspace.

BR = [yR (t1 ), . . . , yR (tK )].

Using the data matrices of BR and B , we have11(47)
y(t)yTR (t) = BBTR , and yR (t)yTR (t) = BR BTR .
KKS :Thus, denoting the interference removed results as B
 S = [d(t1 ), . . . , d(tK )], equation (46) is rewritten asB S = B − BBTR BR BTR −1 BR
B

−1 

BR = B [I − ΠR ] ,= B I − BTR BR BTR(48)
where ΠR indicates the projector onto the row space of BR
defined such that
−1

(49)
BR .
ΠR = BTR BR BTR

Comparison between equations (41) and (48) makes it
clear that ANC is equivalent to the time-domain SSP if therelationship ΠI ≈ ΠR holds where ΠI is the projector onto
the interference subspace defined in equation (40). Actually,
if the sensor noise is small, the row space of BR can reasonably approximate the row space of BI , which is equal to theinterference subspace under the low rank signal assumptionJ  >  Q. However, when the sensor noise is not negligibly smallor when the reference sensor data contains components notcontained in the measurement sensor data, ΠR and ΠI may
have some difference, and results from the ANC algorithm
may contain errors.

5. Interference removal based on the time-domain
SSP
5.1. Adaptive noise canceling(ANC)

Adaptive noise canceling(ANC) is an interference removal
method which makes use of data from reference sensors. The
reference sensors collect data containing interference but not
the signal of interest [6, 7, 18]. It is assumed that the sensor
array is equipped with a total of J additional reference sensors,
and the outputs of the reference sensors are denoted by a J × 1
column vector yR (t). Then, the data model is expressed asy(t) = yS (t) + yI (t) + ε,(42)5.2. Common temporal subspace projection (CTSP)yR (t) = yI (t) + 
ε.
(43)

Common temporal subspace projection (CTSP) also removes
the interference by making use of the reference sensor data [9].
The difference between CTSP and ANC is that CTSP assumes
that the reference sensor data contain components that exist
only in the reference sensor data but not in the measurement
sensor data. The data model assumed in CTSP is expressed as

Note that y(t), yS (t), yI (t) and ε are M × 1 column vectors,
and yR (t), yI (t), and ε are J × 1 column vectors. ANC tries toremove the interference yI (t) from the data vector y(t) by taking out components maximally correlated with the reference
sensor data yR (t). This removal is carried out by regressing
y(t) with yR (t), i.e.

B = BS + BI + Bε ,(50)6Topical Review

J. Neural Eng. 14 (2017) 051001

between two subspaces is described in appendix A.9. Using
this algorithm, the orthonormal basis vectors of the intersection rsp(B) ∩ rsp(BR ) can be obtained. Denoting these basis
vectors as ψ 1 , . . . , ψ r , we can remove the interference and S such thatobtain the estimated signal matrix B

 I + Bw + B
ε,
(51)
BR = B

where a J × K matrix Bw indicates the components that arecontained only in the reference sensor data, and J × K matri I and B ε, are the interference and noise matrices of theces, Breference sensor data. Here, we assume Bw BTS = 0 . According
to the arguments in appendix A.5, equation (50) leads to therelationship, S = B(I − [ψ 1 , . . . , ψ r ][ψ 1 , . . . , ψ r ]T ).
B
(59)The time domain SSP in equation (59) cannot perfectly
remove the interference because the basis vectors ψ 1 , . . . , ψ rspan only a part of the interference subspace. Nonetheless, we
can still expect that the method can reduce the interference
if rsp(B) ∩ rsp(BR ) is a reasonable approximation of KI . The
algorithm that performs the interference removal in a manner
described above is called common temporal subspace projection (CTSP) [9].

rsp(B) = rsp(BS + BI + Bε ) ⊂ rsp(BS )(52)+ rsp(BI ) + rsp(Bε ) = KS + KI + Kε ,
where Kε indicates Kε = rsp(Bε ). Equation (51) leads to I + Bw + B ε ) ⊂ rsp(BI )rsp(BR ) = rsp(Bε, ε ) = K I + Kw + K(53)+ rsp(Bw ) + rsp(B ε indicate Kw = rsp(Bw ) and K
 ε = rsp(B ε ).
where Kw and K
 I ) = KI
The relationships rsp(BS ) = KS and rsp(BI ) = rsp(Bare used here. Then, according to the arguments in appendix
A.6, equations (52) and (53) lead to the following relationships among sets of basis vectors:

5.3. Dual signal subspace projection (DSSP)
5.3.1. Pseudo signal subspace projector. Dual signal sub-space projection (DSSP) removes the interference without
using either the reference sensor data or control data suchas the empty-room-noise data [11]. The algorithm assumes
the data model in equation (50) with the assumption that theinterference sources are located outside the source space. The
DSSP algorithm uses the so-called pseudo signal subspace
projector, and to derive it, voxels are defined over the source
space Ω, in which the voxel locations are denoted r1 , . . . , rN .
The pseudo signal subspace ĔS is defined such thatS B ⊂ SS ∪ S I ∪ S ε
(54)SBR ⊂ SI ∪ Sw ∪ Sε(55)where sets of basis vectors of rsp(B) and rsp(BR ) are respectively denoted by SB , SBR , i.e. rsp(B) = span(SB ) andrsp(BR ) = span(SBR )7. Also, sets of basis vectors of KS , KI ,
 ε are respectively denoted by SS , SI , Sw , Sε,
Kw , Kε, and K
and Sε.
We assume that the noise time courses are orthogonal to each other, and they are orthogonal to the signal andinterference time courses, resulting in the relationships:
Sε ∩ Sε = ∅,Sε ∩ SS = Sε ∩ SI = Sε ∩ Sw = ∅,andSε ∩ SS = Sε ∩ SI = Sε ∩ Sw = ∅. Here ∅ indicates theempty set. We also use SS ∩ Sw = ∅ , which results from theorthogonality assumption Bw BTS = 0 . Then, using the proof
in appendix A.7 and the distributive property, we can deriveĔS = csp( [L(r1 ), . . . , L(rN )] ).
(60)If the voxel interval is sufficiently small and voxel discretization errors are negligible, we have the relationship ĔS ⊃ ES .
Therefore, a vector contained in the signal subspace is also contained in the pseudo signal subspace, i.e. if x ∈ ES, then x ∈ ĔS.
Let us derive the orthonormal basis vectors of thepseudo signal subspace ĔS. To do so, defining a matrix F
as F = [L(r1 ), . . . , L(rN )], we compute the singular value
decomposition of F ,MF=
λj ej f Tj ,(61)SB ∩ SBR ⊂ (SS ∪ SI ∪ Sε ) ∩ (SI ∪ Sw ∪ Sε )= [SS ∩ (SI ∪ Sw ∪ Sε )] ∪ [SI ∩ (SI ∪ Sw ∪ Sε )] ∪ [Sε ∩ (SI ∪ Sw ∪ Sε )]j=1

S I ⊃ S B ∩ S BR
(57)where we assume the relationship M  <  N. If the singularvalues λ1 , . . . , λτ are distinctively large, and other singularvalues λτ +1 , . . . , λM are nearly equal to zero, the singular vectors e1 , . . . , eτ forms orthonormal basis vectors of the pseudo
signal subspace ĔS and the projector onto ĔS is obtained using

KI = span(SI ) ⊃ span(SB ∩ SBR ) = span(SB ) ∩ span(SBR )= rsp(B) ∩ rsp(BR ).
(58)

Note that, since ĔS ⊃ ES , the orthogonal projector (I − P̆S )projects out the signal vector, i.e. (I − P̆S )yS (t) = (I − P̆S )BS = 0.

= (SS ∩ SI ) ∪ SI ∪ ∅ = (SS ∩ SI ) ∪ SI = SI .


Namely,(56)

holds. According to the arguments in appendix A.8, we canobtain the relationship,P̆S = [e1 , . . . , eτ ][e1 , . . . , eτ ]T .
(62)The equation above claims that the intersection betweenrsp(B) and rsp(BR ) forms a subset of the interference subspace
KI . An algorithm that derives basis vectors of the intersection

5.3.2. DSSP algorithm. The DSSP algorithm applyies P̆S
and I − P̆S to the data matrix B to create two kinds of data
sets:

7

(63)
P̆S B = BS + P̆S BI + P̆S Bε ,Here, an expression of X = span(SX ) indicates that a set of basis vectors
SX spans the subspace X .
7

Topical Review

J. Neural Eng. 14 (2017) 051001

interference subspace KI . The basis vectors of the intersection are derived using the algorithm described inappendix A.9. Once the orthonornal basis vectors of theintersection rsp(P̆S B) ∩ rsp((I − P̆S )B) are obtained, we
can compute the projector onto the intersection, andtime-domain SSP can be implemented. If the intersection
rsp(P̆S B) ∩ rsp((I − P̆S )B) is a reasonable approximation of
KI , this time-domain SSP will be able to remove interferences effectively. The method of removing the interference
in a manner as described above is called dual signal subspace projection (DSSP) [11].
Note that since the intersection rsp(P̆S B) ∩ rsp((I − P̆S )B)
is only a subset of the interference subspace KI , the method
cannot perfectly remove interferences. However, we canobserve that rsp(P̆S B) ∩ rsp((I − P̆S )B) becomes a better
approximation of KI , if the interference BI is significantly
greater than the signal BS (and the sensor noise Bε), resulting in a situation that the interference terms are dominated inP̆S B and (I − P̆S )B. This may explain our empirical findings,
which are a little counter-intuitive, that the method works better for larger interferences.

(I − P̆S )B = (I − P̆S )BI + (I − P̆S )Bε .
(64)In equation (64), the signal of interest is suppressed by multiplying (I − P̆S ) with the data matrix, and a virtual reference time series is estimated8. Here, we use P̆S BS = BS and(I − P̆S )BS = 0. According to the arguments in appendix A.5,
the following relationships hold:rsp(P̆S B) ⊂ rsp(BS ) + rsp(P̆S BI ) + rsp(P̆S Bε ),(65)rsp((I − P̆S )B) ⊂ rsp((I − P̆S )BI ) + rsp((I − P̆S )Bε ).
(66)The arguments in appendix A.3 prove that the relationships,rsp(P̆S BI ) = KI and rsp((I − P̆S )BI ) = KI hold. Thus, using
rsp(BS ) = KS , equations (65) and (66) lead torsp(P̆S B) ⊂ KS + KI + K̆ε ,(67)rsp((I − P̆S )B) ⊂ KI + K̆ε ,(68)where we use the notations, rsp(P̆S Bε ) = K̆ε andrsp((I − P̆S )Bε ) = K̆ε . Let us define the sets of basis
vectors of rsp(P̆S B) and rsp((I − P̆S )B) as SPS and SPI ,(i.e. rsp(P̆S B) = span(SPS ) and rsp((I − P̆S )B) = span(SPI )).
According to the arguments in appendix A.6, equations (67)and (68) can be converted into relationships among the sets ofbasis vectors, which are written as

5.4. Spatio-temporal signal space separation (tSSS)

In section 3.3, we argue that one prerequisite of the signal
space separation (SSS) algorithm is that all interference
sources be located in the external region. However, thisrequirement is not always fulfilled, as interference sources canalso be located fairly close to the source space. In this case, the
SSS extractor ΓS cannot adequately remove the interference.
Spatio-temporal signal space separation algorithm (tSSS) has
been developed for such situations [10].
This algorithm first applies the SSS extractors ΓS and ΓI to
the data matrix B to create two kinds of data sets:

SPS ⊂ SS ∪ SI ∪ S̆ε ,(69)SPI ⊂ SI ∪ S̆ε ,(70)where basis vector sets of K̆ε and K̆ε are denoted S̆ε and S̆ε .
Since the relationshipP̆S Bε ((I − P̆S )Bε )T = P̆S Bε BTε (I − P̆S ) = ρ2 P̆S (I − P̆S ) = 0,

holds, S̆ε and S̆εhas no intersection, i.e.
S̆ε ∩ S̆ε = ∅. Also, assuming that BS BTε = BI BTε = 0, the relationships S̆ε ∩ SS = S̆ε ∩ SI = ∅ and S̆ε ∩ SS = S̆ε ∩ SI = ∅hold. Therefore, using the proof in appendix A.7 and the distributive property, we can derive the following relationship: SPS ∩ SPI ⊂ SS ∪ SI ∪ S̆ε ∩ SI ∪ S̆ε

 
 

= SS ∩ (SI ∪ S̆ε ) ∪ SI ∩ (SI ∪ S̆ε ) ∪ S̆ε ∩ (SI ∪ S̆ε )ΓS B = BS + ΓS BI + ΓS Bε ,(74)ΓI B = ΓI BI + ΓI Bε ,(75)where we use the relationships, ΓS BS = BS , and ΓI BS = 0 . We
thus obtain,
rsp(ΓS B) ⊂ KS + rsp(ΓS BI ) + rsp(ΓS Bε ),
(76)

(71)

rsp(ΓI B) ⊂ rsp(ΓI BI ) + rsp(ΓI Bε ),(77)S I ⊃ S PS ∩ S PI .
(72)where we use rsp(BS ) = KS .
When an interference source is located close to the source
space, the lead field vector of this source, ξ p, may have comp­onents expanded by the columns of C , as well as components
expanded by the columns of D, resulting in,= (SS ∩ SI ) ∪ SI ∪ ∅ = (SS ∩ SI ) ∪ SI = SI .

Namely, we have

According to the arguments in appendix A.8, we derive

KI = span(SI ) ⊃ span(SPS ∩ SPI ) = span(SPS ) ∩ span(SPI )= rsp(P̆S B) ∩ rsp((I − P̆S )B).

ξ p = Cα + Dφ ,
(78)

(73)


The quation above shows that the intersection betweenrsp(P̆S B) and rsp((I − P̆S )B) forms a subset of the8

where α and φ are vectors containing the expansion coefficients. Therefore, by multiplying the SSS signal extractor toξ p, we haveΓS ξ p = ΓS Cα + Dφ = Cα = π p ,(79)This procedure is called the generalized sidelobe canceller in [19].

8Topical Review

J. Neural Eng. 14 (2017) 051001

Figure 1. An arrangement of the 275-channel whole-head sensor array of the OmegaTM (VMS Medtech, Coquitlam, Canada)

neuromagnetometer used in our numerical experiments. Locations of sensors are indicated by filled circles. Locations and orientations of
six reference sensors consisting of two sets of vector magnetometers are shown. Here, (a), (b) and (c), respectively, show coronal, axial,
and sagittal views of the sensor arrangement. Note that the arrangement of the reference sensors was assumed solely for our numericalexperiments, and differs from the true arrangement in the Omega neuromagnetometer system.
Table 1. Locations of interference sources assumed in numericalexperiments.

Source number

Location (cm)

1
2
3
4

(100, −30, −500)(5, −50, 600)(−27, −495, 21)(755, −342, 37)

12(10, 30, 30)(5, −20, −30)

Distance fromsignal source (cm)
510
600
500
830
41
39

However, these extractors never change the time courses
of interference sources. Thus, using the same arguments as
in appendix A.3, we can prove that rsp(ΓS BI ) = KI andrsp(ΓI BI ) = KI hold, and equations (76) and (77) becomersp(ΓS B) ⊂ KS + KI + Kε ,(81)Figure 2. (a) Time course of the signal source assumed in our
numerical experiments. The time t is expressed with the unit of time
point ranging from t  =  −1200 to 1200. The broken vertical line
indicates the time point at t  =  520. The field maps at this time point
are shown in the following figures. (b) Time courses of the fourinterference sources assumed in our numerical experiments.

rsp(ΓI B) ⊂ KI + Kε .
(82)Sets of basis vectors of rsp(ΓS B) and rsp(ΓI B) are respectively
denoted by SΓS and SΓI . Then, according to the arguments inappendix A.6, we can obtainSΓS ⊂ SS ∪ SI ∪ Sε ,(83)which shows that ΓS changes the interference lead field fromξ p to π p . Similarly, by applying ΓI to ξ p, we have p,ΓI ξ p = ΓI Cα + Dφ = Dφ = π(80)which shows that the lead field ξ p is changed to π p.

SΓI ⊂ SI ∪ Sε ,(84)where sets of basis vectors of KS , KI , Kε , and Kε are respectively
denoted by SS , SI , Sε , and Sε . Asumuing that relationships,
9Topical Review

J. Neural Eng. 14 (2017) 051001

Figure 3. (a) Time courses of the signal magnetic field with Gaussian noise with the signal-to-noise ratio (SNR) equal to 32. (b) The
magnetic field map of the signal magnetic field at t  =  520. (c) Time courses of interference-overlapped sensor data with the interference-tosignal ratio (ISR) equal to 6. (d) The magnetic field map of the interference overlapped sensor data at t  =  520.

Figure 4. (a) Simulated empty-room noise data used in the SSP interference removal experiments. (b) Sensor time courses of the SSP
interference-removal results. (c) Magnetic field map of the SSP interference removal results. The map at t  =  520 is shown.

10

Topical Review

J. Neural Eng. 14 (2017) 051001

SS ∩ Sε = SI ∩ Sε = ∅, and SI ∩ Sε = Sε ∩ Sε = ∅ hold,we can finally obtain
S I ⊃ SΓ S ∩ S Γ I ,(85)and according to the arguments in appendix A.8, we derive therelationship:
KI ⊃ rsp(ΓS B) ∩ rsp(ΓI B).
(86)The basis vectors of the intersection rsp(ΓS B) ∩ rsp(ΓI B) are
obtained using the algorithm in appendix A.9, and they can beused for forming the projector onto the interference subspace.
Making use of the interference subspace projector obtainedin this manner, the tSSS algorithm performs the time-domain
SSP for interference removal. Again, although these basis
vectors span only a subset of the interference subspace, themethod can effectively remove the interference if the intersection rsp(ΓS B) ∩ rsp(ΓI B) is a reasonable approximation of KI .
One such case is that the interference terms are dominated in
the data sets ΓS B and ΓI B.

Figure 5. (a) Time courses of reference sensor data used in ANC
and CTSP interference-removal experiments of which results
are shown in figure 6. (b) Time courses of reference sensor datacontaining a large additive fluctuation. The data were used in ANC
and CTSP interference removal experiments of which results areshown in figure 7.

5.5. Sensor noise suppression (SNS)The sensor noise suppression (SNS) algorithm has been developed to suppress sensor noise [8]. The algorithm assumes thedata model

It can be seen in equation (93) that the sensor noise suppression algorithm is the time-domain SSP, assuming that rsp(Θj )approximates the time-domain signal subspace KS and the pro
jector Θj ΘTj Θj −1 ΘTj approximates the projector onto KS .
When the number of sensors are large, rsp(Θj ) should
reasonably approximate rsp(B). When the sensor noise isnegligibly small, rsp(B) then approximates rsp(BS ), which isequal to KS , resulting in the projector Θj ΘTj Θj −1 ΘTj thatapproximates the projector onto KS . However, when the sensor noise is not small, (which is exactly the case where the
SNS algorithm is needed,) rsp(B) differs from rsp(BS ), and
the projector Θj ΘTj Θj −1 ΘTj may differ from the signal
subspace projector, resulting in the low performance of the
SNS algorithm in such cases.

B = BS + Bε .
(87)A key assumption of the SNS algorithm is that β j, the jth row
of B , lies within a span of all the rows of B except the jthrow, i.e.
β j ∈ rsp( [β T1 , . . . , β Tj−1 , β Tj+1 , . . . , β TM ]T ).
(88)Thus, β j is expressed as the linear sum of the other rows:β Tj =
ωi β Ti = Θj ω j ,
(89)
i=j

where



Θj = β T1 , . . . , β Tj−1 , β Tj+1 , . . . , β TM ,(90)

6. Numerical examplesTω j = [ω1 , . . . , ωj−1 , ωj+1 , . . . , ωM ] .
(91)In the equations above, the notationi=j indicatesthe summation from i  =  1 to M, except i  =  j, and ωi(i = 1, . . . , j − 1, j + 1, . . . , M ) are weights of the linear
summation.
The optimum weight, ω j , can be obtained using leastsquares fitting:

6.1. Data generation

A series of computer simulations were carried out to provide
illustrative examples of the results of our arguments in thepreceding sections. A sensor alignment of the 275-channel
whole-head sensor array from the OmegaTM (VMS Medtech,
Coquitlam, Canada) neuromagnetometer was used. The configuration of the sensor array arranged on the helmet surface
is shown in figure 1, in which the filled circles indicate thelocations of sensors. In our experiments, all these sensors
were assumed to be magnetometers. We assumed six reference sensors, consisting of two sets of vector magnetometersin addition to the measurement sensors on the helmet, and thereference sensor arrangement is shown also in figure 19.

 j = argmin β Tj − Θj ω j 2 .
ωωjThe solution is−1 T T(92)
 j = ΘTj ΘjΘj β j .
ωThis ω
 j is substituted into equation (89), and the denoised row
β j is obtained as = β Θj ΘT Θj −1 ΘT .
(93)
β
j
j
j
j

9

Note that this arrangement of the reference sensors is assumed solely forour numerical experiments, and differs from the true arrangement in the
Omega neuromagnetometer system.

11

Topical Review

J. Neural Eng. 14 (2017) 051001

Figure 6. Results of ANC and CTSP interference removal experiments. The reference sensor data used in the experiments are shown infigure 5(a). (a) Sensor time courses of ANC interference removal results. (b) Magnetic field map of the ANC interference-removal results.
(c) Sensor time courses of CTSP interference removal results. (d) Magnetic field map of the CTSP interference-removal results. The fieldmaps at t  =  520 are shown for (b) and (d).

A single dipole source, which generated the signal magnetic
field, was assumed to exist at (0,0,10) with the orientation of(1,1,0). The location was 7 cm below the center of the sensor
array. The time course assigned to this source is shown in figure 2(a). In our numerical experiments, the time t is expressedwith the unit of time points ranging from t  =  −1200 to 1200.
We assumed four interference sources whose coordinates arelisted in the top four rows of table 1, which indicates that theseinterference sources were 500–800 cm away from the signal
source. The time courses assigned to these four interference
sources are shown in figure 2(b).
Sensor time courses were generated by projecting thesource time course in figure 2(a) onto the sensor timecourses through the lead field computed using Sarvas’ form­ula [20]. The sensor time courses of the signal magnetic
field (after adding sensor noise) are shown in figure 3(a).
The magnetic field map across sensors at t  =  520 is shownin figure 3(b). (The instant at t  =  520 is shown by the broken vertical line in figure 2(a).) Note that the sensor timecourses in figure 3(a) and the field map in figure 3(b) canwork as the ground truth in interference removal experiments described below.
Also, the time courses of the interference sources in figure 2(b) were converted to the sensor time courses of the interference using the lead field computed from the Biot–Savart
law [18, 20]. The use of the Biot–Savart law here is due toour assumption that these interference sources are non-biological sources. The interference magnetic field was overlapped

Figure 7. Results of ANC and CTSP interference removal
experiments when the reference sensor data shown in figure 5(b)were used. (a) Sensor time courses of ANC interference-removal
results. (b) Sensor time courses of CTSP interference-removal
results.

onto the signal magnetic field with the interference-to-signal
ratio (ISR) equal to 6. Here, the ISR is defined as BI /BS 
where X indicates the Frobenius norm of a matrix X . Theinterference-overlapped sensor time courses are shown in figure 3(c). Since the interference is much stronger than the signal magnetic field, the sensor time courses are dominated by12Topical Review

J. Neural Eng. 14 (2017) 051001

Figure 8. Results of SSS interference-removal experiments. (a) Sensor time cources of SSS interference removal results. (b) Magnetic field
map of the SSS interference removal results. The map at t  =  520 is shown

Figure 9. (a) Sensor time courses of interference-overlapped data generated with nearby interference sources. (b) Magnetic field map of
the interference-overlapped data in (a). (c) Sensor time courses of SSS interference removal results obtained from the sensor data in (a). (d)Magnetic field map of the interference removal results obtained from the sensor data in (a). The maps at t  =  520 are shown for (b) and (d).

the interference. The magnetic field map of the interferenceoverlapped sensor data is shown in figure 3(d).

resultant field map in figure 4(c) and the map of the signalonly magnetic field in figure 3(b) indicates that a considerable
amount of signal distortion occurs in the SSP interferenceremoved results.

6.2. Results of interference removal experiments

6.2.2. Experiments on adaptive noise canceling (ANC) andcommon temporal subspace projection (CTSP). Adaptive

6.2.1. Experiments on signal space projection (SSP). We
applied the signal space projecton (SSP) algorithm describedin section 3.2 to the interference-overlapped sensor data infigure 3(c). The SSP algorithm requires control data, (whichcontain only interference) to estimate the interference subspace. Such control data, generated using a set of different
interference-source time courses, are shown in figure 4(a).
The results of SSP interference removal are shown in figure 4(b), indicating that the interference has been largely
removed. The field map of the interference-removed sensor
data is shown in figure 4(c). The comparison between thenoise canceling (ANC) and common temporal subspace projection (CTSP) were applied to the same interference-overlapped sensor data in figure 3(c). In these experiments, thereference sensor data shown in figure 5(a) were used, and theinterference removal results are shown in figure 6. According
to the resultant sensor time courses in figures 6(a) and (c),these methods removed most of the interference. The maps ofthe interference-removed results in figures 6(b) and (d) indicate that almost no signal distortion was caused.
13

Topical Review

J. Neural Eng. 14 (2017) 051001

Figure 10. Results of tSSS interference removal experiments. (a) Sensor time courses of interference removal results obtained using theinterference-overlapped data in figure 3(c). (b) Magnetic field map of the interference removal results in (a). (c) Sensor time courses ofinterference removal results obtained using the nearby-interference data in figure 9(a). (d) Magnetic field map of the interference removal
results. The maps at t  =  520 are shown for (b) and (d).

the external region, was nearly perfectly satisfied in thesenumerical experiments.
We then moved two of the interference sources closer tothe signal source by assigning the coordinates in the last two
rows of table 1 to these sources. The distances between the
new locations of the two interference sources and that of thesignal source were approximately 40 cm, and the spatial separation between the signal and intereference sources could beinsufficient. The sensor data overlapped by the interference,
generated from the nearby interference sources, are shownin figure 9(a), and the field map is shown in figure 9(b). The
SSS algorithm was applied to these interference-overlapped
data, and the results are shown in figures 9(c) and (d). In theseexperiments, the SSS method failed to remove the intereference. This is because the two interference sources werelocated fairly close to the signal source, resulting in the violation of the prerequisite for the SSS method.

Next, we performed experiments in which the reference
sensor data contained, in addition to the interference, a fluctuation that did not exist in the data from the measurement
sensors arranged on the helmet. In these experiments, the reference sensor data shown in figure 5(b) were used, and theresults of ANC interference removal are shown in figure 7(a).
These results show that a fairly large amount of interference
remained due to the incorrect reference sensor data. The
results of CTSP interference removal are shown in figure 7(b).
Here, the interefernce was nearly completely removed despite
the fact that the reference sensor data contained an additional
component. This is because the data model of the CTSP algorithm allows components to exist only in the reference data
but not in the measured data.
6.2.3. Experiments on signal space separation (SSS). The

signl space separation (SSS) algorithm was tested on theinterference-overlapped sensor data in figure 3(c). This algorithm requires neither control noise data nor reference sensor
data. It solely relies on the spatial separation between signal
and intereference sources. The resultant interference-removed
sensor time courses are shown in figure 8(a), and the field mapof these sensor data at t  =  520 is shown in figure 8(b). These
results show that the SSS algorithm nearly perfectly removed
the interference. This is because the prerequisite of the SSS
method, the conditions that signal sources be located within
the internal region, and interference sources be located within

6.2.4. Experiments on spatio-temporal signal space separation (tSSS) and dual signal subspace projection (DSSP). The

spatio-temporal signal space separation (tSSS) algorithm wasapplied to the interference-overlapped sensor data in figure 3(c), and to the nearby interference data in figure 9(a).
The interference-removed results are shown in figure 10.
The results demonstrate that the tSSS algorithm can effectively remove the interferences not only caused by distantsources but also caused by nearby sources. The field maps in14

Topical Review

J. Neural Eng. 14 (2017) 051001

Figure 11. Results of DSSP interference removal experiments. (a) Sensor time courses of interference removal results obtained using

the interference-overlapped data in figure 3(c). (b) Magnetic field map the interference removal results in (a). (c) Sensor time courses ofinterference removal results obtained using the nearby-interference data in figure 9(a). (d) Magnetic field map of the interference removal
results in (c). The maps at t  =  520 are shown for (b) and (d).

Figure 12. Results of SNS denoising experiments. (a) Sensor time courses of signal plus Gaussian noise data generated for SNS denoisingexperiments. The signal-to-noise ratio (SNR) was set at 10. (b) Magnetic field map of the generated sensor data. (c) Sensor time courses ofdenoised results. (d) Magnetic field map of the denoised results. The maps at t  =  520 are shown for (b) and (d).
15

Topical Review

J. Neural Eng. 14 (2017) 051001

figures 10(b) and (d) demonstrate that almost no signal dist­ortion was caused in the interference removal process.
The dual signal subspace projection (DSSP) algorithm wasapplied to the sensor data in figure 3(c) and the data in figure 9(a). The interference-removal results are shown in figure 11. The results demonstrate that the DSSP algorithm isalso effective in the case of nearby interference sources, aswell as in the case of distant interference sources. The fieldmaps in figures 11(b) and (d) demonstrate that the DSSP interference removal results are free of signal distortion.

Assuming that x ∈ ES, the vector x is expressed as a linear
sum of the lead field vectors lq (q = 1, . . . , Q), such thatQx=ωq l q .
(A.1)
q=1

In order to show that x ∈ csp([ yS1 , . . . , ySK ]), we show that x inequation (A.1) is also expressed as a linear sum of ySj:Kx=
αj ySj .
(A.2)
j=1

6.2.5. Experiments on sensor noise suppression (SNS). The

Actually, substituting equation (8), into (A.2), and using equation (A.1), we get


Q
Q
Q
K
K





j
j

α
α
s
l
=
ωq l q ,sl=(A.3)jj qqq q

sensor noise suppression (SNS) algorithm was also tested
with noisy sensor data. Gaussian noise was added to thesignal magnetic data in figure 3(a). Here, the signal to noise
ratio (SNR) was set to 10 where the SNR was defined asBS /Bε . The time courses of resultant noisy sensor data
are shown in figure 12(a) and the field map is in (b). The
SNS denoised results are shown in figure 12(c). Comparison
between the sensor time courses before SNS denoising in figure 12(a) and after denoising in (c) shows a significant reduction in sensor noise through SNS denoising. The denoised
magnetic field map in figure 12(d) indicates that no observable distortion is caused.

j=1

q=1

j=1

q=1

where sq (tj ) is denoted sqj . Comparing the coefficients of thevector lq on the left and right sides gives following set of Q
linear equations:
α1 s11 + · · · + αK sK1 = ω1 ,α1 s12 + · · · + αK sK2 = ω2 ,
..
..
.
.

7. Summary

α1 s1Q + · · · + αK sKQ = ωQ .

This paper reviews subspace-based interference removal
methods with an emphasis on the time domain signal subspace. We first provide a concise review on the conventional
spatial-domain signal subspace and the signal space projection (SSP) method. We then extend the notion of signal subspace to the time domain by proposing a novel definition of
the time-domain signal subspace, and introduce time domain
signal space projection. We show that many existing interference removal methods can be interpreted as a form of the time
domain SSP. These methods include adaptive noise canceling, sensor noise suppression, common temporal subspace
projection, spatio-temporal signal space separation and therecently proposed dual signal subspace projection. We also
show that the difference between these methods is primarily due to their manner of deriving the interference subspace.
Numerical examples that illustrate the results of our arguments are provided.

Assuming that K  >  Q, for an arbitrary set of ω1 , . . . , ωQ, a
set of α1 , . . . , αK exists, and this fact leads to the conclusionthat equation (A.2) holds for an arbitrary x ( x ∈ ES). Thus,
the relationship csp([yS1 , . . . , ySK ]) ⊃ ES holds. Note that weassume that the source time courses s1j , . . . , sKj ( j = 1, . . . , Q)
are linearly independent.
We next show that ES ⊃ csp([yS1 , . . . , ySK ]). To do so,
we show that if x ∈ csp([ yS1 , . . . , ySK ]), x ∈ ES holds. If
x ∈ csp([yS1 , . . . , ySK ]), equation (A.2) holds. Substituting
equation (8) into equation (A.2), we obtain


Q
Q
K
K
K






x=
αj ySj =αjαj sqj  lq ,sqj lq =(A.4)
j=1

j=1

q=1

q=1

j=1

which shows that the vector x is expressed as a linear sum oflq (q = 1, . . . , Q), namely, x ∈ ES. This leads to the conclusionthat ES ⊃ csp([yS1 , . . . , ySK ]). In summary, since both directions
csp([yS1 , . . . , ySK ]) ⊃ ES , and ES ⊃ csp([yS1 , . . . , ySK ]) hold, therelationship ES = csp(BS ) holds.

ORCID iDs
Kensuke Sekihara

q=1

https://orcid.org/0000-0002-2048-3797

A.2. Proof of equation (19): KS = rsp(B S )KS = rsp(BS ) can be shown just as ES = csp(BS ) was shownin appendix A.1. We first show that if x ∈ KS , the relationshipx ∈ rsp([β S1 , . . . , β SM ]) holds. If x ∈ KS , the relationshipAppendix A. Supplementary mathematicalarguments
A.1. Proof of equation (18): ES = csp(B S )Qx=ωq sq .
(A.5)

We first prove that csp(BS ) = csp([yS1 , . . . , ySK ]) ⊃ ES . To
do so, we show that if x ∈ ES, x ∈ csp([yS1 , . . . , ySK ]) holds.

q=1

16

Topical Review

J. Neural Eng. 14 (2017) 051001

holds. If this vector x sartisfies x ∈ rsp([β S1 , . . . , β SM ]), x is
expressed as

A.4. Asymptotic equivalence between singular vectors of B
and B SMx=αj β Sj .
(A.6)Let us respectively denote sample covariance matrices of B ,BS , and Bε by R , RS , and Rε , which are obtained such that
R = K1 BBT , RS = K1 BS BTS , and Rε = K1 Bε BTε. Assuming that
BS and Bε are uncorrelated, from equation (6) we have

j=1

Using equation (13) and the arguments in the preceding
section, under the assumption M  >  Q, for an arbitrary setof ω1 , . . . , ωQ, a set of α1 , . . . , αM exists. Therefore, equation (A.6) holds, and rsp([β S1 , . . . , β SM ]) ⊃ KS also holds.
We next show that if x ∈ rsp([β S1 , . . . , β SM ]), then x ∈ KS
holds. If x ∈ rsp([β S1 , . . . , β SM ]), equation (A.6) holds.
Substituting equation (13) into equation (A.6), we get


Q
Q
M
M
M






x=
αj β Sj =αjαj lqj  sq .
lqj sq =(A.7)
j=1

j=1

q=1

R = RS + Rε .
(A.10)Assuming the noise model in equation (3) and infinite numberof data samples, the noise covariance matrix is given by
(A.11)Rε = 2 I.
Sample covariance matrices RS and R are expressed using the
eigen decomposition as:QRS =(γjS )2 uSj (uSj )T ,
(A.12)

j=1

q=1

j=1

The equation above shows that x ∈ KS holds. Threrefore,KS ⊃ rsp([β S1 , . . . , β SM ]) holds. In summary, since both directions rsp([β S1 , . . . , β SM ]) ⊂ KS and KS ⊂ rsp([β S1 , . . . , β SM ])hold, the relationship KS = rsp(BS ) holds.

MR=γj2 uj uTj .
(A.13)
j=1

Substituting equations (A.12), (A.13) and (A.11) into equation (A.10), we can obtain

A.3. Proof of KI = rsp(XB I ) where X = P̆S or X = (I − P̆S )MUsing the same arguments as in the preceding section,
we can prove that KI = rsp(XBI ) holds where X = P̆S or
X = (I − P̆S ) and P̆S is the the pseudo signal subspace projector defined in equations (62). Using equation (24), multiplying X with the interference vector yI (t) gives

j=1



=

P

p=1

 P




=


σp (t1 )ξ̆ p , . . . ,(γjS )2 uSj (uSj )T + 2

j=1

M


uSj (uSj )T .

j=1

γj2 = (γjS )2 + 2 for j = 1, . . . , Q,(A.15)γj2 = 2 for j = Q + 1, . . . , M,(A.16)p=1

uj = uSj for j = 1, . . . , M.
(A.17)

where we use the notation ξ̆ p = Xξ p . That is, the matrix X
changes ξ p to ξ̆ p, but the multiplication by X never affects
the time course σp (t). Therefore, we get the following
relationship:


P
P


XBI = X σp (t1 )ξ p , . . . ,σp (tK )ξ p 
p=1

Q


(A.14)
Therefore, we have the following asymptotic relationships:

P
P


Xy
(t)
=
σ
(t)Xξ
=
σp (t)ξ̆ p ,p
(A.8)pI
p=1

γj2 uj uTj =

p=1

P

p=1

Equation (A.17) indicates the asymptotic equivalence between
the singular vectors of B and BS .
A.5. Row space of sum of matrices

Defining two subspaces as X and Y , the sum of X and Y is
defined such that [21],

σp (tK )ξ̆ p ˘p

p=1 [σp (t1 ), . . . , σp (tK )]ξ1

..
.
P
˘p
p=1 [σp (t1 ), . . . , σp (tK )]ξM

 P

 
=
 

˘p
p=1 ξ1 σ p

P

..
.

˘p
p=1 ξM σ p

X + Y = {x + y | x ∈ X and y ∈ Y}.
(A.18)That is, the sum of X and Y is defined as a set of all possible
sums between x ( x ∈ X ) and y ( y ∈ Y ). Then, we show that,for two arbitrary matrices X and Y , the relationship




,


rsp(X + Y) ⊂ rsp(X) + rsp(Y).
(A.19)holds. To show this, let us denote row vectors of X , Y , andZ = X + Y by x, y, and z, respectively. By definition, for z( z ∈ rsp(Z)), the relationship,

(A.9)
p
p
p
˘
˘
˘
where ξ1 , . . . , ξM are the elements of ξ̆ p: ξ̆ p = [ξ1 , . . . , ξ˘Mp ]T .
I
Equation (A.9) shows that denoting the jth row of XBI as β̆ j ,I
we get the relationship β̆ j ∈ KI . Thus, using the arguments inappendix A.2, we can prove that the relationship KI = rsp(XBI )holds.

z=x+yx ∈ rsp(X) andholds wherez ∈ rsp(X) + rsp(Y) holds, andrsp(X + Y) ⊂ rsp(X) + rsp(Y).
17y ∈ rsp(Y).
we have

Thus,
shown

Topical Review

J. Neural Eng. 14 (2017) 051001

An orthonormal set of basis vectors of the intersection
is obtained as a set of the principal vectors whose principal
angles are equal to zero. To find those principal vectors, wedefine matrices whose columns consist of the basis vectors
such that
U = xT1 , . . . , xTµ ,(A.22)A.6. Basis vectors that span sum of subspaces

Next, we derive the basis vectors of X + Y . Sets of thebasis vectors of X and Y are respectively denoted by
SX = {x1 , . . . , xµ } and SY = {y1 , . . . , yν }. Let us assume thatx ∈ X , y ∈ Y , and z ∈ X + Y . We then have
z=x+y=µ
j=1

Therefore, the relationshipcj xj +

ν


dj yj .


V = yT1 , . . . , yTν .
(A.23)

j=1

The results of singular-value decomposition of a matrix UT V
are expressed ascos(θ1 ) · · ·0..  T T ,T..
U
V = Q  ...
(A.24).
. 
0. . . cos(θν )z ∈ span([x1 , . . . , xµ , y1 , . . . , yν ])(A.20)holds. Defining a set of the basis vecors of X + Y as SX+Y ,equation (A.20) indicates that SX+Y is equal to

SX+Y = {x1 , . . . , xµ , y1 , . . . , yν }= {x1 , . . . , xµ } ∪ {y1 , . . . , yν } = SX ∪ SY .
(A.21)
The equation above shows that the basis vectors of the sumof subspaces are obtained as the union of the basis vector setsthat span each of the subspaces.

where Q and T are matrices whose columns consist of singular
vectors, and we assume that µ > ν . In equation (A.24), singular values of the matrix UT V are equal to the cosines of theprincipal angles between the two subspaces csp([xT1 , . . . , xTµ ])(= rsp(X)) and csp([yT1 , . . . , yTν ]) (= rsp(Y)). The intersection
has the property that the principal angles are equal to zero.
Thus, by observing the relation

A.7. Proof: If X ⊂ X  and Y ⊂ Y  , then (X ∩ Y) ⊂ (X  ∩ Y  )holdscos(θ1 ) = cos(θ2 ) = · · · = cos(θr )≈ 1 > cos(θr+1 )  · · ·  cos(θν ),If x ∈ X ∩ Y , the relationships, x ∈ X and x ∈ Y , hold. Since
X ⊂ X  and Y ⊂ Y  , we have x ∈ X  and x ∈ Y . Therefore,x ∈ X  ∩ Y  , and we have shown (X ∩ Y) ⊂ (X  ∩ Y  ).

the dimension of csp(U) ∩ csp(V), (namely, the dimension ofrsp(X) ∩ rsp(Y)) is determined to be r. The principal vectors
are then obtained either as the first r columns of the matrix
UQ or the first r columns of the matrix VT . Defining the firstr columns of UQ as ψ T1 , . . . , ψ Tr , the vectors ψ 1 , . . . , ψ r form
an orthonormal basis set for the intersection rsp(X) ∩ rsp(Y).

A.8. Proof of span(SX ∩ SY ) = span(SX ) ∩ span(SY )Two subspaces are denoted by X and Y , and sets of their basis
vectors by SX = {x1 , . . . , xµ } and SY = {y1 , . . . , yν }. Let us
also define a set of basis vectors SX ∩ SY = {z1 , . . . , zr } where
the dimension of the intersection is r. If x ∈ span(SX ∩ SY ),r
we can write x as a linear sum: x = j=1 wj zj . The factthat zj ∈ SX leads to x ∈ X = span(SX ). Also, the fact thatzj ∈ SY leads to x ∈ Y = span(SY ). Therefore, x ∈ X ∩ Y =span(SX ) ∩ span(SY ) holds, and we can show span(X ∩ Y) ⊂span(X ) ∩ span(X ).
The other direction of the proof is as follows. If
x ∈ span(SX ) ∩ span(SY ), we have x ∈ span(SX ) andx ∈ span(SY ). Thus, x is expressed as a linear sum ofbasis vectors that belong to both SX and SY , namely, a linear sum of {z1 , . . . , zr }. Therefore, x ∈ span(SX ∩ SY ), andthis indicates span(SX ) ∩ span(SY ) ⊂ span(SX ∩ SY ). Since
we have shown the both directions, we have shown
span(SX ∩ SY ) = span(SX ) ∩ span(SY ).

Appendix B. Signal space separation (SSS) method
B.1. Derivation of SSS basis vectors

This section presents an overview of the signal space separation (SSS) method [12–14] with an emphasis on the derivation of SSS signal and interference extractors. One ofbasic assumptions of the SSS method is that the sensors are
installed in a source-free region, which is referred to as thesensor region in this paper. Then, the magnetic field at r, B(r)is expressed using the spherical polar coordinate r = (r, θ, φ)as
B(r) = − µ0

A.9. Derivation of basis vectors that span intersectionof two row spaces

− µ0

∞ 



=1 m=−∞ =1 m=−α,mν ,m (θ, φ)r+2β,m r−1 ω ,m (θ, φ),

(B.1)


where µ0 indicates the magnetic permeability of free space.
In equation (B.1), ν ,m (θ, φ) and ω ,m (θ, φ) are the modified
vector spherical harmonics [12, 23].
In the right-hand-side of equation (B.1), the first term represents the magnetic field generated from sources locatedLet us assume that X and Y are low-rank data matrices. We
define the basis vectors of rsp(X) as SX = {x1 , . . . , xµ } where
μ is the dimension of rsp(X), and the basis vectors of rsp(Y)as SY = {y1 , . . . , yν } where ν is the dimension of rsp(Y). The
procedure used to find a set of basis vectors of rsp(X) ∩ rsp(Y)is described below. The procedure is according to [22].
18Topical Review

J. Neural Eng. 14 (2017) 051001

closer to the origin than the sensors. The region closer to theorigin than the sensors is referred to as the internal region.
The second term represents the magnetic field from sources
located farther from the origin than the sensors. The region
farther from the origin than the sensors is referred to as theexternal region. If an appropriate choice of the origin can
result in the signal sources of interest to be located within
the internal region and all interference sources to be locatedwithin the external region, equation (B.1) provides a natural
separation between the signal and interference. This is the key
idea of the SSS method.
Let us derive the SSS basis vectors. To do so, we denote theoutput of the jth sensor by yj, its location by rj and its orientation by ζ j . Then, we have

C = [c1,−1 , c1,0 , c1,1 , . . . , cLC ,LC ] ,(B.9)
D = [d1,−1 , d1,0 , d1,1 , . . . , dLD ,LD ] ,(B.10)T

α = [α1,−1 , α1,0 , α1,1 , . . . , αLC ,LC ] ,(B.11)T

φ = [β1,−1 , β1,0 , β1,1 , . . . , βLD ,LD ] ,(B.12)

we obtain αy = Cα + Dφ = [C, D]= Sx
(B.13)φ

where S = [C, D] and x = [αT , φT ]T . Here, C is an M × NC
matrix, and D is an M × ND matrix where2
NC = LC2 + 2LC , and ND = LD
+ 2LD .
(B.14)jjyj = B(rj ) · ζ j = yint
+ yext,(B.2)

where the notation · represents the inner product between twojjvectors, and yint
and yext represent the two components of thejth sensor outputs generated from the internal and external
regions, respectively. These components are expressed as



∞ 

ν ,m (θj , φj ) · ζ jjyint = −α,m,(B.3)rj+2
=1 m=−

B.2. SSS signal and interference extractors

When the key assumption that sources of interest are locatedin the internal region and all interference sources are locatedin the external region holds, yint represents the signal ofinterest and yext represents the interference. Assuming that
M > NC + ND, and using equation (B.13), the least squares
estimate of x is obtained as:


∞ 



j
yext = −β,m rj−1 ω ,m (θj , φj ) · ζ j ,(B.4)
=1 m=−
x = (ST S)−1 ST y.
(B.15)T

where we set µ0 = 1 for simplicity. Let us define theinternal and external components of the data vector y asTT
and yext = [y1ext , . . . , yM
yint = [y1int , . . . , yM
ext ] , which areint ]expressed such thatyint = Cα.
(B.16)

We now derive the SSS signal and interference extractors, and rewrite equation (B.16) using the signal extractor.
To do so, let us define an operation to make a new column
vector [ai , . . . , aj ]T by using the ith to jth components of
a = [a1 , . . . , aM ]T as [a][i:j] (namely, [a][i:j] = [ai , . . . , aj ]T ).
From equation (B.15), we have
 = (ST S)−1 ST y [1:N ] .
α(B.17)C∞ yint =
α,m c,m ,(B.5)
=1 m=−∞ yext =
β,m d,m ,(B.6)
=1 m=−

where column vectors c,m and d,m are given by 1
[ν ,m (θ1 , φ1 ) · ζ 1 ]
r1+2




..
c,m = 
 and
.


1
+2 [ν ,m (θM , φM ) · ζ M ]
rM

 −1
r1 [ω ,m (θ1 , φ1 ) · ζ 1 ]


..
.
d,m = 
(B.7)
.


−1
rM [ω ,m (θM , φM ) · ζ M ]
Truncating the summation with respect to the index  to LC
for yint and LD for yext, we finally obtainy = yint + yext =

LC 



=1 m=−α,m c,m +

LD 



T

 ]T and using α
 ,φx = [αConsidering that  , the signal comp­onent yint is estimated as

When the condition number of a matrix ST S is small, i.e. thecolumn spaces of C and D are well separated, the relationship = (ST S)−1 ST y [1:N ] ≈ (ST S + κI)−1 ST y [1:N ]α(B.18)C
C

holds, where κ is a small positive constant10. Using the matrix
inversion formula(ST S + κI)−1 ST y = ST (SST + κI)−1 y,(B.19)we get = (ST S + κI)−1 ST y [1:N ] = ST (SST + κI)−1 y [1:N ]
α
C
C



= CT (SST + κI)−1 y ≈ CT (SST )−1 y = CT (CCT + DDT )−1 y.

(B.20)

β,m d,m .

=1 m=−10This equation simply claims that when the condition number of ST S issmall, nearly identical inverse matrices can be obtained either with or without the regularization term κI .

(B.8)
Thus, defining
19Topical Review

J. Neural Eng. 14 (2017) 051001

Using equation (B.16), we finally obtain[7] Adachi Y, Shimogawara M, Higuchi M, Haruta Y andOchiai M 2001 Reduction of non-periodic environmental
magnetic noise in meg measurement by continuously
adjusted least squares method IEEE Trans. Appl.
Supercond. 11 669–72[8] De Cheveign A and Simon J Z 2008 Sensor noise suppression
J. Neurosci. Methods 168 195–202[9] Watanabe T, Kawabata Y, Ukegawa D, Kawabata S, Adachi Y
and Sekihara K 2013 Removal of stimulus-induced artifacts
in functional spinal cord imaging 35th Annual Int. Conf.
of the IEEE Engineering in Medicine and Biology Society(IEEE) pp 3391–4[10] Taulu S and Simola J 2006 Spatiotemporal signal space
separation method for rejecting nearby interference in MEG
measurements Phys. Med. Biol. 51 1759–68[11] Sekihara K, Kawabata Y, Ushio S, Sumiya S, Kawabata S,
Adachi Y and Nagarajan S S 2016 Dual signal subspace
projection (DSSP): a novel algorithm for removing large
interference in biomagnetic measurements J. Neural Eng.
13 036007[12] Taulu S and Kajola M 2005 Presentation of electromagnetic
multichannel data: the signal space separation method
J. Appl. Phys. 97 124905[13] Taulu S, Simola J and Kajola M 2005 Applications of thesignal space separation method IEEE Trans. Signal Process.
53 3359–72[14] Taulu S, Kajola M and Simola J 2004 Suppression ofinterference and artifacts by the signal space separation
method Brain Topogr. 16 269–75[15] Sekihara K and Nagarajan S S 2008 Adaptive Spatial Filters
for Electromagnetic Brain Imaging (Berlin: Springer)[16] Scharf L L 1991 Statistical Signal Processing: Detection,
Estimation and Time Series Analysis (New York:
Addison-Wesley)[17] Nurminen J, Taulu S and Okada Y 2008 Effects of sensor
calibration, balancing and parametrization on the signal
space separation method Phys. Med. Biol. 53 1975[18] Hämäläinen M, Hari R, Ilmoniemi R J, Knuutila J andLounasmaa O V 1993 Magnetoencephalography-theory
instrumentation and applications to noninvasive studies of
the working human brain Rev. Mod. Phys. 65 413–97[19] Mosher J C, Hamalainen M S, Pantazis D, Hui, H B,
Burgess R C and Leahy R M 2009 Generalized sidelobe
canceller for magnetoencephalography arrays IEEE Int.
Symp. on Biomedical Imaging: from Nano to Macro (IEEE)
pp 149–52[20] Sarvas J 1987 Basic mathematical and electromagnetic
concepts of the biomagnetic inverse problem Phys. Med.
Biol. 32 11–22[21] Goldberg J L 1991 Matrix Theory with Applications(New York: McGraw-Hill)[22] Golub G H and Van Loan C F 2012 Matrix Computations
vol 3 (Baltimore, MD: Johns Hopkins University Press)[23] Hill E 1954 The theory of vector spherical harmonics
Am. J. Phys. 22 211–4

yint ≈ CCT (CCT + DDT )−1 y.
(B.21)

Accordingly, we can conclude that the interference removal
by SSS is achieved by multiplyingΓS = CCT (CCT + DDT )−1
(B.22)

with the data vector y. That is, the matrix ΓS plays a role asa filter that passes the signal of interest and blocks the interference. In this paper, ΓS is called the SSS signal extractor.
Note that, since relationships (ΓS )2 = ΓS and (ΓS )T = ΓS do
not hold, ΓS is not a projector. In exactly the same manner, we
can derive ≈ DDT (CCT + DDT )−1 y.
yext ≈ Dφ
(B.23)

Therefore, defining ΓI as

ΓI = DDT (CCT + DDT )−1,
(B.24)
ΓI plays a role as a filter that passes the interference but blocks
the signal of interest. This ΓI is called the SSS interference
extractor in this paper.

References[1] Paulraj A, Ottersten B, Roy R, Swindlehurst A, Xu G and
Kailath T 1993 Subspace methods for directions-ofarrival estimation Handbook of Statistics ed N K Bose andC R Rao (Amsterdam: Elsevier) pp 693–739[2] Uusitalo M and Ilmoniemi R 1997 Signal-space projection
method for separating MEG or EEG into components
Med. Biol. Eng. Comput. 35 135–40[3] Tesche C, Uusitalo M, Ilmoniemi R, Huotilainen M, Kajola M
and Salonen O 1995 Signal-space projections of MEG
data characterize both distributed and well-localizedneuronal sources Electroencephalogr. Clin. Neurophysiol.
95 189–200[4] Nolte G and Curio G 1999 The effect of artifact rejection bysignal-space projection on source localization accuracy inmeg measurements IEEE Trans. Biomed. Eng. 46 400–8[5] Sekihara K, Poeppel D, Marantz A and Miyashita Y
2000 Neuromagnetic inverse modeling: application of
eigenstructure-based approaches to extracting corticalactivities from meg data Image, Language, Brain: Papers
from the First Mind Articulation Project Symp. (Cambridge,
MA: MIT Press) p 197[6] Widrow B, Glover J R, McCool J M, Kaunitz J, Williams C S,
Hearn R H, Zeidler J R, Dong J E and Goodlin R C 1975
Adaptive noise cancelling: principles and applications
Proc. IEEE 63 1692–716

20