# Signal Space Separation Method for a Biomagnetic Sensor Array Arranged on a Flat Plane for Magnetocardiographic Applications: A Computer Simulation Study

Authors: Kensuke Sekihara
Journal: Journal of Healthcare Engineering
 (2018)

---

## Abstract

Although the signal space separation (SSS) method can successfully suppress interference/artifacts overlapped onto magnetoencephalography (MEG) signals, the method is considered inapplicable to data from nonhelmet-type sensor arrays, such as the flat sensor arrays typically used in magnetocardiographic (MCG) applications. This paper shows that the SSS method is still effective for data measured from a (nonhelmet-type) array of sensors arranged on a flat plane. By using computer simulations, it is shown that the optimum location of the origin can be determined by assessing the dependence of signal and noise gains of the SSS extractor on the origin location. The optimum values of the parametersLC and LD, which, respectively, indicate the truncation values of the multipole-order ℓ of the internal and external subspaces, are also determined by evaluating dependences of the signal, noise, and interference gains (i.e., the shield factor) on these parameters. The shield factor exceeds 104 for interferences originating from fairly distant sources. However, the shield factor drops to approximately 100 when calibration errors of 0.1% exist and to 30 when calibration errors of 1% exist. The shielding capability can be significantly improved using vector sensors, which measure the x, y, and z components of the magnetic field. With 1% calibration errors, a vector sensor array still maintains a shield factor of approximately 500. It is found that the SSS application to data from flat sensor arrays causes a distortion in the signal magnetic field, but it is shown that the distortion can be corrected by using an SSS-modified sensor lead field in the voxel space analysis.


---

## Body

## 1. Introduction

Development of a sensor system that can measure biomagnetic signals in room-temperature environments has gained great interest. One promising candidate among room-temperature sensors for biomagnetic systems is magnetoresistive (MR) sensors [1–4], which can lead to the development of novel low-initial-cost and maintenance-free biomagnetic systems. A potential near-future application of such systems is a low-cost magnetocardiography (MCG) system using MR sensors [5]. Such low-cost and maintenance-free MCG systems could replace the 12-lead electrocardiogram (ECG) now routinely used in daily clinical examinations.However, to develop such low-cost systems, one major problem is the removal of ambient noise magnetic fields that exist in urban hospital environments. Biomagnetic signals are many orders of magnitude weaker than these ambient interference magnetic fields, called environmental noise. To reduce the influence of such environmental noise, biomagnetic measurements have traditionally relied on two kinds of hardware-based solutions: one is magnetically shielded rooms (MSRs) and the other is gradiometer sensors [6].According to a purely technical point of view, use of an MSR would be advantageous even for MR sensor systems. However, the use of a high- or a medium-quality MSR may invalidate our goal, which is developing low-initial-cost biomagnetic systems, because MSRs are, in general, very costly and the use of even a medium-quality MSR causes a considerable increase of the total initial cost of the system. Gradiometers can significantly reduce the environmental noise, and the reduction ratio of a typical first-order gradiometer is believed to reach almost two orders of magnitude [6], although the noise reduction capability strongly depends on the precision in sensor manufacturing. However, for MR sensor systems, it is currently unknown whether gradiometers can be incorporated into the MR sensor hardware design. No gradiometer-type MR sensors have been developed so far.Therefore, to attain the goal of developing low-cost biomagnetic systems, it may not be possible to rely on traditional hardware-based methods. This paper addresses a third option, developing software shielding methods, namely, efficient signal processing methods for environmental noise cancellation. A number of signal processing methods have been developed for this purpose, and arguments for and against those methods can be found in [7]. This paper focuses on a method called signal space separation (SSS), which was originally proposed for environmental noise cancellation for magnetoencephalography (MEG) SQUID sensor arrays [8–10].The SSS method has an excellent characteristic that it imposes almost no prerequisites on the data or sources (i.e., it does not use any restrictive data or source models). One mild prerequisite of the SSS method, which can naturally be fulfilled, is that the region in which sensors are installed should be source-free (i.e., no current sources in that region). Under this assumption, the method decomposes the measured data into two components originated from the so-called “internal” and “external” regions. The internal region refers to a region that is closer to the origin than the sensors are, and the external region refers to the one that is farther from the origin than the sensors are.The method can efficiently suppress interferences overlapped onto biomagnetic signals if clever choices of the origin location can make the internal and external regions match the signal and interference regions, respectively. It is not difficult to find such origin locations for helmet-type sensor arrays used in MEG. However, the SSS method has not been considered applicable to data from nonhelmet-type sensor arrays, such as sensor arrays arranged on a flat plane, which are usually used for MCG systems, because it does not seem possible to find an appropriate origin location for those nonhelmet-type arrays.This paper presents a computer simulation-based investigation that explores the possibility of applying the SSS method to data measured from an array of sensors arranged on a flat plane. The goal of the investigation is to show that the SSS method is still effective for data from such flat sensor arrays. A series of computer simulations are performed assuming flat sensor arrays whose sensor arrangement is typical in MCG applications [11–15]. Using the results of these computer simulations, this paper seeks optimal values for numerical parameters used in the SSS method, showing that their choices are crucial for the effective use of the SSS method when applied to flat sensor arrays. It is found that the application of the SSS method to flat sensor data causes a distortion of signals but that this distortion can be corrected by using the SSS-modified lead field in the voxel space analysis.This paper is organized as follows. In Section2, the SSS method is described in detail, including the analysis of the problems caused when the SSS method is applied to flat sensor data. Section 3 presents computer simulation-based investigation, which shows the effectiveness of the SSS method for data from flat sensor arrays. This section includes empirical determinations of SSS parameters crucial to attain optimal performance of the SSS method. Section 4 summarizes the findings of the investigation.


## 2. Signal Space Separation Method

### 2.1. Data Model

Biomagnetic measurement is conducted using a sensor array, which simultaneously measures the signal with multiple sensors. Let us define the measurement of themth sensor as ym. The measurement from the whole sensor array is expressed as a column vector y:y=y1,y2,…,yMT. Here, M is the number of sensors, and the superscript T indicates the matrix transpose. Throughout this paper, plain italics indicates scalars, lowercase boldface indicates vectors, and uppercase boldface indicates matrices.The location in the three-dimensional space is represented byr:r=x,y,z. The source magnitude at r is denoted by a scalar sr. The source vector is denoted by s(r), and the source orientation is denoted by η=ηx,ηy,ηzT. We thus have the relationship: sr=srη. 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, denoted by lmxr,lmyr and lmzr. Let us define an M×3 matrix Lr whose mth row is equal to lmxr,lmyr,lmzr. This matrix Lr, referred to as the lead field matrix, represents the sensitivity of the sensor array at r. When the unit-magnitude source at r is oriented in the η direction, the outputs of the sensor array are expressed as lr=Lrη. This column vector lr, referred to as the lead field vector, represents the sensitivity of the sensor array in the direction of η at the location r.The outputs of the sensor arrayy are expressed as the sum of a magnetic signal b and additive sensor noise, represented by a random vector ε
(1)y=b+ε.Here, the magnetic signalb is expressed as
(2)b=bS+bI.Here,bS, called the signal vector, represents the biomagnetic signal that is the target of the measurements, and bI, called the interference vector, represents the interference overlapped onto the signal bS. In this paper, the interference bI represents so-called environmental noise, and sources of environmental noise are assumed to be located much farther from the sensors than the sources of interest are. Sources of environmental noise are, in general, located from several meters (in cases of the noise sources such as electronic appliances in a laboratory) to several kilometers (in cases of urban environmental noise sources such as the subway noise) distant from the magnetically shielded room.The sources that generatebS are confined to a region called the source space (e.g., the source space is the cardiac region for MCG and the brain region for MEG measurements). Let us assume that a total of Q discrete sources exist in the source space. Their locations are denoted by r1,…,rQ, their orientations by η1,…,ηQ, and their magnitudes by s1,…,sQ. Then, the source distribution is expressed as
(3)sr,t=∑q=1Qsqηqδr−rq,where δr indicates the Dirac delta function. The signal vector bS is expressed as
(4)bS=∫ΩLr∑q=1Qsqηqδr−rqdr=∑q=1Qsqlq,where lq represents the lead field vector of the qth source obtained such that lq=Lrqηq.


### 2.2. Derivation of SSS Basis Vectors

One fundamental assumption of the SSS method is that the sensors are installed in a source-free region, which is referred to as the sensor region. Then, the magnetic field atr, B(r), is expressed using the spherical polar coordinate r=r,θ,ϕ by
(5)Br=−μ0∑ℓ=1∞∑m=−ℓℓαℓ,mνℓ,mθ,ϕrℓ+2−μ0∑ℓ=1∞∑m=−ℓℓβℓ,mrℓ−1ωℓ,mθ,ϕ,where μ0 indicates the magnetic permeability of free space. In (5), νℓ,mθ,ϕ and ωℓ,mθ,ϕ are the modified vector spherical harmonics [8, 16]. The index parameter ℓ is called the multipole-order or multipole parameter. In the right-hand side of (5), the first term represents the magnetic field generated from sources located closer to the origin than the sensors are.The second term represents the magnetic field from sources located farther from the origin than the sensors are. The region closer to the origin than the sensors is referred to as the internal region, and the region farther from the origin than the sensors is referred to as the external region. Let us define the polar radial coordinate of the sensor nearest to the origin asrDmin and the radial coordinates of the sensor farthest from the origin as rDmax. The internal region is formally defined as the region with r<rDmin and the external region as the region with r>rDmax. The region with rDmin<r<rDmax is called the intermediate region.Let us derive the SSS basis vectors. To do so, the magnetic signal detected by thejth sensor is denoted by bj and the location and the normal vector of the jth sensor by rj and ζj. Then, we have
(6)bj=Brj⋅ζj=bintj+bextj,where the notation “·” indicates taking the inner product between two vectors (note that when the area of the pickup coils is taken into consideration, the sensor signal bj is obtained as a surface integral of Brj·ζj over the area of the jth pick-up coil). Here, bintj and bextj, respectively, represent magnetic components originating from the internal and external regions. These components are expressed as
(7)bintj=−∑ℓ=1∞∑m=−ℓℓαℓ,mνℓ,mθj,ϕj⋅ζjrjℓ+2,bextj=−∑ℓ=1∞∑m=−ℓℓβℓ,mrjℓ−1ωℓ,mθj,ϕj⋅ζj,where we set μ0=1 for simplicity. Let us define the internal and external components of the vector b as bint=bint1,…,bintMT and bext=bext1,…,bextMT, which are expressed such that
(8)bint=∑ℓ=1∞∑m=−ℓℓαℓ,mcℓ,m,bext=∑ℓ=1∞∑m=−ℓℓβℓ,mdℓ,m,where column vectors cℓ,m and dℓ,m are given by
(9)cℓ,m=1r1ℓ+2νℓ,mθ1,ϕ1⋅ζ1⋮1rMℓ+2νℓ,mθM,ϕM⋅ζM,dℓ,m=r1ℓ−1ωℓ,mθ1,ϕ1⋅ζ1⋮rMℓ−1ωℓ,mθM,ϕM⋅ζM.Truncating the summation with respect to the multiple orderℓ to LC for bint and LD for bext, we finally obtain
(10)b=bint+bext=∑ℓ=1LC∑m=−ℓℓαℓ,mcℓ,m+∑ℓ=1LD∑m=−ℓℓβℓ,mdℓ,m.Thus, defining(11)C=c1,−1,c1,0,c1,1,…,cLC,LC,D=d1,−1,d1,0,d1,1,…,dLD,LD,α=α1,−1,α1,0,α1,1,…,αLC,LCT,β=β1,−1,β1,0,β1,1,…,βLD,LDT,we obtain
(12)b=bint+bext=Cα+Dβ=C,Dαβ=Sx,where S=C,D and x=αT,βTT. Here, C is an M×NC matrix, and D is an M×ND matrix, where
(13)NC=LC2+2LC,ND=LD2+2LD.Note that the truncation valuesLC and LD correspond to the highest spatial frequencies possibly contained in bint and bext [9, 17], respectively. Therefore, setting these parameters at too low values may result in an insufficient representation of the signal vectors bint and bext. The effects of LC and LD for data from the 306-channel Elekta Neuromag have been investigated and values of LC=8 and LD=3 were found to be sufficient for such data sets in [8, 9].


