Journal of Neural Engineering
J. Neural Eng. 13 (2016) 036007 (19pp)

doi:10.1088/1741-2560/13/3/036007

Dual signal subspace projection (DSSP): a
novel algorithm for removing large
interference in biomagnetic measurements
Kensuke Sekihara1,5, Yuya Kawabata2, Shuta Ushio1, Satoshi Sumiya1,
Shigenori Kawabata1, Yoshiaki Adachi3 and Srikantan S Nagarajan4
1

Department of Advanced Technology in Medicine, Tokyo Medical and Dental University, 1-5-45
Yushima, Bunkyo-ku, Tokyo 113-8519, Japan2
Department of Systems Design and Engineering, Tokyo Metropolitan University, Asahigaoka 6-6, Hino,
Tokyo 191-0065, Japan
3
Applied Electronics Laboratory, Kanazawa Institute of Technology, Kanazawa-shi, 920-1331, Japan4
Biomagnetic Imaging Laboratory, University of California, San Francisco, 513 Parnassus Avenue, S362,
San Francisco, CA 94143, USA5
Signal Analysis Inc. Hachiouji, Tokyo 192-0031, Japan
E-mail: k-sekihara@nifty.com
Received 14 November 2015, revised 12 March 2016
Accepted for publication 18 March 2016
Published 11 April 2016
Abstract

Objective. In functional electrophysiological imaging, signals are often contaminated byinterference that can be of considerable magnitude compared to the signals of interest. This paper
proposes a novel algorithm for removing such interferences that does not require separate noise
measurements. Approach. The algorithm is based on a dual deﬁnition of the signal subspace inthe spatial- and time-domains. Since the algorithm makes use of this duality, it is named the dual
signal subspace projection (DSSP). The DSSP algorithm ﬁrst projects the columns of themeasured data matrix onto the inside and outside of the spatial-domain signal subspace, creating
a set of two preprocessed data matrices. The intersection of the row spans of these two matrices
is estimated as the time-domain interference subspace. The original data matrix is projected onto
the subspace that is orthogonal to this interference subspace. Main results. The DSSP algorithm
is validated by using the computer simulation, and using two sets of real biomagnetic data: spinal
cord evoked ﬁeld data measured from a healthy volunteer and magnetoencephalography datafrom a patient with a vagus nerve stimulator. Signiﬁcance. The proposed DSSP algorithm iseffective for removing overlapped interference in a wide variety of biomagnetic measurements.
Keywords: biomagnetism, interference removal, MEG, spinal cord evoked ﬁeld, signal
subspace, interference subspace(Some ﬁgures may appear in colour only in the online journal)
1. Introduction

developed to deal with such interference, these algorithms
often rely on availability of separate measurements that capture the statistical properties of the interference [1–3].
Therefore, if such separate measurements are not available,
those existing algorithms would not be effective for removing
overlapped interference. We present two examples of suchcases in this paper: imaging of the spinal cord electrophisiological activitiy and MEG imaging for epilepsy patients
implanted with vagus nerve stimulators.

Functional electrophysiological imaging can be achieved
using a variety of modalities including electroencephalography (EEG), magnetoencephalography (MEG) or
magnetospinography (MSG). In all of these modalities,
measurement signals are often contaminated by interference
that can be of considerable magnitude compared to the signalsof interest. Although quite a few algorithms have been
1741-2560/16/036007+19$33.00

1© 2016 IOP Publishing Ltd Printed in the UK

J. Neural Eng. 13 (2016) 036007

K Sekihara et al

There has been growing interest in dynamic source
imaging of spinal cord electrophysiological activity using its
evoked magnetic ﬁeld. In spite of the fact that spinal corddisorders are very common [4], there are no effective methods
for accurate diagnosis of spinal cord lesions, primarily
because spinal cord abnormalities found in a patient’s anatomical image (such as MRI or x-ray images) often do notcorrespond with the patient’s clinical symptoms. Thus,
functional imaging of spinal cord nerve activity can be a
promising diagnostic tool for spinal cord disorders.
Biomagnetometers optimized for measuring the spinal cord
evoked magnetic ﬁeld (SCEF) have been developed [5–7], alongwith efﬁcient source reconstruction algorithms suited to functional spinal cord imaging [8, 9]. However, one serious problem,
the removal of large stimulus-induced artifacts5, remains to besolved. Such artifacts exist for 8–10 ms immediately after the
stimulus onset. These artifacts overlie the SCEF signal, andcan distort imaging results of spinal cord activity as shown in
our experiments described in sections 3.1 and 4.
Similarly, in MEG measurements, extremely large artifacts sometimes contaminate the measured data and it can bequite difﬁcult to visualize the underlying signal of interest
within this artifact. One striking example is in MEG recordings obtained from epilepsy patients implanted with vagal
nerve stimulators (VNS) where the artifact from the stimulator and the lead-wires can completely contaminate therecordings such that it is extremely difﬁcult to see interictal
epileptiform activity or stimulus evoked responses frompatients’ primary sensory cortices.
This paper proposes a novel algorithm to remove these
artifacts from biomagnetic measurements. The algorithm isbased on the two kinds of deﬁnitions of the signal subspace, as
the spatial-domain signal subspace and the time domain signal
subspace. Since the algorithm makes use of this duality, it isnamed the dual signal subspace projection (DSSP) algorithm.
The DSSP algorithm ﬁrst projects the columns of the measured
data matrix onto the inside and outside of the spatial-domain
signal subspace, creating two ‘projected’ data matrices.
The intersection of the row spans of those two ‘projected’matrices is then taken to be an estimate of the time-domaininterference subspace, and artifact removal is carried out on
the basis of this estimated interference subspace. In this paper,
the theory behind the DSSP algorithm is provided in sufﬁcient details to explain why the interference subspace can be
estimated in such a manner. The proposed algorithm is ﬁrst
validated by using the computer generated data that simulates
the SCEF measurements and by using the SCEF data takenfrom a healthy volunteer.
It should be mentioned that estimating the interferencesubspace as the intersection between two time-domain subspaces was ﬁrst proposed in [10]. That method is called
temporal signal space separation (tSSS). However, there is a
key difference between tSSS and DSSP in that tSSS uses thevector spherical harmonics expansions, which are not used in

DSSP. The application of the tSSS method, therefore, islimited to MEG applications in which the data is acquired byan array of sensors that are arranged on a surface of a spherelike helmet. The tSSS method is, therefore, not applicable toour SCEF measurements, in which the sensors are arranged
on a nearly-ﬂat surface.
In contrast, the proposed DSSP algorithm is a versatilealgorithm. Its application is not limited to the artifact removal
in SCEF measurements, and it is applicable to a wide varietyof biomagnetic data. To demonstrate this versatility, we present the results of applying the algorithm to MEG dataacquired from a patient with a vagus nerve stimulator.
Although the MEG data contained a large amount of interference generated from the stimulator, primary somatosensory
cortices were successfully localized by interference removalusing the proposed DSSP algorithm.

2. Dual signal subspace projection
2.1. Data model

We assume that the biomagnetic measurements are conductedby using an array of M sensors. We denote the output of the mth
sensor at time t by bm (t ) and the column vector containing theoutputs of all sensors is b (t ): b (t ) = [b1 (t ), ¼, bM (t )]T . Theb (t ) is called the data vector. We use a data model in which themeasured data consists of the signal magnetic ﬁeld generated
from signal sources of interest, the interference magnetic ﬁeld
generated from interference sources, and the sensor noise. Thatis, the data vector is expressed asb (t ) = bs (t ) + bI (t ) + e ,(1 )where bs (t ), the signal vector, represents the signal magnetic
ﬁeld, bI (t ), the interference vector, represents the interference
magnetic ﬁeld, and e, the noise vector, represents the sensor
noise. The spatio-temporal data matrices of b (t ), bs (t ), andbI (t ) are, respectively, denoted by B, BS, and BI , such thatB = [b (t1), ¼, b (tK )] ,(2 )BS = [bs (t1), ¼, bs (tK )] ,(3 )BI = [bI (t1), ¼, bI (tK )] ,(4 )where the data acquisition takes place at t1, ¼, tK and M < K
is assumed.
The M×K spatio-temporal data matrix B is modeled as
B = BS + BI + Be ,(5 )where Be is a matrix whose columns consist of the noise
vectors e at t1, ¼, tK . The spatio-temporal data model inequation (5) is used in our analysis.
2.2. Sensor lead field and signal subspace

