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

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 parameters LC 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.


Introduction
Development of a sensor system that can measure biomagnetic signals in room-temperature environments has gained great interest. One promising candidate among roomtemperature sensors for biomagnetic systems is magnetoresistive (MR) sensors [1][2][3][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][9][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 socalled "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][12][13][14][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 Section 2, 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.

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 the mth sensor as y m . The measurement from the whole sensor array is expressed as a column vector y y = y 1 , y 2 , … , y M T . 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 by r r = x, y, z . The source magnitude at r is denoted by a scalar s r . The source vector is denoted by s(r), and the source orientation is denoted by η = η x , η y , η z T . We thus have the relationship: s r = s r η. 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 l x m r , l y m r and l z m r . Let us define an M × 3 matrix L r whose mth row is equal to l x m r , l y m r , l z m r . This matrix L r , 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 l r = L r η. This column vector l r , 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 array y are expressed as the sum of a magnetic signal b and additive sensor noise, represented by a random vector ε Here, the magnetic signal b is expressed as Here, b S , called the signal vector, represents the biomagnetic signal that is the target of the measurements, and b I , called the interference vector, represents the interference overlapped onto the signal b S . In this paper, the interference b I 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 generate b S 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 r 1 , … , r Q , their orientations by η 1 , … , η Q , and their magnitudes by s 1 , … , s Q . Then, the source distribution is expressed as where δ r indicates the Dirac delta function. The signal vector b S is expressed as where l q represents the lead field vector of the qth source obtained such that l q = L r q η q .

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 at r, B(r), is expressed using the spherical polar coordinate r = r, θ, ϕ by 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 righthand 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 as r min D and the radial coordinates of the sensor farthest from the origin as r max D . The internal region is formally defined as the region with r < r min D and the external region as the region with r > r max D . The region with r min D < r < r max D is called the intermediate region.
Let us derive the SSS basis vectors. To do so, the magnetic signal detected by the jth sensor is denoted by b j and the location and the normal vector of the jth sensor by r j and ζ j . Then, we have 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 b j is obtained as a surface integral of B r j · ζ j over the area of the jth pick-up coil). Here, b j int and b j ext , respectively, represent magnetic components originating from the internal and external regions. These components are expressed as where we set μ 0 = 1 for simplicity. Let us define the internal and external components of the vector b as where column vectors c ℓ,m and d ℓ,m are given by Truncating the summation with respect to the multiple order ℓ to L C for b int and L D for b ext , we finally obtain Thus, defining where S = C, D and x = α T , β T T . Here, C is an M × N C matrix, and D is an M × N D matrix, where Note that the truncation values L C and L D correspond to the highest spatial frequencies possibly contained in b int and b ext [9,17], respectively. Therefore, setting these parameters at too low values may result in an insufficient representation of the signal vectors b int and b ext . The effects of L C and L D for data from the 306-channel Elekta Neuromag have been investigated and values of L C = 8 and L D = 3 were found to be sufficient for such data sets in [8,9]. 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 a i , … , a j T by using the ith to jth components of a = a 1 , … , a M T as a i j (namely, . From (14), we havê With a small positive constant κ, the relationship holds. Then, using the matrix inversion formula we get Using (15) and (20), we obtain Thus, the internal component b int can be extracted by multiplying with the magnetic-field data b. That is, the matrix P int acts as a projector that passes the internal components and blocks the external ones (note that since P int 2 = P int and P int T = P int do not hold, P int is not actually a projector). Therefore, we call P int the SSS signal extractor in this paper.
In exactly the same manner, we can derivê and the SSS external-signal extractor P ext is derived as This P ext passes the external components but blocks the internal ones.

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.) 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 b S ∈ span C , 25 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 b s belongs to the internal subspace, b s is expressed as a linear combination of the column vectors of C where c j is the jth column of C, α j is the jth expansion coefficient, and α is a column vector containing the coefficients (i.e., α = α 1 , … , α N C T ). Thus, denoting an N D × 1 column vector whose elements are all zero by 0, we can derive the relationship The equation above indicates that the SSS signal extractor P int passes the signal vector b S with no distortion.
The assumption that the interference sources are located within the external region leads to where span D is referred to as the external subspace. Since the interference vector b I belongs to the external subspace, the interference vector b I is expressed as where d j is the jth column of D, β j is the jth expansion coefficient, and β is a column vector containing coefficients Again denoting the N C × 1 column vector whose elements are all zero by 0, we have the relationship The equation above indicates that the extractor P int completely blocks the interference vector b I .
Consequently, using (1) and (2), we show The equation above indicates that, by multiplying the extractor P int with the data vector y, the signal vector b S is selectively extracted with no distortion.
In (31), ε ′ ε ′ = P int ε indicates the noise in the SSScleaned data. Assuming that the sensor noise ε is Gaussian with the covariance matrix σ 2 I, the covariance matrix of ε′, Σ ε′ is derived such that where we have εε T = σ 2 I and the notation ⋅ indicates averaging. The equation above shows that P int P T int 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 where Σ ε′ j,j indicates the jth diagonal element of the matrix Σ ε′ . 2.5. SSS Method for Flat Sensor Arrays. In Figure 1(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, l q , may have components expanded by the columns of D, as well as components expanded by the columns of C, resulting in where α ′ and β ′ are column vectors containing the expansion coefficients. Therefore, by applying the SSS extractor P int to l q , we have That is, the extractor P int changes the lead field of this signal source from l q to l q . Assuming that Q signal sources exist, the signal vector b S t is expressed as This signal vector changes to b S t , which is given by It is clear here that applying the SSS extractor P int distorts the signal vector b S t . 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 Figures 1(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 b I t does not have components belonging to span C , and applying P int 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 will never be attained.
We can evaluate the performance of the SSS method for data from flat sensor arrays by checking how close the SSSprocessed results come to (38). Namely, we define the performance measures such that In the equations above, G S is called the signal gain. Ideally, G S is equal to 1, and deviation of G S from 1 is a measure of performance degradation of the SSS method. G I is called the interference gain. Ideally, G I is equal to zero, indicating that the method completely blocks the interference. Thus, a deviation of G I from 0 is a measure of the performance degradation. Note that 1/G I has been often called the shield factor in the previous literature [9,18]. We use G S and G I (or 1/G I , 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.

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 Figure 1(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 Figure 2(a). The sensor array consists of 64 sensors arranged in an 8 × 8 configuration on the plane z = 10 cm; 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][12][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.
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 Figure 2(b). The signal sensor data b S (t) were computed by projecting the source time course through the sensor lead field obtained using the Biot-Savart law. The signal data b S (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, b S t + ε, 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, σ 2 I , 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 Table 1. 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 b I t were computed. The generated interference data are shown in Figure 3(c).

Sensors
Source space  We applied the SSS extractors to b S t and b I t in order to check how large the internal and external components are in each of b S t and b I t . The upper and lower panels, respectively, of Figure 4(a) show P int b S t and P ext b S t . Here, P int b S t and P ext b S t indicate the internal and external components contained in b S t . These results show that the signal vector b S t 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 P int b I t and P ext b I t . Here, P int b I t is nearly equal to zero, and P ext b I t is almost the same as b I t . These results verify our arguments in Section 2.5 that b I t contains almost no internal components. 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 Figure 5(a). The upper and lower panels, respectively, of Figure 5(b) show P int b S t and P ext b S t , and the upper and lower panels, respectively, of Figure 5(c) show P int b I t and P ext b I t . These results clearly confirm that the amount of the external components in b S t , as well as the amount of the internal components in b I t , is very small, explaining why the SSS method works well for data from helmet-type sensor arrays used in MEG.

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 component P int b S t and the external component P ext b S t were computed with the origin set at four different locations of (0, 0, and z ori ) where z ori 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) b S t + ε shown in Figure 2(c) were used. The truncation values L C and L D 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. The amount of signal leakage can be assessed by the signal gain G S defined in (39). Let us derive a quantitative relationship of the signal gain G S versus z ori . To do this, voxels with 0.5 cm intervals were assumed in a 3dimensional 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 b S t were computed by setting the signal source at one of the voxel locations, and the signal gain G S was computed using b S t . The mean signal gain was computed by averaging the signal gains obtained from all voxel locations. Here, the SSS extractor was derived with z ori varied from 0 cm to 10 cm. The mean signal gain versus z ori is plotted with a broken line in Figure 7(a). It can be seen in these plots that G S gradually decreases as the origin becomes closer to the sensor plane.
The noise gain G ε 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 G S /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 z ori approximately equal to 9 cm. We next performed computer simulation to derive the relationship between the interference gain G I and the origin location z ori . To do this, we assumed a sphere with its radius equal to r I 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 b I t were computed by assuming that a single interference source exists on the surface. The interference gain G I 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 z ori . Here, denoting the polar coordinate of the interference source by (r I , θ, ϕ), 100 locations of the interference source were determined as locations with 10 equal-interval θ and 10 equalinterval φ.
The plots of G I with respect to z ori are shown in Figure 7(b). Here, two cases, r I = 5 m and r I = 20 m, are shown. The case r I = 5 m represents cases in which an interference source is located relatively near to the sensors, and the case r I = 20 m 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 z ori ≈ 9 cm for both cases. investigations [8,9], and so, L C and L D were, respectively, set at 7 and 3 in the preceding subsections. In Figure 8, the signal gain G S , noise gain G ε , and their ratio G S /G ε are plotted with respect to z ori for three different values of L C (L C = 5, 6, and 7) and two different values of L D (L D = 2 and 3). 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 L D = 2, and the solid lines indicate those with L D = 3. It can be seen that the noise gain is considerably smaller when L D = 2 than when L D = 3. The signal gain when L D = 2 is greater than when L D = 3. The differences in the signal and noise gains among different values of L C are generally small. These results suggest that L D should be set at 2, rather than 3. In Figure 8(c), the plots of the gain ratio G S /G ε for L D = 2 are shown to have peak maxima when z ori > 8, regardless of the value of L C . The results in Figure 8 suggest that, as far as the signal and noise gains concerned, L D = 2 should be the best choice, and any value of L C = 5, 6, and 7 may be chosen if we set z ori such that z ori > 8.
The interference gain is plotted with respect to z ori for the three values of L C and the two values of L D in Figure 9.
Here, the cases r I = 5 m and r I = 20 m are, respectively, shown in Figures 9(a) and 9(b), where the broken lines indicate the results of L D = 2 and the solid lines indicate those of L D = 3. These plots show that regardless of the values of L C and L D , there is a general tendency that the gain decreases as the origin becomes closer to the sensors, and it reaches a minimum near z ori equal to 9 cm. Namely, the value of z ori ≈ 9 cm gives the minimum interference gain (the maximum shield factor) for all values of L C and L D . On the basis of the results in Figures 8 and 9, we can conclude that the origin parameter z ori 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.

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 G I is investigated by using erroneous sensor locations and  orientations when computing the SSS basis vectors. Here, the sensor location error E loc is assessed using where r d j is the true location of the jth sensor, andr d j is the calibrated location of this sensor. Note thatr d j is assumed to contain an error. Here, r d j is obtained by adding (small) random displacements tor d j . The signal vector b S t and the interference vector b I t were computed using the true location r d j , while the SSS basis vectors were computed using the calibrated locationr d j . In exactly the same manner, the error in the sensor orientation, E ori , is assessed using 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 b S t and the interference vector b I t were computed using ζ j , while the SSS basis vectors were computed using the erroneous orientationζ j . In Figure 10, the interference gain is plotted with respect to r I , the distance to the interference source, for the three values of L C and the two values of L D . Again, the broken lines indicate the results with L D = 2, and the solid lines indicate those with L D = 3. Here, the plots in Figure 10(a) indicate the case of no calibration error (E loc = E ori = 0). The plots in Figures 10(b)-10(d), respectively, indicate the cases E loc = E ori = 0 03%, E loc = E ori = 0 1%, and E loc = E ori = 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.
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/G I for r I > 15 m is greater than 10 4 with L D = 2 but it drops to less than 30 when E loc = E ori = 1%. Regarding the optimal choices of the parameters L C and L D , the choice of L D = 2 mostly gives greater shield factor than L D = 3. There are no significant differences among the choices of L C , but the choice of L C = 6 gives slightly better results when E loc = E ori > 0 1%. Since such errors as E loc = E ori > 0 1% may be considered the practical values of sensor calibration errors, the choices of L C = 6 and L D = 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 L C = 6 and L D = 2 also give the best results.
With the choices of L C = 6, L D = 2, and z ori =9 cm, the interference gain G I is replotted with respect to r I , the distance to the interference source. The results are shown in Figure 11(a). It can be seen that the shield factor (1/G I ) exceeds 10 4 for r I > 15 m 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, L C = 6, L D = 2, and 1% sensor calibration errors, signal and noise gains versus z ori 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.

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 largerscale 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 the x, 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 gain G I versus r I (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 (z ori ) were plotted in Figure 12 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).
The SSS performance for the vector sensor array was investigated. The interference gain G I versus r I was plotted in Figure 13(a). We can observe significant improvements in the shielding capability. The current-density map obtained from the SSS interference-removal results P int y t with the original lead field L r . (f) The current-density map obtained using P int y t with the SSS-modified lead field L r = P int L r . Here, two-dimensional current-density maps on the plane z = 2 cm, which is the plane 8 cm below the sensor plane, are reconstructed using the field map at t = 1200. The signal data b S t were computed assuming two sources, located at (−5 cm, 4 cm, and 3 cm) and (4 cm, −3 cm, and 3 cm).
That is, with 1% calibration errors, the vector sensor array has a shield factor of approximately 500 (G I = 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  Figure 2(c). (c) Source reconstruction results obtained using the SSS interference removal results P int y t with the original lead field L r . (d) Source reconstruction results obtained using the SSS interference removal results P int y t 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.
SSS performance improvements due to the use of tangential sensors. The signal and noise gains with respect to the origin's z coordinate (z ori ) 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 Section 2.5, the SSS-processed signal vector b S t b S t = P int b S t becomes distorted because the SSS application causes the removal of the external components from the signal vector b S t . 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.
The signal data b S t 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 b S t + ε are shown in the upper panel of Figure 14(a). The interference b I (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 y t ; y t = b S t + ε + ξb I t , where a positive constant ξ controls the signal-to-interference ratio (SIR) (The SIR is defined as the ratio b S t / ξb I t in this computer simulation), which was set equal to 0.25 in this computer simulation. The SSS interference suppression results P int y t are shown in the lower panel of Figure 14(b).
The current-density reconstruction was performed using RENS beamformer, proposed in [20,21]. The currentdensity map obtained from the signal plus sensor-noise data, b S t + ε, is shown in Figure 14(c). This map works as the ground truth for the following experiments. The currentdensity map obtained from the interference-overlapped data y t 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 results P int y t with the original lead field L r and with the SSS-modified lead field L r : L r = P int L r . The current-density map obtained with the original lead field L r 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 SSSmodified lead field in the voxel space analysis.

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 Figure 2(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 y t 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 y t , and the results are shown in Figure 15(b). A large influence from the interference can be seen here. The source reconstruction was carried out using the SSS interference suppression results P int y t with the original lead field L r . 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 SSSmodified lead field in the voxel space analysis.

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 twodimensional current-density mapping and threedimensional 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 parameters L C and L D , 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 L D = 6 and L D = 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 10 4 for interference originating from fairly distant sources (r I > 15 m) 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 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, 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 maintenancefree magnetocardiography systems, which in the near future may replace the electrocardiogram now routinely used in hospitals.