### 2.3. SSS Signal Extractors

Equation (12) is the basis for estimating the internal and external components bint and bext from given magnetic signal data b. That is, the least squares estimate x^=α^T,β^TT is obtained as
(14)x^=STS−1STb.Then,bint and bext are estimated as
(15)b^int=Cα^,(16)b^ext=Dβ^.We now derive SSS signal extractors and rewrite (15) and (16) using these extractors. To do so, let us define an operation to make a new column vector ai,…,ajT by using the ith to jth components of a=a1,…,aMT as ai:j (namely, ai:j=ai,…,ajT). From (14), we have
(17)α^=STS−1STb1:NC.With a small positive constantκ, the relationship
(18)α^=STS−1STb1:NC≈STS+κI−1STb1:NC,holds. Then, using the matrix inversion formula
(19)STS+κI−1STb=STSST+κI−1b,we get
(20)α^≈STS+κI−1STb1:NC=STSST+κI−1b1:NC=CTSST+κI−1b≈CTSST−1b=CTCCT+DDT−1b.Using (15) and (20), we obtain
(21)b^int≈CCTCCT+DDT−1b.Thus, the internal componentbint can be extracted by multiplying
(22)Pint=CCTCCT+DDT−1,with the magnetic-field data b. That is, the matrix Pint acts as a projector that passes the internal components and blocks the external ones (note that since Pint2=Pint and PintT=Pint do not hold, Pint is not actually a projector). Therefore, we call Pint the SSS signal extractor in this paper.In exactly the same manner, we can derive(23)b^ext=Dβ^≈DDTCCT+DDT−1b,and the SSS external-signal extractor Pext is derived as
(24)Pext=DDTCCT+DDT−1.ThisPext passes the external components but blocks the internal ones.


### 2.4. Interference Suppression

A key condition for the success of the SSS interference suppression is that 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. A typical configuration between the helmet-type sensor array and the source space is depicted in Figure 1(a). As can be seen in this figure, an appropriate location of the origin can be found so that the internal region covers the whole source space and the external region covers all locations of interference sources. (Note that, in this paper, interference indicates only environmental noise, and its sources are assumed to be located much farther from the sensors than the signal sources.)Figure 1
(a) A typical configuration of the internal region relative to the source space for a helmet-type sensor array. As can be seen, an appropriate location of the origin can be found such that the internal region covers the whole source space. (b) A possible configuration of the internal region relative to the source space in case of a flat sensor array. The internal region cannot entirely cover the source space, which extends into an intermediate region. In these figures, the inner circle with a broken line indicates the boundary between the internal and intermediate regions. The outer circle with a broken line indicates the boundary between the intermediate and external regions.
(a)
(b)When this key condition is met, (10) provides a natural separation between the signal and interference. That is, when the source space is included within the internal region, the relationship
(25)bS∈spanC,holds, where the notation span C indicates the span of the column vectors of C. This span C is referred to as the internal subspace in this paper. That is, since the signal vector bs belongs to the internal subspace, bs is expressed as a linear combination of the column vectors of C
(26)bS=∑j=1NCαjcj=Cα,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,…,αNCT). Thus, denoting an ND×1 column vector whose elements are all zero by 0, we can derive the relationship
(27)PintbS=PintCα=PintSα0=CCTSST−1Sα0=CSTSST−1Sα01:NC=CSTS−1STSα0=Cα=bS.The equation above indicates that the SSS signal extractorPint passes the signal vector bS with no distortion.The assumption that the interference sources are located within the external region leads to(28)bI∈spanD,where spanD is referred to as the external subspace. Since the interference vector bI belongs to the external subspace, the interference vector bI is expressed as
(29)bI=∑j=1NDβjdj=Dβ,where dj is the jth column of D, βj is the jth expansion coefficient, and β is a column vector containing coefficients β=β1,…,βNDT . Again denoting the NC×1 column vector whose elements are all zero by 0, we have the relationship
(30)PintbI=PintDβ=PintS0β=CSTS−1STS0β1:NC=C0=0.The equation above indicates that the extractorPint completely blocks the interference vector bI.Consequently, using (1) and (2), we show
(31)Pinty=PintbS+Pint+bI+Pintε=bS+ε′.The equation above indicates that, by multiplying the extractorPint with the data vector y, the signal vector bS is selectively extracted with no distortion.In (31), ε′:ε′=Pintε indicates the noise in the SSS-cleaned data. Assuming that the sensor noise ε is Gaussian with the covariance matrix σ2I, the covariance matrix of ε′,Σε′ is derived such that
(32)Σε′=ε′ε′T=PintεεTPintT=σ2PintPintT,where we have εεT=σ2I and the notation ⋅ indicates averaging. The equation above shows that PintPintT can be considered as the noise gain of the SSS interference suppression process. Particularly, since the diagonal elements of Σε′ expresses the gain relationship between the variances of the input and output noises, we define the noise gain Gε such that
(33)Gε=1MΣj=1MΣε′j,j,where Σε′j,j indicates the jth diagonal element of the matrix Σε′.


### 2.5. SSS Method for Flat Sensor Arrays

In Figure1(b), a possible configuration of the internal region relative to the source space and sensors is depicted for the case of a flat sensor array. As shown here, the internal region may not cover the entire source space, and the source space is extended into the intermediate region. As a result, the lead field vector of a signal source, lq, may have components expanded by the columns of D, as well as components expanded by the columns of C, resulting in
(34)lq=Cα′+Dβ′,where α′ and β′ are column vectors containing the expansion coefficients.Therefore, by applying the SSS extractorPint to lq, we have
(35)Pintlq=PintCα′+Dβ′=Cα′=l~q.That is, the extractorPint changes the lead field of this signal source from lq to l~q. Assuming that Q signal sources exist, the signal vector bSt is expressed as
(36)bSt=∑q=1Qsqtlq.This signal vector changes tob~St, which is given by
(37)b~St=PintbSt=Pint∑q=1Qsqtlq=∑q=1QsqtPintlq=∑q=1Qsqtl~q.It is clear here that applying the SSS extractorPint distorts the signal vector bSt. This is a problem that occurs when the SSS method is applied to data measured with a flat sensor array. Computer simulation-based investigation of the signal distortion and its correction are given in Section 3.6.As can be seen in Figures1(a) and 1(b), the external region does not differ between helmet and flat sensor arrays. In this paper, we assume that all interference is environmental noise and no interference sources exist in the vicinity of sensors. Since, under this assumption, we assume that all interference sources are located in the external region, the interference vector bIt does not have components belonging to span C, and applying Pint to the data vector removes the interference vector.


### 2.6. Evaluation of the SSS Method’s Performance

As discussed in the preceding sections, a flat sensor array imposes nonideal conditions on the SSS interference suppression, and (25) is never fulfilled. In addition, the existence of sensor calibration errors (which will be considered in Section 3.4) could, to some extent, invalidate the assumption in (28). Therefore, for data from flat sensor arrays, the relationships
(38)PintbS=bS,PintbI=0,will never be attained.We can evaluate the performance of the SSS method for data from flat sensor arrays by checking how close the SSS-processed results come to (38). Namely, we define the performance measures such that
(39)GS=PintbSbS,(40)GI=PintbIbI.In the equations above,GS is called the signal gain. Ideally, GS is equal to 1, and deviation of GS from 1 is a measure of performance degradation of the SSS method. GI is called the interference gain. Ideally, GI is equal to zero, indicating that the method completely blocks the interference. Thus, a deviation of GI from 0 is a measure of the performance degradation. Note that 1/GI has been often called the shield factor in the previous literature [9, 18]. We use GS and GI (or 1/GI, the shield factor), as well as the noise gain Gε in (33) to evaluate the performance of the SSS method when it is applied to data from flat sensor arrays in Section 3.


## 2.1. Data Model

Biomagnetic measurement is conducted using a sensor array, which simultaneously measures the signal with multiple sensors. Let us define the measurement of themth sensor as ym. The measurement from the whole sensor array is expressed as a column vector y:y=y1,y2,…,yMT. Here, M is the number of sensors, and the superscript T indicates the matrix transpose. Throughout this paper, plain italics indicates scalars, lowercase boldface indicates vectors, and uppercase boldface indicates matrices.The location in the three-dimensional space is represented byr:r=x,y,z. The source magnitude at r is denoted by a scalar sr. The source vector is denoted by s(r), and the source orientation is denoted by η=ηx,ηy,ηzT. We thus have the relationship: sr=srη. 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, denoted by lmxr,lmyr and lmzr. Let us define an M×3 matrix Lr whose mth row is equal to lmxr,lmyr,lmzr. This matrix Lr, referred to as the lead field matrix, represents the sensitivity of the sensor array at r. When the unit-magnitude source at r is oriented in the η direction, the outputs of the sensor array are expressed as lr=Lrη. This column vector lr, referred to as the lead field vector, represents the sensitivity of the sensor array in the direction of η at the location r.The outputs of the sensor arrayy are expressed as the sum of a magnetic signal b and additive sensor noise, represented by a random vector ε
(1)y=b+ε.Here, the magnetic signalb is expressed as
(2)b=bS+bI.Here,bS, called the signal vector, represents the biomagnetic signal that is the target of the measurements, and bI, called the interference vector, represents the interference overlapped onto the signal bS. In this paper, the interference bI represents so-called environmental noise, and sources of environmental noise are assumed to be located much farther from the sensors than the sources of interest are. Sources of environmental noise are, in general, located from several meters (in cases of the noise sources such as electronic appliances in a laboratory) to several kilometers (in cases of urban environmental noise sources such as the subway noise) distant from the magnetically shielded room.The sources that generatebS are confined to a region called the source space (e.g., the source space is the cardiac region for MCG and the brain region for MEG measurements). Let us assume that a total of Q discrete sources exist in the source space. Their locations are denoted by r1,…,rQ, their orientations by η1,…,ηQ, and their magnitudes by s1,…,sQ. Then, the source distribution is expressed as
(3)sr,t=∑q=1Qsqηqδr−rq,where δr indicates the Dirac delta function. The signal vector bS is expressed as
(4)bS=∫ΩLr∑q=1Qsqηqδr−rqdr=∑q=1Qsqlq,where lq represents the lead field vector of the qth source obtained such that lq=Lrqηq.


## 2.2. Derivation of SSS Basis Vectors