Let us assume that a unit-magnitude source exists at r . Whenthe source is oriented in the x-, y-, and z-directions, thecorresponding outputs of the jth sensor are, respectively,
denoted l jx (r ), l jy (r ), and l jz (r ). We deﬁne an M ´ 3 matrix
L (r ) in which the jth row is equal to the row vector

5

Although the exact causes of these artifacts are still unknown, we speculatethat they are caused by combined effects of stimulus-induced body electric
currents and transient responses of receiver electronics.
2

J. Neural Eng. 13 (2016) 036007

K Sekihara et al[l jx (r ), l jy (r ), l jz (r )]. This matrix L (r ), which is called the leadﬁeld matrix, represents the sensitivity of the whole sensor
array at the location of r .
Let us assume that, in total, Q sources exist. Theirlocations are denoted by r1, ¼, rQ , their orientations by
h1, ¼, hQ 6, and their intensities by s1 (t ), ¼, sQ (t ). Then, thesignal vector bs (t ) is expressed as
Q

Q

p=1

p=1

bs (t ) = å sp (t ) L (rp ) hp = å sp (t ) lp,

2.3. Time-domain signal and interference subspaces

Now we deﬁne the time-domain signal subspace, which isdifferent from the spatial domain signal subspace deﬁned in
section 2.2. For the arguments on the time-domain signal subspace, we denote row vectors consisting of the signal-source
time courses by sp : sp = [sp (t1), ¼, sp (tK )] ( p = 1, ¼, Q). We
assume that Q is smaller than the number of sensors M. This
assumption is called the low-rank signal assumption. Also, thevectors s1, ¼, sQ are assumed to be linearly independent. Using
equations (3) and (6), we have the following relation(6 )where lp = L (rp ) hp , and lp is called the lead ﬁeld vector of
the pth source. The equation above indicates that the signal
vector bs (t ) is expressed as a linear combination of the leadﬁeld vectors. In other words, the vector bs (t ) lies within the
span of the lead ﬁeld vectors, l1, ¼, lQ , i.e.
bs (t ) Î span {l1, ¼, lQ}.

⎤
⎡ Q
Q
BS = ⎢ å sp (t1) lp, ¼, å sp (tK ) lp⎥⎥⎦⎢⎣ p = 1
p=1
⎡ Q
⎤ ⎡ Q
⎤
⎢ å [sp (t1), ¼, sp (tK )] l1p ⎥ ⎢ å sp l1p ⎥⎢ p=1⎥ ⎢ p=1
⎥
⎥=⎢
⎥.
=⎢


⎢
⎥ ⎢
⎥
⎢ Q
⎢ Q
p⎥
p⎥
⎢ å [sp (t1), ¼, sp (tK )] lM ⎥ ⎢ å sp lM ⎥⎣ p=1
⎦ ⎣ p=1⎦(7)This subspace, span {l1, ¼, lQ}, is referred to as the (spatialdomain) signal subspace [11].
In general, the signal subspace is unknown because thesource locations and the source orientations are unknown.
However, a subspace that approximates the signal subspacecan be obtained from an augmented lead ﬁeld over the source
space, which is the region where a source may exist. To
derive the augmented lead ﬁeld, we deﬁne voxels over asource space. Denoting the locations of the voxels by
r1, ¼, rN , the augmented lead ﬁeld matrix over all voxel
locations is deﬁned as
F = [L (r1), ¼, L (rN )] .

Here, l1p, ¼, lMp are the elements of the lead ﬁeld vector
lp : lp = [l1p, ¼, lMp ].
The row vector whose elements form the jth row of thematrix BS is denoted b Sj . Equation (13) then implies therelationQ

b Sj = å l jp sp .

Ê

span {l1, ¼, lQ} .

(8 )Namely, b Sj is expressed as a linear combination of thevectors s1, ¼, sQ . The row vector b Sj thus lies within the spanof these vectors. That isb Sj Î span {s1, ¼, sQ }.

(9)

The column span of F is called the pseudo signal subspace inthis paper, and the projector onto the pseudo signal subspacecan be derived by applying eigenvalue decomposition to
FF T , such that7⎡ g1 0  0 ⎤ ⎡ T ⎤⎢⎥ e
0 g2  0 ⎥ ⎢ 1 ⎥FF T = [e1, ¼, eM ] ⎢⎢  ⎥.
⎢     ⎥⎢ T ⎥⎢⎣ 0 0 ¼ gM ⎥⎦ ⎣ eM ⎦

S = span {s1, ¼, sQ }.

(16)

It can be shown under the low-rank signal assumption that S
is equal to the row span of BS, i.e.,(10)S = span {b1S, ¼, b SM}.

(17)

The proof is presented in section A.1 in the appendix.
Let us deﬁne the time-domain interference subspace,which plays a key role in the proposed algorithm. We express
the interference vector asPbI (t ) = å sj (t ) xj ,

(18)

j=1

(11)

and the following relationship holds:
Pbs (t ) = bs (t ) .

(15)

This subspace, span {s1, ¼, sQ }, is deﬁned as the time-domain
signal subspace, and denoted S:

We assume that the eigenvalues g1, ¼, gz are distinctively
larger than the other eigenvalues. Then, the set of eigenvectors {e1, ¼, ez } is an orthonormal basis of the column span of
the augmented lead ﬁeld. Accordingly, the pseudo signal
subspace projector P is derived as
P = [e1, ¼, ez ][e1, ¼, ez ]T ,

(14)

p=1

It is clear that (ignoring the voxel discretization error) thecolumn span of F includes the signal subspace, that iscolumn span of F

(13)

where P is the number of interference sources. Their intensities
are denoted by sj (t ) ( j = 1, ¼, P ), and their lead ﬁeld vectors
by xj ( j = 1, ¼, P ). We deﬁne row vectors consisting of theinterference-source time courses as sp : sp = [sp (t1), ¼, sp (tK )]( p = 1, ¼, P ). The vectors s1, ¼, sP are assumed to be linearly
independent. The span of s1, ¼, sP is deﬁned to be the time
domain interference subspace, and is denoted by I , i.e.
I = span {s1, ¼, sP}. The row vector whose elements form

(12)

6

The orientation hj is a 3×1 column vector consisting of its x, y, and z
components.
7
The matrix FF T is called the Gram matrix.
3

J. Neural Eng. 13 (2016) 036007

K Sekihara et al

the jth row of the matrix BI is denoted by b Ij . Using the same
arguments used to derive equation (14), we can derive therelation

equation (22) implies that in = S

(19)

p=1

p
where x1p, ¼, x M
are the elements of the lead ﬁeld vector xp :pp
xp = [x1 , ¼, x M ]. It is clear that b Ij is expressed as a linear
combination of the vectors s1, ¼, sP and that the vector b Ij
lies within the interference subspace (i.e., b Ij Î I ). Under
the low-rank signal assumption (P + Q < M ), it can beshown that span {b1I, ¼, b IM} is also equal to I .

P

P

j=1

j=1

(I - P ) bI (t ) = å sj (t )(I - P ) xj = å sj (t ) x¯ j ,(27)

where x¯ j = (I - P ) xj . Using the same arguments employed
to derive equation (25) and denoting the jth row vector of theImatrix (I - P ) BI by b¯ j , we can show

2.4. Finding the time-domain interference subspacePIpb¯ j = å sp x¯ j ,In the proposed DSSP algorithm, the sensor measurements B
are ﬁrst projected onto the inside and the outside of thepseudo signal subspace. Let us deﬁne Bin and Bout as the
projections of B onto the inside and the outside of the pseudo
signal subspace
Bin = PB,

(20)

Bout = (I - P ) B,

(21)