One fundamental assumption of the SSS method is that the sensors are installed in a source-free region, which is referred to as the sensor region. Then, the magnetic field atr, B(r), is expressed using the spherical polar coordinate r=r,θ,ϕ by
(5)Br=−μ0∑ℓ=1∞∑m=−ℓℓαℓ,mνℓ,mθ,ϕrℓ+2−μ0∑ℓ=1∞∑m=−ℓℓβℓ,mrℓ−1ωℓ,mθ,ϕ,where μ0 indicates the magnetic permeability of free space. In (5), νℓ,mθ,ϕ and ωℓ,mθ,ϕ are the modified vector spherical harmonics [8, 16]. The index parameter ℓ is called the multipole-order or multipole parameter. In the right-hand side of (5), the first term represents the magnetic field generated from sources located closer to the origin than the sensors are.The second term represents the magnetic field from sources located farther from the origin than the sensors are. The region closer to the origin than the sensors is referred to as the internal region, and the region farther from the origin than the sensors is referred to as the external region. Let us define the polar radial coordinate of the sensor nearest to the origin asrDmin and the radial coordinates of the sensor farthest from the origin as rDmax. The internal region is formally defined as the region with r<rDmin and the external region as the region with r>rDmax. The region with rDmin<r<rDmax is called the intermediate region.Let us derive the SSS basis vectors. To do so, the magnetic signal detected by thejth sensor is denoted by bj and the location and the normal vector of the jth sensor by rj and ζj. Then, we have
(6)bj=Brj⋅ζj=bintj+bextj,where the notation “·” indicates taking the inner product between two vectors (note that when the area of the pickup coils is taken into consideration, the sensor signal bj is obtained as a surface integral of Brj·ζj over the area of the jth pick-up coil). Here, bintj and bextj, respectively, represent magnetic components originating from the internal and external regions. These components are expressed as
(7)bintj=−∑ℓ=1∞∑m=−ℓℓαℓ,mνℓ,mθj,ϕj⋅ζjrjℓ+2,bextj=−∑ℓ=1∞∑m=−ℓℓβℓ,mrjℓ−1ωℓ,mθj,ϕj⋅ζj,where we set μ0=1 for simplicity. Let us define the internal and external components of the vector b as bint=bint1,…,bintMT and bext=bext1,…,bextMT, which are expressed such that
(8)bint=∑ℓ=1∞∑m=−ℓℓαℓ,mcℓ,m,bext=∑ℓ=1∞∑m=−ℓℓβℓ,mdℓ,m,where column vectors cℓ,m and dℓ,m are given by
(9)cℓ,m=1r1ℓ+2νℓ,mθ1,ϕ1⋅ζ1⋮1rMℓ+2νℓ,mθM,ϕM⋅ζM,dℓ,m=r1ℓ−1ωℓ,mθ1,ϕ1⋅ζ1⋮rMℓ−1ωℓ,mθM,ϕM⋅ζM.Truncating the summation with respect to the multiple orderℓ to LC for bint and LD for bext, we finally obtain
(10)b=bint+bext=∑ℓ=1LC∑m=−ℓℓαℓ,mcℓ,m+∑ℓ=1LD∑m=−ℓℓβℓ,mdℓ,m.Thus, defining(11)C=c1,−1,c1,0,c1,1,…,cLC,LC,D=d1,−1,d1,0,d1,1,…,dLD,LD,α=α1,−1,α1,0,α1,1,…,αLC,LCT,β=β1,−1,β1,0,β1,1,…,βLD,LDT,we obtain
(12)b=bint+bext=Cα+Dβ=C,Dαβ=Sx,where S=C,D and x=αT,βTT. Here, C is an M×NC matrix, and D is an M×ND matrix, where
(13)NC=LC2+2LC,ND=LD2+2LD.Note that the truncation valuesLC and LD correspond to the highest spatial frequencies possibly contained in bint and bext [9, 17], respectively. Therefore, setting these parameters at too low values may result in an insufficient representation of the signal vectors bint and bext. The effects of LC and LD for data from the 306-channel Elekta Neuromag have been investigated and values of LC=8 and LD=3 were found to be sufficient for such data sets in [8, 9].


## 2.3. SSS Signal Extractors

Equation (12) is the basis for estimating the internal and external components bint and bext from given magnetic signal data b. That is, the least squares estimate x^=α^T,β^TT is obtained as
(14)x^=STS−1STb.Then,bint and bext are estimated as
(15)b^int=Cα^,(16)b^ext=Dβ^.We now derive SSS signal extractors and rewrite (15) and (16) using these extractors. To do so, let us define an operation to make a new column vector ai,…,ajT by using the ith to jth components of a=a1,…,aMT as ai:j (namely, ai:j=ai,…,ajT). From (14), we have
(17)α^=STS−1STb1:NC.With a small positive constantκ, the relationship
(18)α^=STS−1STb1:NC≈STS+κI−1STb1:NC,holds. Then, using the matrix inversion formula
(19)STS+κI−1STb=STSST+κI−1b,we get
(20)α^≈STS+κI−1STb1:NC=STSST+κI−1b1:NC=CTSST+κI−1b≈CTSST−1b=CTCCT+DDT−1b.Using (15) and (20), we obtain
(21)b^int≈CCTCCT+DDT−1b.Thus, the internal componentbint can be extracted by multiplying
(22)Pint=CCTCCT+DDT−1,with the magnetic-field data b. That is, the matrix Pint acts as a projector that passes the internal components and blocks the external ones (note that since Pint2=Pint and PintT=Pint do not hold, Pint is not actually a projector). Therefore, we call Pint the SSS signal extractor in this paper.In exactly the same manner, we can derive(23)b^ext=Dβ^≈DDTCCT+DDT−1b,and the SSS external-signal extractor Pext is derived as
(24)Pext=DDTCCT+DDT−1.ThisPext passes the external components but blocks the internal ones.


## 2.4. Interference Suppression

A key condition for the success of the SSS interference suppression is that 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. A typical configuration between the helmet-type sensor array and the source space is depicted in Figure 1(a). As can be seen in this figure, an appropriate location of the origin can be found so that the internal region covers the whole source space and the external region covers all locations of interference sources. (Note that, in this paper, interference indicates only environmental noise, and its sources are assumed to be located much farther from the sensors than the signal sources.)Figure 1
(a) A typical configuration of the internal region relative to the source space for a helmet-type sensor array. As can be seen, an appropriate location of the origin can be found such that the internal region covers the whole source space. (b) A possible configuration of the internal region relative to the source space in case of a flat sensor array. The internal region cannot entirely cover the source space, which extends into an intermediate region. In these figures, the inner circle with a broken line indicates the boundary between the internal and intermediate regions. The outer circle with a broken line indicates the boundary between the intermediate and external regions.
(a)
(b)When this key condition is met, (10) provides a natural separation between the signal and interference. That is, when the source space is included within the internal region, the relationship
(25)bS∈spanC,holds, where the notation span C indicates the span of the column vectors of C. This span C is referred to as the internal subspace in this paper. That is, since the signal vector bs belongs to the internal subspace, bs is expressed as a linear combination of the column vectors of C
(26)bS=∑j=1NCαjcj=Cα,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,…,αNCT). Thus, denoting an ND×1 column vector whose elements are all zero by 0, we can derive the relationship
(27)PintbS=PintCα=PintSα0=CCTSST−1Sα0=CSTSST−1Sα01:NC=CSTS−1STSα0=Cα=bS.The equation above indicates that the SSS signal extractorPint passes the signal vector bS with no distortion.The assumption that the interference sources are located within the external region leads to(28)bI∈spanD,where spanD is referred to as the external subspace. Since the interference vector bI belongs to the external subspace, the interference vector bI is expressed as
(29)bI=∑j=1NDβjdj=Dβ,where dj is the jth column of D, βj is the jth expansion coefficient, and β is a column vector containing coefficients β=β1,…,βNDT . Again denoting the NC×1 column vector whose elements are all zero by 0, we have the relationship
(30)PintbI=PintDβ=PintS0β=CSTS−1STS0β1:NC=C0=0.The equation above indicates that the extractorPint completely blocks the interference vector bI.Consequently, using (1) and (2), we show
(31)Pinty=PintbS+Pint+bI+Pintε=bS+ε′.The equation above indicates that, by multiplying the extractorPint with the data vector y, the signal vector bS is selectively extracted with no distortion.In (31), ε′:ε′=Pintε indicates the noise in the SSS-cleaned data. Assuming that the sensor noise ε is Gaussian with the covariance matrix σ2I, the covariance matrix of ε′,Σε′ is derived such that
(32)Σε′=ε′ε′T=PintεεTPintT=σ2PintPintT,where we have εεT=σ2I and the notation ⋅ indicates averaging. The equation above shows that PintPintT can be considered as the noise gain of the SSS interference suppression process. Particularly, since the diagonal elements of Σε′ expresses the gain relationship between the variances of the input and output noises, we define the noise gain Gε such that
(33)Gε=1MΣj=1MΣε′j,j,where Σε′j,j indicates the jth diagonal element of the matrix Σε′.


## 2.5. SSS Method for Flat Sensor Arrays

In Figure1(b), a possible configuration of the internal region relative to the source space and sensors is depicted for the case of a flat sensor array. As shown here, the internal region may not cover the entire source space, and the source space is extended into the intermediate region. As a result, the lead field vector of a signal source, lq, may have components expanded by the columns of D, as well as components expanded by the columns of C, resulting in
(34)lq=Cα′+Dβ′,where α′ and β′ are column vectors containing the expansion coefficients.Therefore, by applying the SSS extractorPint to lq, we have
(35)Pintlq=PintCα′+Dβ′=Cα′=l~q.That is, the extractorPint changes the lead field of this signal source from lq to l~q. Assuming that Q signal sources exist, the signal vector bSt is expressed as
(36)bSt=∑q=1Qsqtlq.This signal vector changes tob~St, which is given by
(37)b~St=PintbSt=Pint∑q=1Qsqtlq=∑q=1QsqtPintlq=∑q=1Qsqtl~q.It is clear here that applying the SSS extractorPint distorts the signal vector bSt. This is a problem that occurs when the SSS method is applied to data measured with a flat sensor array. Computer simulation-based investigation of the signal distortion and its correction are given in Section 3.6.As can be seen in Figures1(a) and 1(b), the external region does not differ between helmet and flat sensor arrays. In this paper, we assume that all interference is environmental noise and no interference sources exist in the vicinity of sensors. Since, under this assumption, we assume that all interference sources are located in the external region, the interference vector bIt does not have components belonging to span C, and applying Pint to the data vector removes the interference vector.


## 2.6. Evaluation of the SSS Method’s Performance

As discussed in the preceding sections, a flat sensor array imposes nonideal conditions on the SSS interference suppression, and (25) is never fulfilled. In addition, the existence of sensor calibration errors (which will be considered in Section 3.4) could, to some extent, invalidate the assumption in (28). Therefore, for data from flat sensor arrays, the relationships
(38)PintbS=bS,PintbI=0,will never be attained.We can evaluate the performance of the SSS method for data from flat sensor arrays by checking how close the SSS-processed results come to (38). Namely, we define the performance measures such that
(39)GS=PintbSbS,(40)GI=PintbIbI.In the equations above,GS is called the signal gain. Ideally, GS is equal to 1, and deviation of GS from 1 is a measure of performance degradation of the SSS method. GI is called the interference gain. Ideally, GI is equal to zero, indicating that the method completely blocks the interference. Thus, a deviation of GI from 0 is a measure of the performance degradation. Note that 1/GI has been often called the shield factor in the previous literature [9, 18]. We use GS and GI (or 1/GI, the shield factor), as well as the noise gain Gε in (33) to evaluate the performance of the SSS method when it is applied to data from flat sensor arrays in Section 3.


## 3. Computer Simulation

### 3.1. Problems with Flat Sensor Data for SSS Applications

In order to show that the SSS method is still effective for data from sensors arranged on a flat plane, a series of computer simulation has been performed. First, we clarify the problems caused when the SSS method is applied to flat sensor data by comparing two cases of SSS application: one in which a helmet-type sensor array is used and the other in which a flat sensor array is used. For helmet-type sensor arrays, the key condition for the SSS method can be fulfilled, that is, one can find a proper location of the origin so that the internal region includes the source space and the external region includes the locations of interference sources, as depicted in Figure1(a). Therefore, the SSS method can effectively suppress the interference with no signal distortion. For flat sensor arrays, however, the internal region covers only a part of the source space, which extends into the intermediate region, as depicted in Figure 1(b). Therefore, the signal vector has both external and internal components, and, as a result, the signal vector is distorted through the SSS application.Let us first check this fact. An array of sensors and the coordinate system used in our computer simulation are shown in Figure2(a). The sensor array consists of 64 sensors arranged in an 8 × 8 configuration on the plane z=10cm; the plane on which sensors are arranged is called the sensor plane. The sensor array covers an area of 20 cm × 20 cm, and the sensors measure only the z component of the magnetic field, which is the component normal to the sensor plane. This sensor arrangement is almost the same as the one in the MC-6400MCG system (Hitachi High-Technologies Corporation, Tokyo, Japan), which is used in a number of investigations [11–13]. A similar sensor arrangement is used in the KRISS 64-channel biomagnetometer (Bio-Signal Research Center, KRISS, Daejeon, Korea). Therefore, the arrangement of the sensor array assumed in this computer simulation can be considered as a typical one used in current clinical MCG studies.Figure 2
(a) An array of sensors and the coordinate system used in our computer simulation. The sensor array consists of 64 sensors arranged in an 8 × 8 configuration on the planez=10cm, called the sensor plane. The sensor array covers an area of 20 cm × 20 cm, and the sensors measure the z component of the magnetic field, which is the component normal to the sensor plane. The source space (−10≤x≤10cm, −10≤y≤10cm, −7≤z≤7cm) is shown. (b) The source time course assumed for a source located at (−3 cm, 0 cm, and 2 cm). (c) The signal sensor data bSt are shown in the upper panel, and the signal plus sensor-noise data bSt+ε are shown in the lower panel. Here, sensor noise was added to the signal data with the signal-to-noise ratio (SNR) of 10. The sensor time courses are normalized to the maximum value in each panel, and the ordinate indicates the normalized values.
(a)
(b)
(c)A single source was assumed to exist at (−3 cm, 0 cm, and 2 cm), and the source time course assigned to this source is shown in Figure2(b). The signal sensor data bS(t) were computed by projecting the source time course through the sensor lead field obtained using the Biot-Savart law. The signal data bS(t) are shown in the upper panel of Figure 2(c). Sensor noise was then added to the signal data with a signal-to-noise ratio (SNR) of 10 to generate the signal plus sensor-noise data, bSt+ε, which are shown in the lower panel of Figure 2(c). Here, the sensor noise was assumed to be the white Gaussian noise uncorrelated between different sensor channels, that is, the noise vector had a statistical property of ε∼Nε∣0,σ2I, where σ2 is the variance of the noise in all sensor channels.In order to generate environmental noise, four interference sources were assumed to exist; their coordinates and distances from the center of the sensor array are shown in Table1. The locations of these interference sources with respect to the sensor array are shown in Figure 3(a). Four random time courses shown in Figure 3(b) were assigned to these four interference sources, and the interference data bIt were computed. The generated interference data are shown in Figure 3(c).Table 1
Locations of interference sources assumed in a computer simulation.
Source number
Location (cm)
Distance from the sensor array (cm)
1
(−60, 130, 200)
238
2
(300, −150, 360)
485
3
(−72, −80, −200)
236
4
(−105, −150, −100)
214Figure 3
(a) Locations of four interference sources shown with respect to the sensors. (b) Random (normalized) time courses assigned to the four interference sources. (c) Generated interference sensor databIt. These interference sensor data are computed by projecting the interference-source time courses in (b) through the sensor lead field obtained using the Biot-Savart law. The sensor time courses are normalized, and the ordinate indicates the normalized values.
(a)
(b)
(c)We applied the SSS extractors tobSt and bIt in order to check how large the internal and external components are in each of bSt and bIt. The upper and lower panels, respectively, of Figure 4(a) show PintbSt and PextbSt. Here, PintbSt and PextbSt indicate the internal and external components contained in bSt. These results show that the signal vector bSt contains a considerable amount of external components, confirming the validity of our analysis in Section 2.5. The upper and lower panels, respectively, of Figure 4(b) show PintbIt and PextbIt. Here, PintbIt is nearly equal to zero, and PextbIt is almost the same as bIt. These results verify our arguments in Section 2.5 that bIt contains almost no internal components.Figure 4
Results of experiments that apply the SSS extractors tobSt and bIt computed assuming the flat sensor array in Figure 2(a). (a) The upper and lower panels, respectively, show PintbSt and PextbSt, which indicate the internal and external components in bSt. (b) The upper and lower panels, respectively, show PintbIt and PextbIt, which indicate the internal and external components in bIt. In the sensor time courses in (a), they are normalized to the maximum value from the upper panel, and the sensor time courses in (b), they are normalized to the maximum value from the lower panel.
(a)
(b)We performed the same experiments using an MEG helmet-type sensor array for comparison; the arrangement of the helmet sensors used in this computer simulation is shown in Figure5(a). The upper and lower panels, respectively, of Figure 5(b) show PintbSt and PextbSt, and the upper and lower panels, respectively, of Figure 5(c) show PintbIt and PextbIt. These results clearly confirm that the amount of the external components in bSt, as well as the amount of the internal components in bIt, is very small, explaining why the SSS method works well for data from helmet-type sensor arrays used in MEG.Figure 5
Results of the same experiments in which the SSS extractors are applied tobSt and bIt, assuming a helmet-type sensor array used in MEG. (a) Locations of sensors in the helmet-type array assumed in this computer simulation. The sensor arrangement is from the 275-channel whole head sensor array of the Omega™ (VMS Medtech, Coquitlam, Canada). (b) The upper and lower panels, respectively, show PintbSt and PextbSt. (c) The upper and lower panels, respectively, show PintbIt and PextbIt. In the sensor time courses in (b), they are normalized to the maximum value from the upper panel, and the sensor time courses in (c), they are normalized to the maximum value from the lower panel.
(a)
(b)
(c)


### 3.2. Optimal Location of the Origin

We next explored the optimal location of the origin for SSS application to flat sensor data. The origin location significantly affects the final results of the SSS interference suppression, and, thus it is one of the most important parameters in the SSS implementation. In order to see how the origin location affects the SSS results, the internal componentPintbSt and the external component PextbSt were computed with the origin set at four different locations of (0, 0, and zori) where zori was equal to 9 cm, 6 cm, 3 cm, and 0 cm. (Note that the center of the sensor array is located at (0 cm, 0 cm, and 10 cm)). Here, the signal sensor data (plus sensor noise) bSt+ε shown in Figure 2(c) were used. The truncation values LC and LD were, respectively, set at 7 and 3, which are the values found to be optimal in previous investigations [8, 9]. Results of this experiment are shown in Figure 6, exhibiting a general tendency that the signal leakage becomes larger when the origin becomes closer to the sensor plane (z = 10 cm). However, the results also show that when the origin becomes farther from the sensor plane, the SSS results become significantly noisy. In summary, the location of the origin affects both the amount of signal leakage and the noise in the SSS results.Figure 6
The internal componentPintbSt+ε and the external component PextbSt+ε of the signal vector (plus sensor noise) with four different locations of the origin (0, 0, and zori): (a) zori = 0 cm, (b) zori = 3 cm, (c) zori = 6 cm, and (d) zori = 9 cm. The upper panel indicates the internal components, and the lower indicates the external components PextbSt. The filled circle shows the location of the origin, and the square indicates the location of the source. In each pair of the sensor time courses, they are normalized to the maximum value from the upper panel, and the ordinate indicates the normalized values.The amount of signal leakage can be assessed by the signal gainGS defined in (39). Let us derive a quantitative relationship of the signal gain GS versus zori. To do this, voxels with 0.5 cm intervals were assumed in a 3-dimensional source space (−10 ≤ x ≤ 10 cm, −10 ≤ y ≤ 10 cm, and −7 ≤ z ≤ 7 cm), which is shown in Figure 2(a). The signal sensor data bSt were computed by setting the signal source at one of the voxel locations, and the signal gain GS was computed using bSt. The mean signal gain was computed by averaging the signal gains obtained from all voxel locations. Here, the SSS extractor was derived with zori varied from 0 cm to 10 cm. The mean signal gain versus zori is plotted with a broken line in Figure 7(a). It can be seen in these plots that GS gradually decreases as the origin becomes closer to the sensor plane.Figure 7
(a) Plots of the mean signal gainGS, the noise gain Gε, and their ratio GS/Gε versus the z coordinate of the origin zori. The signal sensor data bSt were computed by setting the signal source at one of the voxel locations, and the signal gain GS was obtained. The mean signal gain was computed by averaging GS obtained from all voxel locations. (b) The mean interference gain GI versus zori. The plot with the broken line shows the case rI=20m, and the plot with the solid line shows the case rI=5m. The interference data bI were computed by setting a source at 100 equally spaced locations on the surface of a sphere with a radius of rI, and the interference gain GI was obtained. The mean interference gain was computed by averaging GI obtained from all 100 source locations. The sphere with a radius rI is shown in the upper part.
(a)
(b)The noise gainGε is also plotted in Figure 7(a), and we can see that the noise gain becomes significantly larger as the origin becomes farther from the sensor plane, confirming the general tendency observed in Figure 6. The ratio GS/Gε is also plotted in Figure 7(a), and it can be seen that the ratio gradually increases as the origin becomes closer to the sensor plane. It reaches maximum at zori approximately equal to 9 cm.We next performed computer simulation to derive the relationship between the interference gainGI and the origin location zori. To do this, we assumed a sphere with its radius equal to rI and its center at the center of the sensor array; such a sphere is shown in the upper part of Figure 7(b). The surface of the sphere was defined as a region where interference sources exist, and the interference data bIt were computed by assuming that a single interference source exists on the surface. The interference gain GI was computed with 100 different locations of the interference source, and the mean interference gain is computed and plotted with respect to the origin parameter zori. Here, denoting the polar coordinate of the interference source by (rI,θ,ϕ), 100 locations of the interference source were determined as locations with 10 equal-interval θ and 10 equal-interval φ.The plots ofGI with respect to zori are shown in Figure 7(b). Here, two cases, rI=5m and rI=20m, are shown. The case rI=5m represents cases in which an interference source is located relatively near to the sensors, and the case rI=20m represents cases in which an interference source is located relatively far from the sensors. These plots indicate that the interference gain generally decreases as the origin becomes closer to the sensors, and it reaches a minimum at zori ≈ 9 cm for both cases.


### 3.3. Choices ofLC and LD, Truncated Values of Multipole-Order ℓ

We explore the optimal values of the parametersLC and LD, which are the truncation values of the multipole-order ℓ introduced in (10). So far, these parameters have been set at the values found to be optimal in previous investigations [8, 9], and so, LC and LD were, respectively, set at 7 and 3 in the preceding subsections. In Figure 8, the signal gain GS, noise gain Gε, and their ratio GS/Gε are plotted with respect to zori for three different values of LC (LC=5,6,and7) and two different values of LD (LD=2and3). The noise gain is plotted in Figure 8(a), the signal gain in Figure 8(b), and the gain ratio in Figure 8(c). In these figures, the broken lines indicate the results with LD=2, and the solid lines indicate those with LD=3.Figure 8
(a) The signal gainGS, (b) noise gain Gε, and (c) their ratio GS/Gε plotted with respect to the origin’s z coordinate, zori, for three values of LC (LC=5,6,and7) and two values of LD (LD=2and3). The broken lines indicate the results with LD=2, and the solid lines indicate those with LD=3.
(a)
(b)
(c)It can be seen that the noise gain is considerably smaller whenLD=2 than when LD=3. The signal gain when LD=2 is greater than when LD=3. The differences in the signal and noise gains among different values of LC are generally small. These results suggest that LD should be set at 2, rather than 3. In Figure 8(c), the plots of the gain ratio GS/Gε for LD=2 are shown to have peak maxima when zori > 8, regardless of the value of LC. The results in Figure 8 suggest that, as far as the signal and noise gains concerned, LD=2 should be the best choice, and any value of LC=5,6,and7 may be chosen if we set zori such that zori > 8.The interference gain is plotted with respect tozori for the three values of LC and the two values of LD in Figure 9. Here, the cases rI=5m and rI=20m are, respectively, shown in Figures 9(a) and 9(b), where the broken lines indicate the results of LD=2 and the solid lines indicate those of LD=3. These plots show that regardless of the values of LC and LD, there is a general tendency that the gain decreases as the origin becomes closer to the sensors, and it reaches a minimum near zori equal to 9 cm. Namely, the value of zori ≈ 9 cm gives the minimum interference gain (the maximum shield factor) for all values of LC and LD. On the basis of the results in Figures 8 and 9, we can conclude that the origin parameter zori should be set at 9, that is, the best origin location is determined as (0, 0, and 9), which is 1 cm below the center of the sensor array. We use this value throughout the experiments described below.Figure 9
The interference gain plotted with respect to the origin’sz coordinate, zori, for three values of LC (LC=5,6,and7) and two values of LD (LD=2and3). (a) The distance to the interference source rI is set at 5 m. (b) The distance to the interference source rI is set at 20 m. The broken lines indicate the results with LD=2, and the solid lines indicate those with LD=3.
(a)
(b)