Bin = PB = P (BS + BI + Be )= BS + PBI + PBe ,ppwhere x¯1 , ¼, x¯M are the elements of the modiﬁed lead ﬁeld
I
vector x¯ p. This equation indicates that the relation b¯ j Î I
holds, and the row span of (I - P )BI is equal to I .
Thus, denoting the jth row vector of Bout by b out
j and theby,equation(23)
impliesthat out
span of b1out, ¼, b outoutMis given by out = IÈ ¯ e,(29)¯ e is the subspace spanned by the row vectors of
where (I - P ) Be . Equations (26) and (29) imply that the intersectionof the two subspaces  in and out is the time-domaininterference subspace I :(22)

Bout = (I - P ) B = (I - P ) BI + (I - P ) Be .

(28)

p=1

where I is the identity matrix. Using equations (5) and (12),
we have(23)

Note here that, due to the ‘dull’ cut-off property of the
projector P , neither PBI nor (I - P ) BI becomes negligibly
small.
We now analyze the subspace spanned by the row vectors of Bin . We can write
8

P

(26)

where the span of the row vectors of PBe is denoted by e.
The proof of (26) is presented in section A.3 in the appendix.
We next analyze the subspace spanned by the row vectors of Bout . We can writeP

b Ij = å x pj sp,

È I È e,

I =  inÇ out ,(30)¯ e = Æ holds9, where Æ indicates the empty set. It isif e Ç ¯ e = Æ, because theactually straightforward to show e Ç relationshipP

PbI (t ) = å sj (t ) Pxj = å sj (t ) x j ,j=1

PBe ((I - P ) Be )T = PBe BeT (I - P ) = r 2P (I - P ) = 0,

(24)

j=1

holds. This means that the rows of PBe are orthogonal to therows of (I - P ) Be . In the equation above, we assume that
Be BeT = r 2I , where r 2 is the power of the sensor noise.

where xj = Pxj . The equation above indicates that although
the projector P modiﬁes the lead ﬁeld vectors, it neverchanges the time courses of the interference sources. There~I
fore, deﬁning bj to be the jth row vector of the matrix PBI ,we can showPp~I
bj = å xj sp,

2.5. Derivation of interference-subspace projector and removalof interference

In this subsection, we present an algorithm to derive an
orthonormal basis of I . We ﬁrst extract  in and out from
Bin and Bout , respectively by applying singular value
decomposition:

(25)

p=1

p
p
where x1 , ¼, xM are the elements of the modiﬁed lead ﬁeld
vector xp . This equation indicates that the row vector of thematrix PBI lies within the time-domain interference subspace,~Ii.e., bj Î I ( j = 1, ¼, M ). We denote a row vector of Bin byb inj and the span of b1in, ¼, b in
M by  in . Then, taking the~ISrelations b j Î S and bj Î I into consideration,⎡j1 0  0 ⎤ T⎢⎥ ⎡u ⎤
0 j2  0 ⎥ ⎢ 1 ⎥⎢
Bin = [ f1, ¼, fM ]⎢     ⎥⎢ T ⎥⎢⎢⎣ 0 0 ¼ j ⎥⎦ ⎣ uM ⎥⎦M(31)¯ e ¹ Æ,¯ e,
Ifdenoting
wehavee Ç 
e = e Ç  in Ç out = I È e and the intersection contains an additional noiserelated subspace.

89

An example of the cut-off property of P computed assuming the sourcesensor geometry used in our computer simulation is shown in section A.2 inthe appendix.
4

J. Neural Eng. 13 (2016) 036007

K Sekihara et al

Figure 1. (a) The source model used in the computer simulation. The model consists of four sources in which each pair of sources is 3 cm

apart. (b) Schematic drawing of the source-sensor geometry used in the computer simulation. The four sources, depicted by the four solid
arrows, travel along the y-axis with a speed of 60 cm s−1 in the y direction indicated by the white arrow. The source space is a twodimensional plane (-8  x  8 cm , -6  y  6 cm), which is located 7 cm below the sensor array.

andThus, by observing the relation⎡ f1 0  0 ⎤⎢⎥ ⎡v T⎤
0 f2  0 ⎥ ⎢ 1 ⎥⎢Bout = [g1, ¼, gM ] .
⎢     ⎥⎢ T ⎥⎢⎢ 0 0 ¼ f ⎥ ⎣ vM ⎥⎦⎣M⎦

cos (q1) = cos (q2) == cos (qr ) » 1,

the dimension of I is determined to be r.
The principal vectors are then obtained either as the ﬁrst rcolumns of the matrix UY or the ﬁrst r columns of the matrix
VZ . Deﬁning the ﬁrst r columns of UY as y1, ¼, yr , therelation(32)

I = span {y1, ¼, yr }For the two sets of singular values, distinctively large
singular values are denoted by j1, ¼, jm and f1, ¼, fn . We
can then derive the subspaces  in and out such that in » span {u1, ¼, um },(33)

 out » span {v1, ¼, vn }.

(34)

holds and the vectors y1, ¼, yr form an orthonormal basis forI . The interference removal is thus carried out by projecting
the measured data B onto the subspace orthogonal to theinterference subspace. Deﬁning a matrix G by
G = [y1, ¼, yr ] ,Additional arguments, presented in section 6, lead to theconclusion that the interference subspace I is equal to theintersection of span {u1, ¼, um } and span {v1, ¼, vn }. Thatis I = span {u1, ¼, um } Ç span {v1, ¼, vn }.
The procedure used to ﬁnd the intersection is described
below. According to [12], an orthonormal basis of the intersection is obtained as a set of the principal vectors whose
principal angles are equal to zero. To ﬁnd those principal
vectors, we ﬁrst deﬁne
U = [u1, ¼, um ] ,(35)

V = [v1, ¼, vn ] .

(36)

S = B (I - GG T) .
B

⎢
⎣


0

(40)

3. Computer simulation
3.1. SCEF measurements distorted by large artifacts

Computer simulations were carried out to investigate thevalidity of the proposed algorithm. We ﬁrst simulated SCEF
measurements distorted by large stimulus-induced artifacts.
For data generation, we assumed the sensor array of a 120channel biomagnetometer [5, 7], which has been speciﬁcallydeveloped for SCEF measurements. The biomagnetometer isequipped with 40 vector sensors10, which are arranged at
8×5 measurement locations covering a 14 cm×9 cm area.

⎡ cos (q1) 0 ⎤⎥ ⎥ ZT,¼ cos (qn )⎦(39)

the interference removal can be accomplished by rightmultiplying B by (I - GGT ). That is, the interferenceS are obtained as
removed measurements B

Singular-value decomposition of a matrix U TV is thenperformed, and the results are expressed as
U TV = Y ⎢

(38)

(37)

where Y and Z are matrices whose columns consist of
singular vectors. In the equation above, we use the fact thatm > n . Equation (37) indicates that the singular values of thematrix U TV are equal to the cosines of the principal angles
between span {u1, ¼, um } and span {v1, ¼, vn }. The intersection
I has the property that the principal angles are equal to zero.

10

A vector sensor can measure the radial component and two tangential
components of the magnetic ﬁeld. The details of the sensor are reported
in [13].
5

J. Neural Eng. 13 (2016) 036007

K Sekihara et al

Figure 2. (a) The sensor time courses of the signal magnetic ﬁeld, computed assuming that four sources travel along the y-axis with the speedof 60 cm s−1. (b) The sensor time courses of the artifact-only measurement, which was taken with the stimulus electrode positioned a few
centimeters away from the median nerve of a subject. (c) The sensor time courses of the combined (signal plus artifact) results. (d) The sensor
time courses of the artifact-removed results. The ordinates of these plots show relative values of ﬁeld intensity. (The values are common for
all four plots.)