### 3.4. Influence of Sensor Calibration Errors

The SSS method is known to be very sensitive to sensor calibration errors, and an accurately calibrated sensor array is needed for effective suppression of interference [18]. Sensor calibration errors are known to severely affect the shielding capability, so the influence of sensor calibration errors on the interference gain GI is investigated by using erroneous sensor locations and orientations when computing the SSS basis vectors. Here, the sensor location error Eloc is assessed using
(41)Eloc=1M∑j=1Mrjd−r^jdr^jd,where rjd is the true location of the jth sensor, and r^jd is the calibrated location of this sensor. Note that r^jd is assumed to contain an error. Here, rjd is obtained by adding (small) random displacements to r^jd. The signal vector bSt and the interference vector bIt were computed using the true location rjd, while the SSS basis vectors were computed using the calibrated location r^jd. In exactly the same manner, the error in the sensor orientation, Eori, is assessed using
(42)Eori=1M∑j=1Mζj−ζ^jζ^j,where ζj is the true normal vector of the jth sensor, and ζbj is the calibrated normal vector, which is ζ^j=0,0,1. Here, ζj is obtained by adding random vectors to ζ^j. The signal vector bSt and the interference vector bIt were computed using ζj, while the SSS basis vectors were computed using the erroneous orientation ζ^j.In Figure10, the interference gain is plotted with respect to rI, the distance to the interference source, for the three values of LC and the two values of LD. Again, the broken lines indicate the results with LD=2, and the solid lines indicate those with LD=3. Here, the plots in Figure 10(a) indicate the case of no calibration error (Eloc=Eori=0). The plots in Figures 10(b)–10(d), respectively, indicate the cases Eloc=Eori=0.03%, Eloc=Eori=0.1%, and Eloc=Eori=1%. In Figures 10(b)–10(d), the mean values obtained over 100 Monte Carlo trials are plotted. That is, a set of random values were assigned as the errors of sensor locations and orientations in each trial, and the mean results from 100 such trials are plotted.Figure 10
The plots of the interference gainGI with respect to the distance to the interference source, rI, for the three values of LC and the two values of LD. (a) No calibration errors (ELoc=EOri=0). (b) Calibration errors of ELoc=EOri=0.03%. (c) Calibration errors of ELoc=EOri=0.1%. (d) Calibration errors of ELoc=EOri=1%. The broken lines indicate the results with LD=2 and the solid lines indicate those with LD=3. The origin was set at (0 cm, 0 cm, and 9 cm) (zori=9cm).
(a)
(b)
(c)
(d)First, we can observe that the interference gain (i.e., the shield factor) is seriously affected by the calibration errors. When there are no calibration errors, the shield factor 1/GI for rI > 15 m is greater than 104 with LD=2 but it drops to less than 30 when Eloc=Eori=1%. Regarding the optimal choices of the parameters LC and LD, the choice of LD=2 mostly gives greater shield factor than LD=3. There are no significant differences among the choices of LC, but the choice of LC=6 gives slightly better results when Eloc=Eori>0.1%. Since such errors as Eloc=Eori>0.1% may be considered the practical values of sensor calibration errors, the choices of LC=6 and LD=2 are determined as the best choices for the truncation of the multipole parameter ℓ. Note that, according to the arguments for signal and noise gains in Section 3.3, the choices of LC=6 and LD=2 also give the best results.With the choices ofLC=6, LD=2, and zori=9 cm, the interference gain GI is replotted with respect to rI, the distance to the interference source. The results are shown in Figure 11(a). It can be seen that the shield factor (1/GI) exceeds 104 for rI>15m when no calibration errors exist but that the shield factor drops to 100 when calibration errors of 0.1% exist, and to 30 when calibration errors of 1% exist. With the parameter settings, LC=6, LD=2, and 1% sensor calibration errors, signal and noise gains versus zori are plotted in Figure 11(b). Here, the plot with the solid line indicates the case of no calibration errors, and the plot with the broken line indicates the case of 1% calibration errors. The two plots nearly overlap, and this fact confirms that the influence of calibration errors on the signal and noise gains is small.Figure 11
(a) The interference gainGI replotted with respect to rI, the distance to interference sources. The choices of LC=6,LD=2, and zori=9cm were used. (b) The signal and noise gains versus the origin’s z coordinate (zori). The plots with the solid line indicate the case of no calibration errors, and the plots with the broken line indicate the case of 1% calibration errors. The multiple-order truncation was set such that LC=6 and LD=2.
(a)
(b)


### 3.5. The Performance of the SSS Method for Different Sensor Arrays

Here, the SSS performance with two different types of flat sensor arrays is tested: one is a sensor array with a larger sensor coverage and the other is a sensor array consisting of vector sensors. In the sensor array with a larger coverage, the sensors are arranged on a 24 cm × 24 cm coverage area and consist of 10 × 10 sensors that measure the magnetic field normal to the sensor plane. This sensor array is a larger-scale version of the sensor array used in the preceding subsections. The vector sensor array consists of 6 × 6 vector sensors arranged on a 20 cm × 20 cm coverage area. Here, a vector sensor indicates a set of three sensors; each measures one of thex, y, and z components of the magnetic field, and therefore this sensor array actually has a total of 108 sensors.Assuming the sensor array with a larger sensor coverage, the interference gainGI versus rI (the distance to the interference sources) was plotted in Figure 12(a). It can be seen that the plots in this figure are almost the same as the plots in Figure 11(a), suggesting that the influence of the number of sensors and sensor coverage on the shielding capability is rather small. The signal and noise gains with respect to the origin’s z coordinate (zori) were plotted in Figure 12(b), where the plot with the solid line indicates the case of no calibration errors and the plot with the broken line indicates the case of 1% calibration errors. Again, we can observe that the plots here are very similar to the plots in Figure 11(b).Figure 12
The SSS performance computed by assuming a larger sensor array in which the sensors are arranged on a 24 cm × 24 cm coverage area and consist of 10 × 10 sensors that measures the magnetic field normal to the sensor plane. (a) The interference gainGI plotted with respect to rI for six different calibration errors. The choices of LC=6,LD=2, and zori = 9 cm were used. (b) The signal and noise gains versus the origin’s z coordinate (zori). The plots with the solid line indicate the case of no calibration errors, and the plots with the broken line indicate the case of 1% calibration errors. The multiple-order truncation was set such that LC=6 and LD=2.
(a)
(b)The SSS performance for the vector sensor array was investigated. The interference gainGI versus rI was plotted in Figure 13(a). We can observe significant improvements in the shielding capability.Figure 13
The SSS performance computed by assuming a vector sensor array in which the array consists of 6 × 6 vector sensors arranged on a 20 cm × 20 cm coverage area. (a) The interference gainGI plotted with respect to rI for six different calibration errors. The choices of LC=6, LD=2, and zori=9cm were used. (b) The signal and noise gains versus the origin’s z coordinate (zori). The plots with the solid line indicates the case of no calibration errors, and the plots with the broken line indicates the case of 1% calibration errors. The multiple-order truncation was set such that LC=6 and LD=2.
(a)
(b)That is, with 1% calibration errors, the vector sensor array has a shield factor of approximately 500 (GI=2×10−3). The two sensor arrays with normal-component-only sensors have shield factors of nearly 30 with 1% calibration errors, according to Figures 11(a) and 12(a). These results are consistent with the previous investigation [19], which have reported SSS performance improvements due to the use of tangential sensors. The signal and noise gains with respect to the origin’s z coordinate (zori) are plotted for the vector sensor array in Figure 13(b). It can be seen that the plots here are very similar to the plots in Figure 11(b) or 12(b), indicating that these gains are not affected by the use of vector sensors.


### 3.6. Correction of Signal Distortion in Source Space Analysis

#### 3.6.1. Two-Dimensional Current-Density Reconstruction Experiments

As discussed in Section2.5, the SSS-processed signal vector b~Stb~St=PintbSt becomes distorted because the SSS application causes the removal of the external components from the signal vector bSt. This distortion, however, can be corrected by using the SSS-modified lead field because the distorted signal vector is expressed as a sum of distorted lead field vectors, as shown in (37). A series of computer simulations were performed to verify this idea. The results for 2-dimensional current-density mapping are presented in Figure 14.Figure 14
(a) Upper panel: the signal plus sensor-noise time coursesbSt+ε. Lower panel: the interference time courses bIt. (b) Upper panel: the interference-overlapped sensor time courses yt. Lower panel: the SSS interference removal results Pintyt. These time courses are normalized to the maximum value in each panel. (c) The current-density map obtained from bSt+ε. (d) The current-density map obtained from interference-overlapped sensor time courses yt. (e) The current-density map obtained from the SSS interference-removal results Pintyt with the original lead field Lr. (f) The current-density map obtained using Pintyt with the SSS-modified lead field L~r=PintLr. Here, two-dimensional current-density maps on the plane z=2cm, which is the plane 8 cm below the sensor plane, are reconstructed using the field map at t=1200. The signal data bSt were computed assuming two sources, located at (−5 cm, 4 cm, and 3 cm) and (4 cm, −3 cm, and 3 cm).
(a)
(b)
(c)
(d)
(e)
(f)The signal databSt were computed assuming the two sources located at (−5 cm, 4 cm, 3 cm) and (4 cm, −3 cm, and 3 cm). The signal plus sensor-noise data bSt+ε are shown in the upper panel of Figure 14(a). The interference bI(t) is generated by assuming the same four sources as in Figure 3(a); the generated interference is shown in the lower panel of Figure 14(a). The upper panel of Figure 14(b) shows the sensor time courses yt; yt=bSt+ε+ξbIt, where a positive constant ξ controls the signal-to-interference ratio (SIR) (The SIR is defined as the ratio bSt/ξbIt in this computer simulation), which was set equal to 0.25 in this computer simulation. The SSS interference suppression results Pintyt are shown in the lower panel of Figure 14(b).The current-density reconstruction was performed using RENS beamformer, proposed in [20, 21]. The current-density map obtained from the signal plus sensor-noise data, bSt+ε, is shown in Figure 14(c). This map works as the ground truth for the following experiments. The current-density map obtained from the interference-overlapped data yt is shown in Figure 14(d). The results far from the ground truth are obtained due to the overlap of the large interference.The current-density reconstruction was performed twice, using the SSS resultsPintyt with the original lead field Lr and with the SSS-modified lead field L~r: L~r=PintLr. The current-density map obtained with the original lead field Lr is shown in Figure 14(e). A considerable amount of distortion can be seen in this current-density map, although the SSS method seems to have nearly perfectly removed the interference according to the SSS-processed time courses (the lower panel of Figure 14(b)). The current-density map obtained with the SSS-modified lead field is shown in Figure 14(f). Here, results very close to the ground truth can be obtained; these results verify the idea that the signal distortion can be compensated for by using the SSS-modified lead field in the voxel space analysis.


#### 3.6.2. Three-Dimensional Source Localization Experiments

Next, three-dimensional source localization experiments were performed. Reconstruction results obtained using the signal plus sensor-noise data shown in Figure2(c) are shown in Figure 15(a). A single source is reconstructed near (−3, 0, and 2), which is the location assumed for the data generation. These results work as the ground truth when evaluating the following results. The interference-overlapped sensor data yt were computed using the interference data shown in Figure 3(c) overlapped onto these signal plus sensor-noise data with the SIR equal to 0.25. The source reconstruction was performed using the interference-overlapped sensor data yt, and the results are shown in Figure 15(b). A large influence from the interference can be seen here.Figure 15
(a) Source reconstruction results obtained using the signal plus sensor-noise time courses in Figure2(c). (b) Source reconstruction results obtained using the sensor time courses yt created by adding the interference data shown in Figure 3(c) onto the signal plus sensor-noise time courses in Figure 2(c). (c) Source reconstruction results obtained using the SSS interference removal results Pintyt with the original lead field Lr. (d) Source reconstruction results obtained using the SSS interference removal results Pintyt with the SSS-modified lead field L~r. The RENS beamformer [20, 21] was applied to the field data at t=1260 for three-dimensional source reconstruction.
(a)
(b)
(c)
(d)The source reconstruction was carried out using the SSS interference suppression resultsPintyt with the original lead field Lr. The results are shown in Figure 15(c), in which a single source is reconstructed but the location of the source differs considerably from the assumed location. The source reconstruction was then carried out using the SSS-modified lead field L~r, and the results are shown in Figure 15(d). These results are almost the same as the ground truth, verifying the effectiveness of the use of the SSS-modified lead field in the voxel space analysis.


## 3.1. Problems with Flat Sensor Data for SSS Applications

In order to show that the SSS method is still effective for data from sensors arranged on a flat plane, a series of computer simulation has been performed. First, we clarify the problems caused when the SSS method is applied to flat sensor data by comparing two cases of SSS application: one in which a helmet-type sensor array is used and the other in which a flat sensor array is used. For helmet-type sensor arrays, the key condition for the SSS method can be fulfilled, that is, one can find a proper location of the origin so that the internal region includes the source space and the external region includes the locations of interference sources, as depicted in Figure1(a). Therefore, the SSS method can effectively suppress the interference with no signal distortion. For flat sensor arrays, however, the internal region covers only a part of the source space, which extends into the intermediate region, as depicted in Figure 1(b). Therefore, the signal vector has both external and internal components, and, as a result, the signal vector is distorted through the SSS application.Let us first check this fact. An array of sensors and the coordinate system used in our computer simulation are shown in Figure2(a). The sensor array consists of 64 sensors arranged in an 8 × 8 configuration on the plane z=10cm; the plane on which sensors are arranged is called the sensor plane. The sensor array covers an area of 20 cm × 20 cm, and the sensors measure only the z component of the magnetic field, which is the component normal to the sensor plane. This sensor arrangement is almost the same as the one in the MC-6400MCG system (Hitachi High-Technologies Corporation, Tokyo, Japan), which is used in a number of investigations [11–13]. A similar sensor arrangement is used in the KRISS 64-channel biomagnetometer (Bio-Signal Research Center, KRISS, Daejeon, Korea). Therefore, the arrangement of the sensor array assumed in this computer simulation can be considered as a typical one used in current clinical MCG studies.Figure 2
(a) An array of sensors and the coordinate system used in our computer simulation. The sensor array consists of 64 sensors arranged in an 8 × 8 configuration on the planez=10cm, called the sensor plane. The sensor array covers an area of 20 cm × 20 cm, and the sensors measure the z component of the magnetic field, which is the component normal to the sensor plane. The source space (−10≤x≤10cm, −10≤y≤10cm, −7≤z≤7cm) is shown. (b) The source time course assumed for a source located at (−3 cm, 0 cm, and 2 cm). (c) The signal sensor data bSt are shown in the upper panel, and the signal plus sensor-noise data bSt+ε are shown in the lower panel. Here, sensor noise was added to the signal data with the signal-to-noise ratio (SNR) of 10. The sensor time courses are normalized to the maximum value in each panel, and the ordinate indicates the normalized values.
(a)
(b)
(c)A single source was assumed to exist at (−3 cm, 0 cm, and 2 cm), and the source time course assigned to this source is shown in Figure2(b). The signal sensor data bS(t) were computed by projecting the source time course through the sensor lead field obtained using the Biot-Savart law. The signal data bS(t) are shown in the upper panel of Figure 2(c). Sensor noise was then added to the signal data with a signal-to-noise ratio (SNR) of 10 to generate the signal plus sensor-noise data, bSt+ε, which are shown in the lower panel of Figure 2(c). Here, the sensor noise was assumed to be the white Gaussian noise uncorrelated between different sensor channels, that is, the noise vector had a statistical property of ε∼Nε∣0,σ2I, where σ2 is the variance of the noise in all sensor channels.In order to generate environmental noise, four interference sources were assumed to exist; their coordinates and distances from the center of the sensor array are shown in Table1. The locations of these interference sources with respect to the sensor array are shown in Figure 3(a). Four random time courses shown in Figure 3(b) were assigned to these four interference sources, and the interference data bIt were computed. The generated interference data are shown in Figure 3(c).Table 1
Locations of interference sources assumed in a computer simulation.
Source number
Location (cm)
Distance from the sensor array (cm)
1
(−60, 130, 200)
238
2
(300, −150, 360)
485
3
(−72, −80, −200)
236
4
(−105, −150, −100)
214Figure 3
(a) Locations of four interference sources shown with respect to the sensors. (b) Random (normalized) time courses assigned to the four interference sources. (c) Generated interference sensor databIt. These interference sensor data are computed by projecting the interference-source time courses in (b) through the sensor lead field obtained using the Biot-Savart law. The sensor time courses are normalized, and the ordinate indicates the normalized values.
(a)
(b)
(c)We applied the SSS extractors tobSt and bIt in order to check how large the internal and external components are in each of bSt and bIt. The upper and lower panels, respectively, of Figure 4(a) show PintbSt and PextbSt. Here, PintbSt and PextbSt indicate the internal and external components contained in bSt. These results show that the signal vector bSt contains a considerable amount of external components, confirming the validity of our analysis in Section 2.5. The upper and lower panels, respectively, of Figure 4(b) show PintbIt and PextbIt. Here, PintbIt is nearly equal to zero, and PextbIt is almost the same as bIt. These results verify our arguments in Section 2.5 that bIt contains almost no internal components.Figure 4
Results of experiments that apply the SSS extractors tobSt and bIt computed assuming the flat sensor array in Figure 2(a). (a) The upper and lower panels, respectively, show PintbSt and PextbSt, which indicate the internal and external components in bSt. (b) The upper and lower panels, respectively, show PintbIt and PextbIt, which indicate the internal and external components in bIt. In the sensor time courses in (a), they are normalized to the maximum value from the upper panel, and the sensor time courses in (b), they are normalized to the maximum value from the lower panel.
(a)
(b)We performed the same experiments using an MEG helmet-type sensor array for comparison; the arrangement of the helmet sensors used in this computer simulation is shown in Figure5(a). The upper and lower panels, respectively, of Figure 5(b) show PintbSt and PextbSt, and the upper and lower panels, respectively, of Figure 5(c) show PintbIt and PextbIt. These results clearly confirm that the amount of the external components in bSt, as well as the amount of the internal components in bIt, is very small, explaining why the SSS method works well for data from helmet-type sensor arrays used in MEG.Figure 5
Results of the same experiments in which the SSS extractors are applied tobSt and bIt, assuming a helmet-type sensor array used in MEG. (a) Locations of sensors in the helmet-type array assumed in this computer simulation. The sensor arrangement is from the 275-channel whole head sensor array of the Omega™ (VMS Medtech, Coquitlam, Canada). (b) The upper and lower panels, respectively, show PintbSt and PextbSt. (c) The upper and lower panels, respectively, show PintbIt and PextbIt. In the sensor time courses in (b), they are normalized to the maximum value from the upper panel, and the sensor time courses in (c), they are normalized to the maximum value from the lower panel.
(a)
(b)
(c)


## 3.2. Optimal Location of the Origin

We next explored the optimal location of the origin for SSS application to flat sensor data. The origin location significantly affects the final results of the SSS interference suppression, and, thus it is one of the most important parameters in the SSS implementation. In order to see how the origin location affects the SSS results, the internal componentPintbSt and the external component PextbSt were computed with the origin set at four different locations of (0, 0, and zori) where zori was equal to 9 cm, 6 cm, 3 cm, and 0 cm. (Note that the center of the sensor array is located at (0 cm, 0 cm, and 10 cm)). Here, the signal sensor data (plus sensor noise) bSt+ε shown in Figure 2(c) were used. The truncation values LC and LD were, respectively, set at 7 and 3, which are the values found to be optimal in previous investigations [8, 9]. Results of this experiment are shown in Figure 6, exhibiting a general tendency that the signal leakage becomes larger when the origin becomes closer to the sensor plane (z = 10 cm). However, the results also show that when the origin becomes farther from the sensor plane, the SSS results become significantly noisy. In summary, the location of the origin affects both the amount of signal leakage and the noise in the SSS results.Figure 6
The internal componentPintbSt+ε and the external component PextbSt+ε of the signal vector (plus sensor noise) with four different locations of the origin (0, 0, and zori): (a) zori = 0 cm, (b) zori = 3 cm, (c) zori = 6 cm, and (d) zori = 9 cm. The upper panel indicates the internal components, and the lower indicates the external components PextbSt. The filled circle shows the location of the origin, and the square indicates the location of the source. In each pair of the sensor time courses, they are normalized to the maximum value from the upper panel, and the ordinate indicates the normalized values.The amount of signal leakage can be assessed by the signal gainGS defined in (39). Let us derive a quantitative relationship of the signal gain GS versus zori. To do this, voxels with 0.5 cm intervals were assumed in a 3-dimensional source space (−10 ≤ x ≤ 10 cm, −10 ≤ y ≤ 10 cm, and −7 ≤ z ≤ 7 cm), which is shown in Figure 2(a). The signal sensor data bSt were computed by setting the signal source at one of the voxel locations, and the signal gain GS was computed using bSt. The mean signal gain was computed by averaging the signal gains obtained from all voxel locations. Here, the SSS extractor was derived with zori varied from 0 cm to 10 cm. The mean signal gain versus zori is plotted with a broken line in Figure 7(a). It can be seen in these plots that GS gradually decreases as the origin becomes closer to the sensor plane.Figure 7
(a) Plots of the mean signal gainGS, the noise gain Gε, and their ratio GS/Gε versus the z coordinate of the origin zori. The signal sensor data bSt were computed by setting the signal source at one of the voxel locations, and the signal gain GS was obtained. The mean signal gain was computed by averaging GS obtained from all voxel locations. (b) The mean interference gain GI versus zori. The plot with the broken line shows the case rI=20m, and the plot with the solid line shows the case rI=5m. The interference data bI were computed by setting a source at 100 equally spaced locations on the surface of a sphere with a radius of rI, and the interference gain GI was obtained. The mean interference gain was computed by averaging GI obtained from all 100 source locations. The sphere with a radius rI is shown in the upper part.
(a)
(b)The noise gainGε is also plotted in Figure 7(a), and we can see that the noise gain becomes significantly larger as the origin becomes farther from the sensor plane, confirming the general tendency observed in Figure 6. The ratio GS/Gε is also plotted in Figure 7(a), and it can be seen that the ratio gradually increases as the origin becomes closer to the sensor plane. It reaches maximum at zori approximately equal to 9 cm.We next performed computer simulation to derive the relationship between the interference gainGI and the origin location zori. To do this, we assumed a sphere with its radius equal to rI and its center at the center of the sensor array; such a sphere is shown in the upper part of Figure 7(b). The surface of the sphere was defined as a region where interference sources exist, and the interference data bIt were computed by assuming that a single interference source exists on the surface. The interference gain GI was computed with 100 different locations of the interference source, and the mean interference gain is computed and plotted with respect to the origin parameter zori. Here, denoting the polar coordinate of the interference source by (rI,θ,ϕ), 100 locations of the interference source were determined as locations with 10 equal-interval θ and 10 equal-interval φ.The plots ofGI with respect to zori are shown in Figure 7(b). Here, two cases, rI=5m and rI=20m, are shown. The case rI=5m represents cases in which an interference source is located relatively near to the sensors, and the case rI=20m represents cases in which an interference source is located relatively far from the sensors. These plots indicate that the interference gain generally decreases as the origin becomes closer to the sensors, and it reaches a minimum at zori ≈ 9 cm for both cases.