The artifact data was obtained from artifact-only measurements, in which we applied exactly the same stimulus as
used in measuring the SCEF, with a stimulus electrodepositioned a few centimeters away from the median nerve of a
subject [17]. Since the location of the electrode was noton the nerve, stimuli did not induce the nerve activity but
elicited only the artifacts. The artifact-only data is shown inﬁgure 2(b).
The results of adding the artifact data onto the computergenerated signal are shown in ﬁgure 2(c), and the ﬁeld contour map of the artifact-contaminated data at a latency of 5 ms
is shown in ﬁgure 3(b). Here, the interference-to-signal ratio(ISR) was set at 12. The ISR is deﬁned as the ratio ofBI F BS F , where XF indicates the Frobenius norm of amatrix X . Both the sensor time courses and the contour map
show that a signiﬁcant amount of distortion arises due tocontamination by the artifacts.
We then applied the DSSP algorithm to this artifactcontaminated data. Here, a two-dimensional region(-8  x  8 cm , -6  y  6 cm), which is located 7 cm
below the sensor plane, was deﬁned as the source space. The
augmented lead ﬁeld was computed with a 0.5 cm voxel grid
over this source space. A plot of the singular values cos (qj ) inequation (37) is shown in ﬁgure 4. The plot shows that theﬁrst six singular values are very close to 1. Thus, thedimension of the intersection, r, was determined to be six bythresholding the singular values at 0.99. The artifact-removed
sensor time courses are shown in ﬁgure 2(d). A ﬁeld contour
map of the artifact-removed data at a latency of 5 ms is shown

In this computer simulation, we used a source model
consisting of four equal-intensity current vectors; such asource model has been found to be physiologically plausible
in our previous studies [8]. The model is shown in ﬁgure 1(a).
In this model, the two sources aligned along the y-axis areintercellular sources propagating in the nerve axon, representing two anti-directional current dipoles [14], called theleading and the trailing dipoles [15, 16]. The other twosources aligned along the x-axis represent the volume current.
The volume current ﬂows from the extracellular milieu to thesite between the two dipole sources in the nerve axon. In thiscomputer simulation, the distance between the two dipole
sources, as well as that between the two volume current
sources, was set to 3 cm.
The SCEF sources move with a speed of 40–80 cm s−1.
In this computer simulation, we assumed that four sources
traveled along the y-direction (from -y to +y ) with a speed of
60 cm s−1. The source-sensor geometry used in the computer
simulation is shown in ﬁgure 1(b). The magnetic ﬁeld generated from the four moving sources was computed, and aplot of the simulated measurements of all 120sensors isshown in ﬁgure 2(a). Here, we used no conductor models, andthe forward solution was computed using the well-known
Biot–Savart law, which is derived based on the quasi-static
approximation of Maxwell’s equations. The time was set atzero when the y coordinate of the leading dipole was equal to−30 cm. The ﬁeld contour map at a latency of 5 ms is shownin ﬁgure 3(a).
6

J. Neural Eng. 13 (2016) 036007

K Sekihara et al

Bx

By

-5

0
x (cm)

-5

5

-5

0
x (cm)y (cm)0y (cm)y (cm)

-5

5

(b)

Bx

-5

y (cm)

0x (cm)

5

0

-5

-5

0
x (cm)

5

Bz

0

-5

-5

0
x (cm)

5

By

5

0

0

positive-5-50x (cm)

5-50x (cm)

5Bz

5y (cm)

-5

5

5

-5

(c)

-5

0

5

y (cm)

5

0

Bz

5

(a)

0
x (cm)By

5y (cm)y (cm)y (cm)

0

-5

Bx

5

5

0

negative-5-50x (cm)

5

Figure 3. The ﬁeld contour maps at the latency of 5 ms for (a) the signal-only sensor time courses, (b) the combined (signal plus artifact)

sensor time courses, and (c) the artifact-removed sensor time courses. The upper-left, upper-right, and the lower-right panels in each ﬁgureshow, respectively, the plots of the x, y, and z components of the magnetic ﬁeld (denoted Bx, By, and Bz). (Note that the x and y components
are the tangential components. The z component is the normal component. ) The color of the contours indicates the relative intensity of themagnetic ﬁeld according to the color bar. The small ﬁlled circles indicate the 40 locations of the sensors.

in ﬁgure 3(c). Comparing between ﬁgures 2 and 3, we can seethat the distortion from the artifacts is signiﬁcantly reduced.
We next assessed the effectiveness of the DSSP algorithm on the basis of source reconstruction results. Time-pointby time-point source reconstruction is essential for imaging
spinal cord activity because of the rapid movements of thesources. We therefore used the recursively applied nullsteering (RENS) beamformer algorithm [9, 18] because of its
applicability to single time-point data. A snapshot of thesource image at 5 ms latency is shown in ﬁgure 5. In thisﬁgure, the results in (a) are obtained by using the artifact-free
sensor data in ﬁgure 3(a). The results from the artifact-contaminated sensor data in ﬁgure 3(b) are shown in ﬁgure 5(b).
Due to the presence of the interference, the reconstruction
results are signiﬁcantly distorted. The results from the artifactremoved sensor data in ﬁgure 3(c) are shown in ﬁgure 5(c). In

Figure 4. The singular values cos (qj ) in equation (37) versus theirorder j.

7

J. Neural Eng. 13 (2016) 036007

K Sekihara et al

Figure 5. The source reconstruction results obtained using the ﬁeld maps at 5 ms, which are shown in ﬁgure 3. (a) The source image from thesignal-only data in ﬁgure 3(a). (b) The source image from the signal-plus-artifact data in ﬁgure 3(b). (c) The source image from the artifactremoved data in ﬁgure 3(c).

We put a single interference source shown as a ﬁlled
circle with the label ‘interference (location A)’ in ﬁgure 6(a).
As shown in this ﬁgure, the interference source was fairly
close to the signal sources and it was located just outside thehead near the head surface. The coordinate of the interference
source was (-1, -1, -6) cm. This interference simulated thenoise caused from some types of brain stimulator. In
section 5, we actually present experiments in which the MEG
data was contaminated by the noise from patient’s VNS. The
time course of the interference source is presented in thebottom panel of ﬁgure 6(b). We generated interference
magnetic ﬁeld with the ISR equal to 100 where the ISR isdeﬁned as BI  BS . The interference magnetic ﬁeld wascomputed and overlapped onto the signal magnetic ﬁeld
computed from the activities of the three signal sources.
To generate the magnetic ﬁelds, source activities were
projected to the sensor time courses through the lead ﬁeld,
which is obtained using the homogeneous spherical headmodel [19] with the center of the sphere set to (0, 0, 4) cm.
The spatio-temporal data with 1200 time points was generated. In ﬁgure 7, the time courses of the signal magnetic ﬁeld

these results, the four sources are clearly reconstructed andthe results are very close to the artifact-free case in ﬁgure 5(a),demonstrating the effectiveness of the proposed method.

3.2. MEG measurement overlapped with large interference
magnetic field

We next present a computer simulation on MEG measurements in which the data was overlapped with a large interference caused near the source space. A sensor alignment of
the 275-channel whole-head sensor array from the OmegaTM(VMS Medtech, Coquitlam, Canada) neuromagnetometer wasused. The coordinate system and source-sensor conﬁguration
used in the computer simulation are depicted in ﬁgure 6(a). A
vertical plane (x=0 cm) was assumed at the middle of thesensor array, and three sources were assumed to exist on thisplane. The coordinates of the sources are (0, -2, 10.3),(0, 2.5, 10.3), and (0, 1, 7.3) cm. The time courses assigned
to the three sources are shown in the top three panels inﬁgure 6(b).
8J. Neural Eng. 13 (2016) 036007

K Sekihara et al

Figure 6. (a) The coordinate system and source-sensor conﬁguration used in the computer simulation in section 3.2. Small three circles show
the locations of the signal sources. The ﬁlled circles with the labels of‘interference (location A)’ and ‘interference (location B)’ indicate thetwo locations of the interference source. The square indicates the y–z cross section of the source space. The large circle indicates the sphereused for the forward calculation. (b) The assumed time courses of the ﬁrst, the second, and the third sources are shown respectively in theﬁrst, the second and the third panels from the top. The bottom panel shows the assumed time course of the interference source. The ordinate
indicates the normalized, relative intensity and the abscissa indicates the time points.