## 3.3. Choices ofLC and LD, Truncated Values of Multipole-Order ℓ

We explore the optimal values of the parametersLC and LD, which are the truncation values of the multipole-order ℓ introduced in (10). So far, these parameters have been set at the values found to be optimal in previous investigations [8, 9], and so, LC and LD were, respectively, set at 7 and 3 in the preceding subsections. In Figure 8, the signal gain GS, noise gain Gε, and their ratio GS/Gε are plotted with respect to zori for three different values of LC (LC=5,6,and7) and two different values of LD (LD=2and3). The noise gain is plotted in Figure 8(a), the signal gain in Figure 8(b), and the gain ratio in Figure 8(c). In these figures, the broken lines indicate the results with LD=2, and the solid lines indicate those with LD=3.Figure 8
(a) The signal gainGS, (b) noise gain Gε, and (c) their ratio GS/Gε plotted with respect to the origin’s z coordinate, zori, for three values of LC (LC=5,6,and7) and two values of LD (LD=2and3). The broken lines indicate the results with LD=2, and the solid lines indicate those with LD=3.
(a)
(b)
(c)It can be seen that the noise gain is considerably smaller whenLD=2 than when LD=3. The signal gain when LD=2 is greater than when LD=3. The differences in the signal and noise gains among different values of LC are generally small. These results suggest that LD should be set at 2, rather than 3. In Figure 8(c), the plots of the gain ratio GS/Gε for LD=2 are shown to have peak maxima when zori > 8, regardless of the value of LC. The results in Figure 8 suggest that, as far as the signal and noise gains concerned, LD=2 should be the best choice, and any value of LC=5,6,and7 may be chosen if we set zori such that zori > 8.The interference gain is plotted with respect tozori for the three values of LC and the two values of LD in Figure 9. Here, the cases rI=5m and rI=20m are, respectively, shown in Figures 9(a) and 9(b), where the broken lines indicate the results of LD=2 and the solid lines indicate those of LD=3. These plots show that regardless of the values of LC and LD, there is a general tendency that the gain decreases as the origin becomes closer to the sensors, and it reaches a minimum near zori equal to 9 cm. Namely, the value of zori ≈ 9 cm gives the minimum interference gain (the maximum shield factor) for all values of LC and LD. On the basis of the results in Figures 8 and 9, we can conclude that the origin parameter zori should be set at 9, that is, the best origin location is determined as (0, 0, and 9), which is 1 cm below the center of the sensor array. We use this value throughout the experiments described below.Figure 9
The interference gain plotted with respect to the origin’sz coordinate, zori, for three values of LC (LC=5,6,and7) and two values of LD (LD=2and3). (a) The distance to the interference source rI is set at 5 m. (b) The distance to the interference source rI is set at 20 m. The broken lines indicate the results with LD=2, and the solid lines indicate those with LD=3.
(a)
(b)


## 3.4. Influence of Sensor Calibration Errors

The SSS method is known to be very sensitive to sensor calibration errors, and an accurately calibrated sensor array is needed for effective suppression of interference [18]. Sensor calibration errors are known to severely affect the shielding capability, so the influence of sensor calibration errors on the interference gain GI is investigated by using erroneous sensor locations and orientations when computing the SSS basis vectors. Here, the sensor location error Eloc is assessed using
(41)Eloc=1M∑j=1Mrjd−r^jdr^jd,where rjd is the true location of the jth sensor, and r^jd is the calibrated location of this sensor. Note that r^jd is assumed to contain an error. Here, rjd is obtained by adding (small) random displacements to r^jd. The signal vector bSt and the interference vector bIt were computed using the true location rjd, while the SSS basis vectors were computed using the calibrated location r^jd. In exactly the same manner, the error in the sensor orientation, Eori, is assessed using
(42)Eori=1M∑j=1Mζj−ζ^jζ^j,where ζj is the true normal vector of the jth sensor, and ζbj is the calibrated normal vector, which is ζ^j=0,0,1. Here, ζj is obtained by adding random vectors to ζ^j. The signal vector bSt and the interference vector bIt were computed using ζj, while the SSS basis vectors were computed using the erroneous orientation ζ^j.In Figure10, the interference gain is plotted with respect to rI, the distance to the interference source, for the three values of LC and the two values of LD. Again, the broken lines indicate the results with LD=2, and the solid lines indicate those with LD=3. Here, the plots in Figure 10(a) indicate the case of no calibration error (Eloc=Eori=0). The plots in Figures 10(b)–10(d), respectively, indicate the cases Eloc=Eori=0.03%, Eloc=Eori=0.1%, and Eloc=Eori=1%. In Figures 10(b)–10(d), the mean values obtained over 100 Monte Carlo trials are plotted. That is, a set of random values were assigned as the errors of sensor locations and orientations in each trial, and the mean results from 100 such trials are plotted.Figure 10
The plots of the interference gainGI with respect to the distance to the interference source, rI, for the three values of LC and the two values of LD. (a) No calibration errors (ELoc=EOri=0). (b) Calibration errors of ELoc=EOri=0.03%. (c) Calibration errors of ELoc=EOri=0.1%. (d) Calibration errors of ELoc=EOri=1%. The broken lines indicate the results with LD=2 and the solid lines indicate those with LD=3. The origin was set at (0 cm, 0 cm, and 9 cm) (zori=9cm).
(a)
(b)
(c)
(d)First, we can observe that the interference gain (i.e., the shield factor) is seriously affected by the calibration errors. When there are no calibration errors, the shield factor 1/GI for rI > 15 m is greater than 104 with LD=2 but it drops to less than 30 when Eloc=Eori=1%. Regarding the optimal choices of the parameters LC and LD, the choice of LD=2 mostly gives greater shield factor than LD=3. There are no significant differences among the choices of LC, but the choice of LC=6 gives slightly better results when Eloc=Eori>0.1%. Since such errors as Eloc=Eori>0.1% may be considered the practical values of sensor calibration errors, the choices of LC=6 and LD=2 are determined as the best choices for the truncation of the multipole parameter ℓ. Note that, according to the arguments for signal and noise gains in Section 3.3, the choices of LC=6 and LD=2 also give the best results.With the choices ofLC=6, LD=2, and zori=9 cm, the interference gain GI is replotted with respect to rI, the distance to the interference source. The results are shown in Figure 11(a). It can be seen that the shield factor (1/GI) exceeds 104 for rI>15m when no calibration errors exist but that the shield factor drops to 100 when calibration errors of 0.1% exist, and to 30 when calibration errors of 1% exist. With the parameter settings, LC=6, LD=2, and 1% sensor calibration errors, signal and noise gains versus zori are plotted in Figure 11(b). Here, the plot with the solid line indicates the case of no calibration errors, and the plot with the broken line indicates the case of 1% calibration errors. The two plots nearly overlap, and this fact confirms that the influence of calibration errors on the signal and noise gains is small.Figure 11
(a) The interference gainGI replotted with respect to rI, the distance to interference sources. The choices of LC=6,LD=2, and zori=9cm were used. (b) The signal and noise gains versus the origin’s z coordinate (zori). The plots with the solid line indicate the case of no calibration errors, and the plots with the broken line indicate the case of 1% calibration errors. The multiple-order truncation was set such that LC=6 and LD=2.
(a)
(b)


## 3.5. The Performance of the SSS Method for Different Sensor Arrays

Here, the SSS performance with two different types of flat sensor arrays is tested: one is a sensor array with a larger sensor coverage and the other is a sensor array consisting of vector sensors. In the sensor array with a larger coverage, the sensors are arranged on a 24 cm × 24 cm coverage area and consist of 10 × 10 sensors that measure the magnetic field normal to the sensor plane. This sensor array is a larger-scale version of the sensor array used in the preceding subsections. The vector sensor array consists of 6 × 6 vector sensors arranged on a 20 cm × 20 cm coverage area. Here, a vector sensor indicates a set of three sensors; each measures one of thex, y, and z components of the magnetic field, and therefore this sensor array actually has a total of 108 sensors.Assuming the sensor array with a larger sensor coverage, the interference gainGI versus rI (the distance to the interference sources) was plotted in Figure 12(a). It can be seen that the plots in this figure are almost the same as the plots in Figure 11(a), suggesting that the influence of the number of sensors and sensor coverage on the shielding capability is rather small. The signal and noise gains with respect to the origin’s z coordinate (zori) were plotted in Figure 12(b), where the plot with the solid line indicates the case of no calibration errors and the plot with the broken line indicates the case of 1% calibration errors. Again, we can observe that the plots here are very similar to the plots in Figure 11(b).Figure 12
The SSS performance computed by assuming a larger sensor array in which the sensors are arranged on a 24 cm × 24 cm coverage area and consist of 10 × 10 sensors that measures the magnetic field normal to the sensor plane. (a) The interference gainGI plotted with respect to rI for six different calibration errors. The choices of LC=6,LD=2, and zori = 9 cm were used. (b) The signal and noise gains versus the origin’s z coordinate (zori). The plots with the solid line indicate the case of no calibration errors, and the plots with the broken line indicate the case of 1% calibration errors. The multiple-order truncation was set such that LC=6 and LD=2.
(a)
(b)The SSS performance for the vector sensor array was investigated. The interference gainGI versus rI was plotted in Figure 13(a). We can observe significant improvements in the shielding capability.Figure 13
The SSS performance computed by assuming a vector sensor array in which the array consists of 6 × 6 vector sensors arranged on a 20 cm × 20 cm coverage area. (a) The interference gainGI plotted with respect to rI for six different calibration errors. The choices of LC=6, LD=2, and zori=9cm were used. (b) The signal and noise gains versus the origin’s z coordinate (zori). The plots with the solid line indicates the case of no calibration errors, and the plots with the broken line indicates the case of 1% calibration errors. The multiple-order truncation was set such that LC=6 and LD=2.
(a)
(b)That is, with 1% calibration errors, the vector sensor array has a shield factor of approximately 500 (GI=2×10−3). The two sensor arrays with normal-component-only sensors have shield factors of nearly 30 with 1% calibration errors, according to Figures 11(a) and 12(a). These results are consistent with the previous investigation [19], which have reported SSS performance improvements due to the use of tangential sensors. The signal and noise gains with respect to the origin’s z coordinate (zori) are plotted for the vector sensor array in Figure 13(b). It can be seen that the plots here are very similar to the plots in Figure 11(b) or 12(b), indicating that these gains are not affected by the use of vector sensors.


## 3.6. Correction of Signal Distortion in Source Space Analysis

### 3.6.1. Two-Dimensional Current-Density Reconstruction Experiments

As discussed in Section2.5, the SSS-processed signal vector b~Stb~St=PintbSt becomes distorted because the SSS application causes the removal of the external components from the signal vector bSt. This distortion, however, can be corrected by using the SSS-modified lead field because the distorted signal vector is expressed as a sum of distorted lead field vectors, as shown in (37). A series of computer simulations were performed to verify this idea. The results for 2-dimensional current-density mapping are presented in Figure 14.Figure 14
(a) Upper panel: the signal plus sensor-noise time coursesbSt+ε. Lower panel: the interference time courses bIt. (b) Upper panel: the interference-overlapped sensor time courses yt. Lower panel: the SSS interference removal results Pintyt. These time courses are normalized to the maximum value in each panel. (c) The current-density map obtained from bSt+ε. (d) The current-density map obtained from interference-overlapped sensor time courses yt. (e) The current-density map obtained from the SSS interference-removal results Pintyt with the original lead field Lr. (f) The current-density map obtained using Pintyt with the SSS-modified lead field L~r=PintLr. Here, two-dimensional current-density maps on the plane z=2cm, which is the plane 8 cm below the sensor plane, are reconstructed using the field map at t=1200. The signal data bSt were computed assuming two sources, located at (−5 cm, 4 cm, and 3 cm) and (4 cm, −3 cm, and 3 cm).
(a)
(b)
(c)
(d)
(e)
(f)The signal databSt were computed assuming the two sources located at (−5 cm, 4 cm, 3 cm) and (4 cm, −3 cm, and 3 cm). The signal plus sensor-noise data bSt+ε are shown in the upper panel of Figure 14(a). The interference bI(t) is generated by assuming the same four sources as in Figure 3(a); the generated interference is shown in the lower panel of Figure 14(a). The upper panel of Figure 14(b) shows the sensor time courses yt; yt=bSt+ε+ξbIt, where a positive constant ξ controls the signal-to-interference ratio (SIR) (The SIR is defined as the ratio bSt/ξbIt in this computer simulation), which was set equal to 0.25 in this computer simulation. The SSS interference suppression results Pintyt are shown in the lower panel of Figure 14(b).The current-density reconstruction was performed using RENS beamformer, proposed in [20, 21]. The current-density map obtained from the signal plus sensor-noise data, bSt+ε, is shown in Figure 14(c). This map works as the ground truth for the following experiments. The current-density map obtained from the interference-overlapped data yt is shown in Figure 14(d). The results far from the ground truth are obtained due to the overlap of the large interference.The current-density reconstruction was performed twice, using the SSS resultsPintyt with the original lead field Lr and with the SSS-modified lead field L~r: L~r=PintLr. The current-density map obtained with the original lead field Lr is shown in Figure 14(e). A considerable amount of distortion can be seen in this current-density map, although the SSS method seems to have nearly perfectly removed the interference according to the SSS-processed time courses (the lower panel of Figure 14(b)). The current-density map obtained with the SSS-modified lead field is shown in Figure 14(f). Here, results very close to the ground truth can be obtained; these results verify the idea that the signal distortion can be compensated for by using the SSS-modified lead field in the voxel space analysis.


### 3.6.2. Three-Dimensional Source Localization Experiments

Next, three-dimensional source localization experiments were performed. Reconstruction results obtained using the signal plus sensor-noise data shown in Figure2(c) are shown in Figure 15(a). A single source is reconstructed near (−3, 0, and 2), which is the location assumed for the data generation. These results work as the ground truth when evaluating the following results. The interference-overlapped sensor data yt were computed using the interference data shown in Figure 3(c) overlapped onto these signal plus sensor-noise data with the SIR equal to 0.25. The source reconstruction was performed using the interference-overlapped sensor data yt, and the results are shown in Figure 15(b). A large influence from the interference can be seen here.Figure 15
(a) Source reconstruction results obtained using the signal plus sensor-noise time courses in Figure2(c). (b) Source reconstruction results obtained using the sensor time courses yt created by adding the interference data shown in Figure 3(c) onto the signal plus sensor-noise time courses in Figure 2(c). (c) Source reconstruction results obtained using the SSS interference removal results Pintyt with the original lead field Lr. (d) Source reconstruction results obtained using the SSS interference removal results Pintyt with the SSS-modified lead field L~r. The RENS beamformer [20, 21] was applied to the field data at t=1260 for three-dimensional source reconstruction.
(a)
(b)
(c)
(d)The source reconstruction was carried out using the SSS interference suppression resultsPintyt with the original lead field Lr. The results are shown in Figure 15(c), in which a single source is reconstructed but the location of the source differs considerably from the assumed location. The source reconstruction was then carried out using the SSS-modified lead field L~r, and the results are shown in Figure 15(d). These results are almost the same as the ground truth, verifying the effectiveness of the use of the SSS-modified lead field in the voxel space analysis.


## 3.6.1. Two-Dimensional Current-Density Reconstruction Experiments

As discussed in Section2.5, the SSS-processed signal vector b~Stb~St=PintbSt becomes distorted because the SSS application causes the removal of the external components from the signal vector bSt. This distortion, however, can be corrected by using the SSS-modified lead field because the distorted signal vector is expressed as a sum of distorted lead field vectors, as shown in (37). A series of computer simulations were performed to verify this idea. The results for 2-dimensional current-density mapping are presented in Figure 14.Figure 14
(a) Upper panel: the signal plus sensor-noise time coursesbSt+ε. Lower panel: the interference time courses bIt. (b) Upper panel: the interference-overlapped sensor time courses yt. Lower panel: the SSS interference removal results Pintyt. These time courses are normalized to the maximum value in each panel. (c) The current-density map obtained from bSt+ε. (d) The current-density map obtained from interference-overlapped sensor time courses yt. (e) The current-density map obtained from the SSS interference-removal results Pintyt with the original lead field Lr. (f) The current-density map obtained using Pintyt with the SSS-modified lead field L~r=PintLr. Here, two-dimensional current-density maps on the plane z=2cm, which is the plane 8 cm below the sensor plane, are reconstructed using the field map at t=1200. The signal data bSt were computed assuming two sources, located at (−5 cm, 4 cm, and 3 cm) and (4 cm, −3 cm, and 3 cm).
(a)
(b)
(c)
(d)
(e)
(f)The signal databSt were computed assuming the two sources located at (−5 cm, 4 cm, 3 cm) and (4 cm, −3 cm, and 3 cm). The signal plus sensor-noise data bSt+ε are shown in the upper panel of Figure 14(a). The interference bI(t) is generated by assuming the same four sources as in Figure 3(a); the generated interference is shown in the lower panel of Figure 14(a). The upper panel of Figure 14(b) shows the sensor time courses yt; yt=bSt+ε+ξbIt, where a positive constant ξ controls the signal-to-interference ratio (SIR) (The SIR is defined as the ratio bSt/ξbIt in this computer simulation), which was set equal to 0.25 in this computer simulation. The SSS interference suppression results Pintyt are shown in the lower panel of Figure 14(b).The current-density reconstruction was performed using RENS beamformer, proposed in [20, 21]. The current-density map obtained from the signal plus sensor-noise data, bSt+ε, is shown in Figure 14(c). This map works as the ground truth for the following experiments. The current-density map obtained from the interference-overlapped data yt is shown in Figure 14(d). The results far from the ground truth are obtained due to the overlap of the large interference.The current-density reconstruction was performed twice, using the SSS resultsPintyt with the original lead field Lr and with the SSS-modified lead field L~r: L~r=PintLr. The current-density map obtained with the original lead field Lr is shown in Figure 14(e). A considerable amount of distortion can be seen in this current-density map, although the SSS method seems to have nearly perfectly removed the interference according to the SSS-processed time courses (the lower panel of Figure 14(b)). The current-density map obtained with the SSS-modified lead field is shown in Figure 14(f). Here, results very close to the ground truth can be obtained; these results verify the idea that the signal distortion can be compensated for by using the SSS-modified lead field in the voxel space analysis.


## 3.6.2. Three-Dimensional Source Localization Experiments

Next, three-dimensional source localization experiments were performed. Reconstruction results obtained using the signal plus sensor-noise data shown in Figure2(c) are shown in Figure 15(a). A single source is reconstructed near (−3, 0, and 2), which is the location assumed for the data generation. These results work as the ground truth when evaluating the following results. The interference-overlapped sensor data yt were computed using the interference data shown in Figure 3(c) overlapped onto these signal plus sensor-noise data with the SIR equal to 0.25. The source reconstruction was performed using the interference-overlapped sensor data yt, and the results are shown in Figure 15(b). A large influence from the interference can be seen here.Figure 15
(a) Source reconstruction results obtained using the signal plus sensor-noise time courses in Figure2(c). (b) Source reconstruction results obtained using the sensor time courses yt created by adding the interference data shown in Figure 3(c) onto the signal plus sensor-noise time courses in Figure 2(c). (c) Source reconstruction results obtained using the SSS interference removal results Pintyt with the original lead field Lr. (d) Source reconstruction results obtained using the SSS interference removal results Pintyt with the SSS-modified lead field L~r. The RENS beamformer [20, 21] was applied to the field data at t=1260 for three-dimensional source reconstruction.
(a)
(b)
(c)
(d)The source reconstruction was carried out using the SSS interference suppression resultsPintyt with the original lead field Lr. The results are shown in Figure 15(c), in which a single source is reconstructed but the location of the source differs considerably from the assumed location. The source reconstruction was then carried out using the SSS-modified lead field L~r, and the results are shown in Figure 15(d). These results are almost the same as the ground truth, verifying the effectiveness of the use of the SSS-modified lead field in the voxel space analysis.


## 4. Discussion and Summary

This paper presented computer simulation-based investigation to explore the possibility of applying the SSS method to data measured from an array of sensors arranged on a flat plane, which is commonly used in magnetocardiographic applications. The findings from the investigation are summarized as follows:(1)
When applying the SSS method to data from a flat sensor array, a signal vector has components of the external subspace as well as those of the internal subspace. As a result, the signal is distorted through the SSS interference suppression process.
(2)
The signal distortion can be compensated for by using the SSS-modified lead field in voxel space analysis. The computer simulations using two-dimensional current-density mapping and three-dimensional source localization confirmed that the distortion can be corrected in the voxel space.
(3)
The origin location can significantly affect the results of the SSS method. It is shown that the optimal location of the origin can be determined by assessing the dependence of signal and noise gains of the SSS extractor on the origin location. The optimal location is empirically found to be approximately 1 cm below the sensor plane for typical flat sensor arrays used in MCG applications.
(4)
The optimal values of the parametersLC and LD, the truncation values of the multipole-order ℓ, can also be determined by evaluating dependences of the signal, noise, and interference gains (shield factor) on these parameters. Results of computer simulation suggest LD=6 and LD=2 to be the optimal choices for typical flat sensor arrays currently used in clinical magnetocardiography.
(5)
The sensor calibration errors affect the shielding capability. The shield factor exceeds 104 for interference originating from fairly distant sources (rI>15m) when no calibration errors exist. However, the shield factor drops to approximately 100 when the calibration errors become 0.1% and to 30 when the calibration errors become 1%.
(6)
The shielding capability can significantly be improved by using vector sensors, which measure thex, y, and z components of the magnetic field. With 1% calibration errors, a vector sensor array still maintains a shield factor of approximately 500, while the arrays with sensors measuring only the normal direction have a shield factor of about 30 with the same 1% calibration errors.Finally, it should be emphasized that the SSS interference suppression method is effective even for arrays of sensors arranged on a flat plane. This is the main finding in this paper. The use of such signal processing methods with low-cost sensors, such as magnetoresistive sensors, can lead to the development of low-initial-cost and maintenance-free magnetocardiography systems, which in the near future may replace the electrocardiogram now routinely used in hospitals.



---

*Source: 7689589-2018-04-26.xml*