100 times stronger than the signal magnetic ﬁeld, the sensor
time courses in the bottom panel are dominated by theinterference magnetic ﬁeld.
We set the source space to a region of -4.5
x  4.5 cm , -5  y  5 cm and 5  z  13 cm . Thisregion approximately covers the whole brain and includes thelocations of the three sources, but it does not include thelocation of the interference source. The augmented lead ﬁeld
was computed using a 0.5 cm voxel grid over this region. We
then applied the DSSP algorithm for removing the interference. The dimension of the interference subspace r wasdetermined to be 1 by thresholding the cosine of the principal
angle (equation (37)) by the value of 0.99. The interferenceremoved sensor time courses are shown in the bottom panel inﬁgure 7, which shows that the interference is nearly completely removed.
We next evaluate the effect of interference removal basedon the quality of reconstructed source images. The adaptive
beamformer algorithm [11] was applied for source reconstruction from the sensor data shown in ﬁgure 7. The results areshown in ﬁgure 8. In ﬁgure 8(a), the source reconstruction
results obtained from the interference-free magnetic ﬁeld (thetop panel of ﬁgure 7) are shown. Three sources can clearly beobserved. The reconstruction results obtained using the interference-overlapped magnetic ﬁeld (the middle panel of ﬁgure 7)
are shown in ﬁgure 8(b). Due to the overlap of the large
interference, the results contain signiﬁcant distortion. The
reconstruction results obtained using the artifact-removed sensor
data (the bottom panel of ﬁgure 7) are shown in ﬁgure 8(c) inwhich the distortion is nearly completely removed, demonstrating the effectiveness of the proposed method.

Figure 7. The sensor data corresponding to the signal magnetic ﬁeld

plus sensor noise (top). The interference-added sensor data (middle).
The interference-to-signal ratio (ISR) was set equal to 100 in which
the ISR is deﬁned by the ratio BI  BS . (The time courses of thesensor data in this panel is dominated by the interference-source timecourse, due to the large ISR.) The interference-removed sensor timecourses (bottom).

(plus sensor noise) are shown in the top panel. The timecourses of interference-overlapped magnetic ﬁeld are shownin the middle panel. Since the interference magnetic ﬁeld is
9

J. Neural Eng. 13 (2016) 036007

K Sekihara et al

Figure 8. Results of source reconstruction obtained using the sensor data in ﬁgure 7. (a) Results obtained from the interference-free magnetic

ﬁeld (the top panel of ﬁgure 7). (b) Results using the interference-overlapped magnetic ﬁeld (the middle panel of ﬁgure 7). (c) Results
obtained using the interference-removed sensor data (the bottom panel of ﬁgure 7). The reconstruction region was set equal to the source
space -4.5  x  4.5 cm , -5  y  5 cm , and 5  z  13 cm . The source distribution on the plane x=0 cm is shown. The ﬁlled circle
with the label‘interference (location A)’ in ﬁgure 6(a) indicates the location of the interference source in this simulation.

Figure 9. Results of source reconstruction when the interference source put at (−1, −1, 5), which is just on the border of the source space. (a)Results obtained without the interference removal. (b) Results obtained with the DSSP interference removal. The reconstruction region wasset equal to the source space -4.5  x  4.5 cm , -5  y  5 cm , and 5  z  13 cm . The source distribution on the plane x=0 cm isshown. The ﬁlled circle with the label‘interference (location B)’ in ﬁgure 6(a) indicates the location of the interference source in thissimulation.

(−1, −1, 5), exactly on the border of the source space. Thislocation is indicated by a ﬁlled circle with the label ‘interference (location B)’ in ﬁgure 6(a). The results of sourcereconstruction are shown in ﬁgure 9 where the results without
interference removal are shown in (a) and the results with theremoval in (b). These results clearly indicate that the DSSPAs shown in the computer simulation above, the DSSP
algorithm allows interference sources to be fairly close to thesource space. However the prerequisite exists that the interference sources are located outside the source space. We
performed computer simulation in which this prerequisite was
not fulﬁlled. We here put a single interference source at
10

J. Neural Eng. 13 (2016) 036007

K Sekihara et al

Figure 10. Schematic illustration of the experimental setup formeasuring the spinal cord evoked ﬁeld. The subject lies down in thesupine position with the lower neck placed on the upper surface ofthe protrusion of the cryostat; the protrusion contains 40 vector
sensors directed upward, covering a 14 ´ 9 cm area. The sensors arearranged at 8×5 measurement locations. A stimulus current isapplied to the subject’s median nerve near the elbow.
Figure 11. X-ray image covering the subject neck and the sensors,

used for determining a plane containing the spinal cord. This 2D
plane is deﬁned as the source space and the source reconstruction is
performed on this plane.

algorithm can still remove the inﬂuence of interference but theremoval is achieved at the sacriﬁce of signiﬁcant amount ofsignal source intensity.

4. Experiments using SCEF measurements from ahealthy volunteer
4.1. SCEF measurements

A 120-channel biomagnetometer [5, 7] was used for measuring a human SCEF11. The subject was a healthy malevolunteer. The experiment was approved by the ethics committee of Tokyo Medical and Dental University.
The experimental setup is schematically shown inﬁgure 10. As depicted here, the cryostat of the biomagnetometer has a cylindrical body with a protrusion, and thisprotrusion contains sensors directed upward. The subject lies
down in the supine position, and the subject’s lower neck is
positioned on the upper surface of the protrusion of thecryostat. A stimulus current was applied to the subject’s
median nerve near his left elbow. The stimulus with anintensity of 10mA and a duration of 0.3 ms was repeated 2000
times at a repetition rate of 4 Hz. The data acquisition was
performed with a sampling frequency of 40 kHz. An analog
bandpass ﬁlter with a bandwidth of 100–5000 Hz was applied.
The signal was averaged across all 2000 measured trials.
An x-ray image covering the subject’s neck and thesensors of the biomagnetometer was obtained to identify thelocation of the spinal cord. The reconstruction region wasdetermined to be a curved plane containing the spinal cord.
An x-ray image with the extracted 2D reconstruction region isshown in ﬁgure 11.
The measured SCEF, averaged over 2000 trials, is shownin ﬁgure 12(a). The electric stimulus was given at a latency of
0 ms. Large artifacts are observed particularly in the data

Figure 12. (a) The original sensor time courses. The results averaged
across 2000 trials are shown. (b) Artifact-removed sensor timecourses. The reconstruction plane shown in ﬁgure 11 was used as thesource space over which the augmented lead ﬁeld matrix wascomputed. The dimension of the intersection, r, was set at three,which was determined by applying the 0.99 thresholding criterion.
The fairly sharp peak near 2.5 ms is due to a hardware problem, andshould be disregarded.

before 4 ms, although the peaks between 4 and 8 ms, which
are caused by the spinal cord nerve activity, can still beobserved. We applied the proposed DSSP algorithm toremove these artifacts. The 2D reconstruction region was usedas the source space over which the augmented lead ﬁeld
matrix was computed. The results are shown in ﬁgure 12(b).

11The sensor array of this biomagnetometer was used for the data
generation in our computer simulation.

11

J. Neural Eng. 13 (2016) 036007

K Sekihara et al

When the subject’s left median nerve is stimulated, thenerve activity is known to propagate from the left median
nerve into the spinal cord near the fourth vertebra (c4). Thus,
it should be more plausible that the current vector has atransverse, negative x component, as depicted in ﬁgure 13.
Accordingly, considering the fact that the leading dipoles in(b) and (c) have such negative x components, we can draw theconclusion that the results in ﬁgure 14(b) are physiologicallymore plausible than those in ﬁgure 14(a).
We have conducted measurements exactly the same as
the ones described above but with right median nerve stimulation. The source images from these measurements areshown in ﬁgure 15. Here, the source image from the original
artifact-contaminated sensor data is shown in (a). The image
from the artifact-removed data is shown in (b). Both results
show the leading dipole but the direction of the leading dipole(indicated by the white arrow) is again different between these
two images. The source image obtained with the stimulus
applied at the right median nerve near subject’s wrist is shownin (c). In this ﬁgure, the current vector should have a transverse, positive x component since the x component of thesource should be opposite to that in the left median nerve
stimulation. The images in ﬁgures 15(b) and (c) show that thedipoles have the expected positive x components. Therefore,
we can again conclude that the results in (b) are physiologically more plausible than those in (a), which show the source
vector with an almost upward direction. It should be notedthat the DSSP algorithm can be used with any type of sourcereconstruction methods. In these experiments, the RENS
beamformer was used for reconstructing the spinal cord
activity, primarily because it gives the spatial resolution better
than other existing methods12 applicable to single time
point data.

Figure 13. Schematic illustration of the reconstruction region, whichis indicated by the square. Relative positions of the subject’s neckand median nerve are shown. The arrow indicates a typicalorientation of the current vector of the leading dipole ﬂowing fromthe peripheral nerve into the spinal cord when subject’s left median
nerve is stimulated.

Here, the artifacts are signiﬁcantly reduced and the SCEF
signal, consisting of peaks between 4 and 8 ms, is clearlyobserved.

4.2. Validation based on reconstruction of spinal cord nerve
activity

The spinal cord source activity was reconstructed by using the
RENS beamformer. Schematic illustration in ﬁgure 13 shows
relative positions of the subject’s neck and median nerve withrespect to the reconstruction region, which had an area of
16 ´ 12 cm with voxel dimensions of 0.5 cm in the x- andy-directions. The reconstructed source images at a latency of
5.8 ms are shown in ﬁgure 14. In this ﬁgure, the source image
from the original artifact-contaminated sensor data is shownin (a). The image from the artifact-removed data is shown in(b). Both images show the leading dipoles but their directions(indicated by the white arrows) are signiﬁcantly different.
To determine which results are physiologically moreplausible, we used a source image obtained with the stimulus
applied at the median nerve near the subject’s wrist; the imageis shown in ﬁgure 14(c). The signal obtained with wrist stimulation is known to be less affected by artifacts. This isbecause, with the stimulation near the subject’s wrist, it takes
4–5 ms more for the nerve activation to reach the neck region,and thus the spinal cord signal is much less contaminated bystimulus-induced artifacts due to their rapid decay. Therefore,
the results from the wrist stimulation should serve as thephysiological ‘ground truth’ for the comparison. The leading
dipoles have transverse (the negative x) components inﬁgures 14(b) and (c), while the direction is almost upward in(a), i.e., the current vector has almost no x component in (a).

5. Experiments using MEG data measured from apatient with a vagus nerve stimulator
The application of the DSSP algorithm is not limited toartifact removal in SCEF measurements, and the algorithm
can be applied to a wide variety of biomagnetic measurements. To demonstrate the versatility of the algorithm, we
present results from applying the DSSP algorithm to MEG
data obtained from a patient with a VNS. Generally, MEG
data from VNS patients contains a large amount of interference generated from the stimulator located near thepatient’s chest area. The measurements were conducted using
the 275-channel whole-head sensor array of the OmegaTM(VMS Medtech, Coquitlam, Canada) neuromagnetometer.
In the upper panel of ﬁgure 16(a), the somatosensory
MEG data measured with tactile stimulation applied to thepatient’s left index ﬁnger is shown. Since the tactile stimulation was used, a large peak should be observed around thelatency of 50 ms. However, such a peak is not observed due tothe presence of interference from the VNS stimulator.
12

Such existing methods include the minimum-norm [20, 21] andsLORETA methods [22].

12J. Neural Eng. 13 (2016) 036007

K Sekihara et al

Figure 14. Results of source reconstruction when a stimulus applied to the subject’s left median nerve. (a) Source image at the latency of 5.8

ms obtained from the original sensor data in ﬁgure 12(a). (b) Source image at the latency of 5.8 ms from the artifact-removed data inﬁgure 12(b). (c) Source image at 9.85 ms obtained with the stimulus applied to the median nerve near the subject’s wrist. (This image wasobtained without using DSSP artifact removal.) The white arrows indicate the directions of the leading dipoles. The relative intensity of thereconstructed source is color-coded according to the color bar. The sketch of the spine was drawn from the overlaid x-ray image used foraligning the sensor coordinates to the subject’s neck position.

of the primary somatosensory cortices, demonstrating that the
DSSP algorithm can be used for the interference removal in
MEG. It should be mentioned that, without interference
removal, a single strong fake source was detected near thecenter of the head for the both cases of left- and right-index
ﬁnger stimulations, although these results are not presented inthis paper.

The sensor time course with the interference removed by
the proposed DSSP method is shown in the lower panel ofﬁgure 16(a). Here, a peak around the latency of 50 ms isclearly observed. The MEG data taken from the same patientwith the tactile stimulation applied to the patient’s right index
ﬁnger is shown in ﬁgure 16(b). The upper panel shows theoriginal sensor time courses and the lower panel shows thesensor time courses processed with the DSSP algorithm.
Again, the sensor time courses processed by the DSSP
algorithm show a peak around the latency of 50 ms, but sucha peak is not observed in the original sensor time courses.
The source localization results are shown in ﬁgure 17.
The results for the case of the left index ﬁnger stimulation are
shown in ﬁgure 17(a), and those for the right index ﬁnger
stimulation are shown in ﬁgure 17(b). In both cases, thesource activity is successfully localized near the primarysomatosensory area in the contra-lateral hemisphere. These
reconstruction results show that the DSSP algorithm removes
the inﬂuence of the VNS stimulator and enables the mapping

6. Discussion and summary
Since equations (33) and (34) are important in the derivationof the proposed method, we should further clarify themeaning of these equations. Theoretically, if Q and P are
known and μ and ν are respectively set to Q+P and P, thefollowing relationships hold:S

13È I = span {u1, ¼, uQ+P},(41)

I = span {v1, ¼, vP }.

(42)

J. Neural Eng. 13 (2016) 036007

K Sekihara et al

Figure 15. Results of source reconstruction when a stimulus applied to the subject’s right median nerve. (a) Source image at the latency of 5.8

ms obtained from the original sensor data. (b) Source image at the latency of 5.8 ms from the artifact-removed data. (c) Source image at
9.85 ms obtained with the stimulus applied to the median nerve near the subject’s wrist. See ﬁgure 14 for legend.

The proof can be found, for example13, in [23] or [11].
Therefore, we could directly derive the interference subspace
I using equation (42) if we knew the number of interference
sources P.
However, information on the number of sources is difﬁcult to obtain in real-life applications. In addition, the singular values in equations (31) and (32) may not have sharp
cutoffs due to a low signal-to-noise ratio condition or an
insufﬁcient number of time points. Accordingly, there isalways considerable ambiguity in determining μ and ν, thethreshold values of the singular value spectra. A general
strategy for such cases is to overestimate μ and ν. However,span {u1, ¼, um } and span {v1, ¼, vn } then unavoidably containa part of the noise subspace. Thus, to be exact with overestimated μ and ν, equations (33) and (34) should be rewrittenasSÈ I È ¢e = span {u1, ¼, um },IÈ ¯ ¢e = span {v1, ¼, vn },(44)

where ¢e is the part of e that is contained in span {u1, ¼, um }¯ e that is contained in span {v1, ¼, vn }.
¯ ¢e the part of and The arguments above indicate that the choices of theparameters μ and ν are not crucial for implementing thealgorithm, because the excess dimensions of span {u1, ¼, um }and span {v1, ¼, vn } do not affect the extraction of the interference subspace, since the intersection between
span {u1, ¼, um } and span {v1, ¼, vn } is computed. In all theexperiments described in previous sections, the parameters μand ν were always set at 20. In the computer simulation insection 3.1, the ﬁeld map and the reconstructed source image
were recomputed with μ and ν set to 40. The results areshown in ﬁgures 18(a) and (b), which are almost the same as
the results obtained with μ and ν set at 20 (shown inﬁgures 3(c) and 5(c)).
Unlike μ and ν, the choice of the dimension of theinterference subspace, r, can be crucial in the implementationof the algorithm. The inﬂuence of overestimating r is analyzed as follows. Even when the parameter r is overestimated,

(43)

13
Equation (41) can be proved, for example, by replacing B in section13.1
of [11] with BinT and by following the same arguments in that section.

14

J. Neural Eng. 13 (2016) 036007

K Sekihara et al

Figure 16. Somatosensory MEG data obtained from a patient with a vagus nerve stimulator (VNS). (a) Sensor time courses when the tactile

stimulus was applied to the subject’s left index ﬁnger. (b) Sensor time courses when the tactile stimulus was applied to the subject’s rightindex ﬁnger. Original sensor time courses are shown in the upper panels, and sensor time courses processed by the proposed DSSP algorithm
for interference removal are shown in the lower panels.

since the interference time-course vectors sTp ( p = 1, ¼, P )still lie within the column span of G , the orthogonal projector(I - GGT ) suppresses the interference BI . Therefore, only
when the projector (I - GGT ) also suppresses the signal ofinterest, the overestimation becomes problematic. This scenario happens if the excess dimensions of the interferencesubspace include a part of the signal subspace.
Part of the signal subspace may be included in the estimated interference subspace when the distance between thesignal subspace, span {s1, ¼, sQ }, and the interference subspace, span {s1, ¼, sP}, is small, i.e., when the two subspaces
are close. The distance between the two subspaces can bedeﬁned using the subspace angles between the two subspaces[24], and if all the inner products between the basis vectors
siT sj (i = 1, ¼, Q; j = 1, ¼, P ) are small, the distance
between the two subspaces is large. If some of these innerproducts have large values, the distance between the two
subspaces can be small.
Therefore, when the time courses of the signal sources sp( p = 1, ¼, Q) and those of the interference sources sp( p = 1, ¼, P ) are very different and the correlations between
them are small, the distance between the signal and interference subspaces is large. In such cases, the overestimationof r may not be problematic, because the excess columns of G
are nearly orthogonal to the signal time-course vector sp( p = 1, ¼, Q), so the projector (I - GGT ) hardly affects BS.
The results of overestimating r at 12 in our computer simulation are shown in ﬁgures 18(c) and (d). Comparison of theseresults with the results with r=6 (ﬁgures 3(c) and 5(c))
shows that overestimation of r causes no serious problems.

This result can be explained by the arguments above, because
in our computer simulation, sp and sp were very different.
Many algorithms have been proposed by our group andothers for removing artifacts and interferences from sensor
measurements. One method, called common-mode subspace
projection (CSP), was reported by our group as a conference
proceedings paper [17]. The CSP algorithm uses a part of the
DSSP algorithm, and estimates the interference subspace as
an intersection between the row spans of two spatio-temporal
data matrices. It is, however, different from the DSSP algorithm in that the CSP algorithm requires the artifact-onlymeasurements to be conducted separately from the signal
measurements. The proposed DSSP algorithm does not
require such separate, additional measurements, and thus thetotal measurement time for the DSSP algorithm is half of thatfor the CSP method, which is advantageous in clinicalenvironments, because the burden for a patient is nearly
proportional to the total time for measurements.
Signal space projection (SSP) is a well-known method
for interference removal in MEG and EEG measurements[25, 26]. In the SSP method, the (spatial-domain) interferencesubspace projector PI is ﬁrst computed, and the sensor datab (t ) is then projected onto the subspace orthogonal to theinterference subspace by computing (I - PI ) b (t ). In order toderive PI , the SSP method requires the information on thelocations of interference sources. Therefore, the method
cannot be applied to the artifact removal problem in the spinal
cord imaging, because the spatial locations of the artifacts
cannot be determined. The proposed DSSP algorithm, in
15

J. Neural Eng. 13 (2016) 036007

K Sekihara et al

Figure 17. Results of source localization using the interference-removed sensor time courses shown in the bottom panels of ﬁgure 16. The
relative voxel intensity is color-coded according to the color-bar, and overlaid onto the patient’s MRI. The cross hair indicates the point withthe maximum reconstruction intensity, and three cross sectional MR images at this point are shown. (a) Results from the data taken with left
index ﬁnger stimulation. (b) Results from the data taken with right index ﬁnger stimulation.

contrast, does not require such information, and this fact, ingeneral, makes the DSSP algorithm much more useful thanthe SSP method, because such information is unavailable inmany real-life applications.
For interferences from brain stimulators such as the VNS
interference, the location of the interference source may bepredetermined. However, even in this case, a very accurate
determination is needed to implement SSP when the interference is much stronger than the signal. This fact can be seenusing our computer simulation in section 3.2. We applied SSP
to the interference-overlapped data shown in the top panel ofﬁgure 7. The results are shown in ﬁgure 19. In this ﬁgure, thetop panel shows the results obtained when the exact locationof the interference source is used. The middle and bottom
panels, respectively, show the results obtained when thelocation of the interference source has 1 and 3 mm errors.
These results show that even a 1 mm localization error affects
results of interference removal, suggesting that SSP requires avery accurate localization of an interference source. This fact
may make the application of SSP to removing large artifacts(such as the VNS interference) somewhat impractical.
In summary, this paper has proposed a novel algorithm toremove large interferences overlapped on biomagnetic data.
Since the algorithm makes use of the duality in deﬁning thesignal subspace, it is named DSSP. We have presented adetailed mathematical description of the proposed DSSP
algorithm, which we have validated by computer simulationand two kinds of biomagnetic data: spinal cord data obtained
from a healthy volunteer and MEG data taken from a patientwith a vagus nerve stimulator. The proposed DSSP algorithm

is suited to remove interferences from sources located fairly
close to the boundary of the source space. The eye-blink andcardiac artifacts encountered in MEG measurements are also
such interferences. The investigation on applying the DSSP
algorithm to the removal of these artifacts is planned andresults will be published.

Acknowledgments
This work has been supported by grants from the Japanese
Ministry of Education, Science, Culture and Sports(26282149, 26462231) and by a grant from SECOM Scienceand Technology Foundation.

Appendix
A.1. The row span of B S is equal to S

Here, we prove that, under the low-rank signal assumption(Q < M ), S is equal to the row span of BS. That isS = span {b1S, ¼, b SM}.

(45)

Equation (45) is equivalent to claiming that any row vector x( x Î S) can be expressed as a linear combination of the row
vectors, b1S, ¼, b SM . The claim is proved as follows. Sincex Î S, x can be expressed as a linear combination of sp :Q

x = å gp sp .
p=1

16

(46)

J. Neural Eng. 13 (2016) 036007

K Sekihara et al

Bx

5

By

5

y (cm)y (cm)

6

0

4

0

-5

-5
-5

0

5

x (cm)-50

5x (cm)Bz

5y (cm)

20-2

y (cm)

-4
-6
-8

0

(a) -5

-5

0

-4

-2

-6

-4

-2

0

2

4

6

8

0

2

4

6

8

x (cm)(b)5

x (cm)

-6

Bx

By

5

5

y (cm)y (cm)

6

0

4

0

-5

-5
-5

5

0

x (cm)-550x (cm)

0

-2

Bz

5

y (cm)2y (cm)

-4
-6
-8

0

x (cm)

-5

(c)

-5

5

0

x (cm)(d)Figure 18. (a) The ﬁeld contour maps at the latency of 5 ms for the artifact-removed results when μ and ν are set at 40. (b) The sourcereconstruction results using the ﬁeld maps shown in (a). (c) The ﬁeld contour maps at the latency of 5 ms for the artifact-removed results
when r was set at 12 (μ and ν were set at 20). (d) The source reconstruction results using the ﬁeld maps shown in (c).

If this vector x is also expressed as a linear combination ofb Sj , we havelinear equations:a1l11 +  + aM l1M = g1

M

x = å aj b Sj .

a1l21 +  + aM l2M = g2
a1lQ1 +  + aM lQM = gQ.

(47)

j=1

Substituting equation (14) into (47) and using equation (46),
we get
M

Q

Q

j=1

p=1

p=1

å aj å lpj sp = å gp sp .

Under the low-rank signal assumption M > Q, for any arbitrary
set of g1, ¼, gQ , we have a set of a1, ¼, aM that fulﬁlls the above
equations. Therefore, equation (47) always holds for any x( x Î S), and thus span {b1S, ¼, b SM} is equal to S.
Using the same arguments, it can be shown that the row
span of BI , span {b1I, ¼, b IM}, is equal to the interference sub~I~I
space I . Also, both the row span of PBI , span {b1 , ¼, bM},

(48)

Comparing the coefﬁcients of the vector sp in the left- andright-hand sides of the equation above, we have a set of Q
17

J. Neural Eng. 13 (2016) 036007

K Sekihara et al

Figure 19. Results of applying the signal space projection (SSP) tothe interference-overlapped data shown in the top panel of ﬁgure 7.
Top panel in this ﬁgure shows the results obtained when the accurate
location of the interference source is available. The moddle and thebottom panels, respectively, show results obtained when the locationhas 1 and 3 mm errors.
I

Figure 20. Cutoff property of the pseudo signal subspace projector P

computed using the geometry in our computer simulation in section
3.1. The ratio R = PlI (r )2 lI (r )2 is plotted with respect to y, thedistance from the center of the source space. The upper panel shows
the relative locations of the source space and the interference source.

I

and the row span of (I - P ) BI , span {b¯ 1 , ¼, b¯ M}, are equal
to I .

DSSP algorithm ﬁrst applies the pseudo subspace projector P
to create Bin and I− P to create Bout . Both Bin and Boutshould contain signiﬁcant amount of the interference, becauseof this ‘dull’ cutoff of the pseudo signal subspace projector P .
The DSSP algorithm makes use of this property of P and
detects the interference subspace by ﬁnding common components in the row spans of Bin and Bout .

A.2. Cutoff property of the pseudo signal subspace projector

Here, we present an example of the cutoff property of the pseudo
signal subspace projector. The example was derived assuming
the source-sensor geometry used in our computer simulation insection 3.1 where the augmented lead ﬁeld was computed over
the two-dimensional source space deﬁned by -8  x  8 cm ,-6  y  6 cm , and z=0 cm. The sensors are arranged on aplane of z=7 cm. Assuming that a source generating interference signal existed at r , the ratio R = PlI (r )2 lI (r )2 wascomputed, where lI (r ) indicates the lead ﬁeld vector of theinterference source. The ratio R expresses the decrease of thetotal power of the interference after processing by the projector
P . In this example, the source is located on the y-axis (i.e.,r = (0, y, 0)), and the conﬁguration of the source space and theinterference source is shown in the upper panel of ﬁgure 20. The
power decrease ratio R versus the y coordinate, which is the
distance of the interference source from the origin, is shown in
the lower part of ﬁgure 20.
The plot shows that the cutoff property of the pseudo
signal subspace projector P is not sharp. It passes almost
100% of the signal power if the interference source is locatedwithin the source space (y  6 cm) and attenuates only 5%–
20% of the signal power when the source is located in thevicinity of the source space (10  y  50 cm). The proposed

A.3. Proof of equation (26)Here, we present a proof of equation (26) in its general form.
That is, we deﬁne two subspaces,  and  , with basis vectors
x1, ¼, xD and y1 , ¼, y D¢ , respectively. We show that, if
x Î  = span {x1, ¼, xD } and y Î  = span {y1 , ¼, y D¢}, therelationship
z=x+yÎÈ

holds. The proof is straightforward. Since x Î span {x1, ¼, xD }and y Î span {y1 , ¼, y D¢}, we can write
D

D¢

j=1

j=1

z = x + y = å cj xj + å dj yj .

(49)

Therefore, we getz Î span {x1, ¼, xD , y1 , ¼, y D¢ }= span {x1, ¼, xD } È span {y1 , ¼, y D¢ }=  È .

18(50)

J. Neural Eng. 13 (2016) 036007

K Sekihara et al[13] Adachi Y, Kawai J, Miyamoto M, Uehara G, Kawabata S,
Okubo H, Fukuoka Y and Komori H 2003 Threedimensionally conﬁgured squid vector gradiometers forbiomagnetic measurement Supercond. Sci. Technol. 16 1442
[14] Wikso J P 1989 Biomagnetic sources and their models
Advances in Biomagnetism ed SJ Williamson et al (NewYork: Plenum) pp 1–8[15] Fukuoka Y, Komori H, Kawabata S, Ohkubo H,
Shinomiya K and Terasaki O 2002 Imaging of neuralconduction block by neuromagnetic recording Clin.
Neurophysiol. 113 1985–92[16] Fukuoka Y, Komori H, Kawabata S, Ohkubo H and Shinomiya K
2004 Visualization of incomplete conduction block by
neuromagnetic recording Clin. Neurophysiol. 115 2113–22[17] Watanabe T, Kawabata Y, Ukegawa D, Kawabata S,
Adachi Y and Sekihara K 2013 Removal of stimulusinduced artifacts in functional spinal cord imaging 2013 35th
Annual Int. Conf. IEEE Engineering in Medicine andBiology Society (EMBC) (Piscataway, NJ: IEEE) pp 3391–4[18] Sekihara K and Nagarajan S S 2015 Electromagnetic Brain
Imaging: A Bayesian Perspective (Berlin: Springer)[19] Sarvas J 1987 Basic mathematical and electromagnetic conceptsof the biomagnetic inverse problem Phys. Med. Biol. 32 11–22[20] Hämäläinen M S and Ilmoniemi R J 1984 Interpretingmeasured magnetic ﬁelds of the brain: estimates of current
distributions Technical Report TKK-F-A559 (Helsinki
University of Technology)[21] Hämäläinen M S and Ilmoniemi R J 1994 Interpreting
magnetic ﬁelds of the brain: minimum norm estimates Med.
Biol. Eng. Comput. 32 35–42[22] Pascual-Marqui R D 2002 Standardized low resolution brain
electromagnetic tomography (sLORETA): technical details
Methods Findings Exp. Clin. Pharmacology 24 5–12
[23] Scharf L L 1991 Statistical Signal Processing: Detection,
Estimation, and Time Series Analysis (New York: AddisonWesley Publishing Company)[24] Meyer C D 2000 Matrix Analysis and Applied Linear Algebra(Philadelphia, PA: Society for Industrial and Applied
Mathematics)[25] Uusitalo M and Ilmoniemi R 1997 Signal-space projection
method for separating meg or eeg into components Med.
Biol. Eng. Comput. 35 135–40[26] Nolte G and Curio G 1999 The effect of signal-space
projection on dipole localization accuracy IEEE Trans.
Biomed. Eng. 46 400–8

References[1] Nagarajan S S, Attias H T, Hild K E and Sekihara K 2007 A
probabilistic algorithm for robust interference suppression
in bioelectromagnetic sensor data Stat. Med.
26 3886–910[2] Nagarajan S S, Attias H T, Hild K E II and Sekihara K 2006 A
graphical model for estimating stimulus-evoked brain
responses from magnetoencephalography data withlarge background brain activity NeuroImage
30 400–16[3] de Cheveigné A and Parra L C 2014 Joint decorrelation, aversatile tool for multichannel data analysis Neuroimage 98
487–505[4] Shinomiya K, Furuya K, Sato R, Okamoto A, Kurosa Y andFuchioka M 1998 Electrophysiologic diagnosis of cervical
OPLL myelopathy using evoked spinal cord potentials Spine
13 1225–33[5] Adachi Y, Oyama D, Kawabata S, Sekihara K, Haruta Y and
Uehara G 2013 Magnetospinography: instruments andapplication to functional imaging of spinal cords IEICE
Trans. Electron. 96 326–33[6] Kawabata S, Komori H, Mochida K, Ohkubo H andShinomiya K 2003 Visualization of conductive spinal cord
activity using a bio-magnetometer Spine 27 475–9[7] Adachi Y, Kawai J, Miyamoto M, Ogata H, Tomori M,
Kawabata S, Sato T and Uehara G 2009 A SQUID system formeasurement of spinal cord evoked ﬁeld of supine subjects
IEEE Trans. Appl. Supercond. 19 861–6[8] Sato T, Adachi Y, Tomori M, Ishii S, Kawabata S andSekihara K 2009 Functional imaging of spinal cord electricalactivity from its evoked magnetic ﬁeld IEEE Trans. Biomed.
Eng. 56 2452–60[9] Kumihashi I and Sekihara K 2010 Array-gain constraint
minimum-norm spatial ﬁlter with recursively updated grammatrix for biomagnetic source imaging IEEE Trans. Biomed.
Eng. 57 1358–65[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 and Nagarajan S S 2008 Adaptive Spatial Filters
for Electromagnetic Brain Imaging (Berlin: Springer)[12] Golub G H and Van Loan C F 2012 Matrix Computations vol 3(Baltimore, MD: Johns Hopkins University Press)

19