Journal of Healthcare Engineering
Volume 2018 (2018), Article ID 7689589, 19 pages
https://doi.org/10.1155/2018/7689589

Research Article

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

Kensuke Sekihara1,2

¹Signal Analysis Inc., Hachioji, Tokyo, Japan
²Department of Advanced Technology in Medicine, Tokyo Medical and Dental University, 1-5-45 Yushima, Bunkyo-ku, Tokyo 113-8519, Japan

Correspondence should be addressed to Kensuke Sekihara; k-sekihara@nifty.com

Received 15 December 2017; Accepted 22 February 2018; Published 26 April 2018

Academic Editor: Vincenzo Positano

Copyright © 2018 Kensuke Sekihara. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 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 10⁴ 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.

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

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 the mth sensor as ym. The measurement from the whole sensor array is expressed as a column vector . Here, is the number of sensors, and the superscript 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 . The source magnitude at is denoted by a scalar . The source vector is denoted by s(r), and the source orientation is denoted by . We thus have the relationship: . Let us assume that a unit-magnitude source exists at . When this unit-magnitude source is directed in the x, y, and z directions, the outputs of the mth sensor are, respectively, denoted by and . Let us define an matrix whose mth row is equal to . This matrix , referred to as the lead field matrix, represents the sensitivity of the sensor array at r. When the unit-magnitude source at is oriented in the direction, the outputs of the sensor array are expressed as . This column vector , referred to as the lead field vector, represents the sensitivity of the sensor array in the direction of at the location .

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 is expressed as

Here, , called the signal vector, represents the biomagnetic signal that is the target of the measurements, and , called the interference vector, represents the interference overlapped onto the signal . In this paper, the interference 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 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 discrete sources exist in the source space. Their locations are denoted by , their orientations by , and their magnitudes by . Then, the source distribution is expressed as

where indicates the Dirac delta function. The signal vector is expressed as

where represents the lead field vector of the qth source obtained such that .

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 at r, B(r), is expressed using the spherical polar coordinate by

where indicates the magnetic permeability of free space. In (5), and 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 as and the radial coordinates of the sensor farthest from the origin as . The internal region is formally defined as the region with and the external region as the region with . The region with 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 and the location and the normal vector of the jth sensor by and . 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 is obtained as a surface integral of over the area of the jth pick-up coil). Here, and , respectively, represent magnetic components originating from the internal and external regions. These components are expressed as

where we set for simplicity. Let us define the internal and external components of the vector as and , which are expressed such that

where column vectors and are given by

Truncating the summation with respect to the multiple order to LC for bint and LD for bext, we finally obtain

Thus, defining

we obtain

where and . Here, is an matrix, and is an matrix, where

Note that the truncation values and correspond to the highest spatial frequencies possibly contained in and [9, 17], respectively. Therefore, setting these parameters at too low values may result in an insufficient representation of the signal vectors and . The effects of and for data from the 306-channel Elekta Neuromag have been investigated and values of and 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 and from given magnetic signal data b. That is, the least squares estimate is obtained as

Then, and are estimated as

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 by using the ith to jth components of as (namely, ). From (14), we have

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 can be extracted by multiplying

with the magnetic-field data . That is, the matrix acts as a projector that passes the internal components and blocks the external ones (note that since and do not hold, is not actually a projector). Therefore, we call the SSS signal extractor in this paper.

In exactly the same manner, we can derive

and the SSS external-signal extractor is derived as

This 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.)

(a)

(b)

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.

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

holds, where the notation span indicates the span of the column vectors of . This span is referred to as the internal subspace in this paper. That is, since the signal vector belongs to the internal subspace, is expressed as a linear combination of the column vectors of

where is the jth column of , is the jth expansion coefficient, and is a column vector containing the coefficients (i.e., ). Thus, denoting an column vector whose elements are all zero by 0, we can derive the relationship

The equation above indicates that the SSS signal extractor passes the signal vector with no distortion.

The assumption that the interference sources are located within the external region leads to

where span is referred to as the external subspace. Since the interference vector belongs to the external subspace, the interference vector is expressed as

where is the jth column of , is the jth expansion coefficient, and is a column vector containing coefficients . Again denoting the column vector whose elements are all zero by 0, we have the relationship

The equation above indicates that the extractor completely blocks the interference vector .

Consequently, using (1) and (2), we show

The equation above indicates that, by multiplying the extractor with the data vector y, the signal vector is selectively extracted with no distortion.

In (31), indicates the noise in the SSS-cleaned data. Assuming that the sensor noise is Gaussian with the covariance matrix , the covariance matrix of is derived such that

where we have and the notation indicates averaging. The equation above shows that 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 such that

where 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, , may have components expanded by the columns of , as well as components expanded by the columns of , resulting in

where and are column vectors containing the expansion coefficients.

Therefore, by applying the SSS extractor to , we have

That is, the extractor changes the lead field of this signal source from to . Assuming that signal sources exist, the signal vector is expressed as

This signal vector changes to , which is given by

It is clear here that applying the SSS extractor distorts the signal vector . 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 does not have components belonging to span , and applying 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 SSS-processed results come to (38). Namely, we define the performance measures such that

In the equations above, is called the signal gain. Ideally, is equal to 1, and deviation of from 1 is a measure of performance degradation of the SSS method. is called the interference gain. Ideally, is equal to zero, indicating that the method completely blocks the interference. Thus, a deviation of 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 and (or 1/GI, the shield factor), as well as the noise gain 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 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 ; 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.

(a)

(b)

(c)

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 plane , 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 (, , ) 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 are shown in the upper panel, and the signal plus sensor-noise data 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 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, 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 , where 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 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)
214

(a)

(b)

(c)

Figure 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 data . 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.

We applied the SSS extractors to and in order to check how large the internal and external components are in each of and . The upper and lower panels, respectively, of Figure 4(a) show and . Here, and indicate the internal and external components contained in . These results show that the signal vector 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 and . Here, is nearly equal to zero, and is almost the same as . These results verify our arguments in Section 2.5 that contains almost no internal components.

(a)

(b)

Figure 4: Results of experiments that apply the SSS extractors to and computed assuming the flat sensor array in Figure 2(a). (a) The upper and lower panels, respectively, show and , which indicate the internal and external components in . (b) The upper and lower panels, respectively, show and , which indicate the internal and external components in . 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.

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 and , and the upper and lower panels, respectively, of Figure 5(c) show and . These results clearly confirm that the amount of the external components in , as well as the amount of the internal components in , is very small, explaining why the SSS method works well for data from helmet-type sensor arrays used in MEG.

(a)

(b)

(c)

Figure 5: Results of the same experiments in which the SSS extractors are applied to and , 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 and . (c) The upper and lower panels, respectively, show and . 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.

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 component and the external component 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) shown in Figure 2(c) were used. The truncation values and 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 component and the external component 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 . 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 gain defined in (39). Let us derive a quantitative relationship of the signal gain 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 were computed by setting the signal source at one of the voxel locations, and the signal gain was computed using . 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.

(a)

(b)

Figure 7: (a) Plots of the mean signal gain , the noise gain , and their ratio versus the z coordinate of the origin zori. The signal sensor data 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 versus zori. The plot with the broken line shows the case , and the plot with the solid line shows the case . 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 , and the interference gain was obtained. The mean interference gain was computed by averaging obtained from all 100 source locations. The sphere with a radius is shown in the upper part.

The noise gain 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 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 gain GI and the origin location zori. To do this, we assumed a sphere with its radius equal to 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 were computed by assuming that a single interference source exists on the surface. The interference gain 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 (), 100 locations of the interference source were determined as locations with 10 equal-interval and 10 equal-interval φ.

The plots of GI with respect to zori are shown in Figure 7(b). Here, two cases, and , are shown. The case represents cases in which an interference source is located relatively near to the sensors, and the case 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 of and , Truncated Values of Multipole-Order

We explore the optimal values of the parameters and , 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, and were, respectively, set at 7 and 3 in the preceding subsections. In Figure 8, the signal gain GS, noise gain Gε, and their ratio are plotted with respect to zori for three different values of () and two different values of (). 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 , and the solid lines indicate those with .

(a)

(b)

(c)

Figure 8: (a) The signal gain , (b) noise gain , and (c) their ratio plotted with respect to the origin’s z coordinate, zori, for three values of () and two values of (). The broken lines indicate the results with , and the solid lines indicate those with .

It can be seen that the noise gain is considerably smaller when than when . The signal gain when is greater than when . The differences in the signal and noise gains among different values of LC are generally small. These results suggest that should be set at 2, rather than 3. In Figure 8(c), the plots of the gain ratio for are shown to have peak maxima when zori > 8, regardless of the value of . The results in Figure 8 suggest that, as far as the signal and noise gains concerned, should be the best choice, and any value of may be chosen if we set zori such that zori > 8.

The interference gain is plotted with respect to zori for the three values of and the two values of in Figure 9. Here, the cases and are, respectively, shown in Figures 9(a) and 9(b), where the broken lines indicate the results of and the solid lines indicate those of . These plots show that regardless of the values of and , 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 and . 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.

(a)

(b)

Figure 9: The interference gain plotted with respect to the origin’s z coordinate, zori, for three values of () and two values of (). (a) The distance to the interference source is set at 5 m. (b) The distance to the interference source is set at 20 m. The broken lines indicate the results with , and the solid lines indicate those with .

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 is investigated by using erroneous sensor locations and orientations when computing the SSS basis vectors. Here, the sensor location error is assessed using

where is the true location of the jth sensor, and is the calibrated location of this sensor. Note that is assumed to contain an error. Here, is obtained by adding (small) random displacements to . The signal vector and the interference vector were computed using the true location , while the SSS basis vectors were computed using the calibrated location . In exactly the same manner, the error in the sensor orientation, , is assessed using

where ζ j is the true normal vector of the jth sensor, and ζbj is the calibrated normal vector, which is . Here, is obtained by adding random vectors to . The signal vector and the interference vector were computed using , while the SSS basis vectors were computed using the erroneous orientation .

In Figure 10, the interference gain is plotted with respect to , the distance to the interference source, for the three values of and the two values of . Again, the broken lines indicate the results with , and the solid lines indicate those with . Here, the plots in Figure 10(a) indicate the case of no calibration error (). The plots in Figures 10(b)–10(d), respectively, indicate the cases , , and . 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.

(a)

(b)

(c)

(d)

Figure 10: The plots of the interference gain GI with respect to the distance to the interference source, , for the three values of and the two values of . (a) No calibration errors (). (b) Calibration errors of . (c) Calibration errors of . (d) Calibration errors of . The broken lines indicate the results with and the solid lines indicate those with . The origin was set at (0 cm, 0 cm, and 9 cm) ().

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 10⁴ with but it drops to less than 30 when . Regarding the optimal choices of the parameters and , the choice of mostly gives greater shield factor than . There are no significant differences among the choices of , but the choice of gives slightly better results when . Since such errors as may be considered the practical values of sensor calibration errors, the choices of and 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 and also give the best results.

With the choices of , , and zori=9 cm, the interference gain is replotted with respect to , 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 10⁴ for 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, , , 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.

(a)

(b)

Figure 11: (a) The interference gain replotted with respect to , the distance to interference sources. The choices of , and 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 and .

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 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 versus (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).

(a)

(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 gain plotted with respect to for six different calibration errors. The choices of , 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 and .

The SSS performance for the vector sensor array was investigated. The interference gain GI versus rI was plotted in Figure 13(a). We can observe significant improvements in the shielding capability.

(a)

(b)

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 gain GI plotted with respect to rI for six different calibration errors. The choices of , , and 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 and .

That is, with 1% calibration errors, the vector sensor array has a shield factor of approximately 500 (). 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 Section 2.5, the SSS-processed signal vector becomes distorted because the SSS application causes the removal of the external components from the signal vector . 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.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 14: (a) Upper panel: the signal plus sensor-noise time courses . Lower panel: the interference time courses . (b) Upper panel: the interference-overlapped sensor time courses . Lower panel: the SSS interference removal results . These time courses are normalized to the maximum value in each panel. (c) The current-density map obtained from . (d) The current-density map obtained from interference-overlapped sensor time courses . (e) The current-density map obtained from the SSS interference-removal results with the original lead field . (f) The current-density map obtained using with the SSS-modified lead field . Here, two-dimensional current-density maps on the plane , which is the plane 8 cm below the sensor plane, are reconstructed using the field map at . The signal data were computed assuming two sources, located at (−5 cm, 4 cm, and 3 cm) and (4 cm, −3 cm, and 3 cm).

The signal data 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 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 ; , where a positive constant ξ controls the signal-to-interference ratio (SIR) (The SIR is defined as the ratio in this computer simulation), which was set equal to 0.25 in this computer simulation. The SSS interference suppression results 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, , 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 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 with the original lead field and with the SSS-modified lead field : . The current-density map obtained with the original lead field 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 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 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 , and the results are shown in Figure 15(b). A large influence from the interference can be seen here.

(a)

(b)

(c)

(d)

Figure 15: (a) Source reconstruction results obtained using the signal plus sensor-noise time courses in Figure 2(c). (b) Source reconstruction results obtained using the sensor time courses 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 with the original lead field . (d) Source reconstruction results obtained using the SSS interference removal results with the SSS-modified lead field . The RENS beamformer [20, 21] was applied to the field data at for three-dimensional source reconstruction.

The source reconstruction was carried out using the SSS interference suppression results with the original lead field . 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 , 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 parameters and , 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 and 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⁴ for interference originating from fairly distant sources () 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 maintenance-free magnetocardiography systems, which in the near future may replace the electrocardiogram now routinely used in hospitals.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

References

P. P. Freitas, R. Ferreira, S. Cardoso, and F. Cardoso, “Magnetoresistive sensors,” Journal of Physics: Condensed Matter, vol. 19, no. 16, article 165221, 2007.

P. Campiglio, L. Caruso, E. Paul et al., “GMR-based sensors arrays for biomagnetic source imaging applications,” IEEE Transactions on Magnetics, vol. 48, no. 11, pp. 3501–3504, 2012.

M. Pannetier-Lecoeur, C. Fermon, P. Campiglio, Q. Herreros, and G. Jasmin-Lebras, “Spin electronics based magnetic sensors for biomagnetic measurements,” in Magnetoencephalography, pp. 1001–1005, Springer, Berlin, Heidelberg, 2014.

S. Kawabata, S. Ushio, T. Shibuya, Y. Adachi, K. Sekihara, and A. Okawa, “Measurement of magnetomyography using an array of magnetoresistive(MR) sensor,” in Abstract Book for the 20th International Conference on Biomagnetism, Biomag 2016, p. 80, Seoul, Korea, 2016.

Y. Shirai, K. Hirano, T. Shibuya, S. Ushio, and S. Kawabata, “Measurement of magnetocardiogram using magnetoresistive sensor,” in Abstract Book for the 20th International Conference on Biomagnetism, Biomag 2016, p. 80, Seoul, Korea, 2016.

M. Hämäläinen, R. Hari, R. J. Ilmoniemi, J. Knuutila, and O. V. Lounasmaa, “Magnetoencephalography-theory, instrumentation, and applications to noninvasive studies of the working human brain,” Reviews of Modern Physics, vol. 65, no. 2, pp. 413–497, 1993.

K. Sekihara and S. S. Nagarajan, “Subspace-based interference removal methods for a multichannel biomagnetic sensor array,” Journal of Neural Engineering, vol. 14, no. 5, article 051001, 2017.

S. Taulu and M. Kajola, “Presentation of electromagnetic multichannel data: the signal space separation method,” Journal of Applied Physics, vol. 97, no. 12, article 124905, 2005.

S. Taulu, J. Simola, and M. Kajola, “Applications of the signal space separation method,” IEEE Transactions on Signal Processing, vol. 53, no. 9, pp. 3359–3372, 2005.

S. Taulu, M. Kajola, and J. Simola, “Suppression of interference and artifacts by the signal space separation method,” Brain Topography, vol. 16, no. 4, pp. 269–275, 2004.

A. Kandori, K. Ogata, Y. Watanabe et al., “Space-time database for standardization of adult magnetocardiogram-making standard MCG parameters,” Pacing and Clinical Electrophysiology, vol. 31, no. 4, pp. 422–431, 2008.

K. Ogata, A. Kandori, T. Miyashita, K. Sekihara, and K. Tsukada, “A comparison of two-dimensional techniques for converting magnetocardiogram maps into effective current source distributions,” Review of Scientific Instruments, vol. 82, no. 1, article 014302, 2011.

Y. Ito, K. Shiga, K. Yoshida et al., “Development of a magnetocardiography-based algorithm for discrimination between ventricular arrhythmias originating from the right ventricular outflow tract and those originating from the aortic sinus cusp: a pilot study,” Heart Rhythm, vol. 11, no. 9, pp. 1605–1612, 2014.

W.-D. Bang, K. Kim, Y.-H. Lee et al., “Repolarization heterogeneity of magnetocardiography predicts long-term prognosis in patients with acute myocardial infarction,” Yonsei Medical Journal, vol. 57, no. 6, pp. 1339–1346, 2016.

D. Kim and H. Ah, “Current status and future of cardiac mapping in atrial fibrillation,” in Atrial Fibrillation-Basic Research and Clinical Applications, InTech, 2012.

E. L. Hill, “The theory of vector spherical harmonics,” American Journal of Physics, vol. 22, no. 4, pp. 211–214, 1954.

T. E. Özkurt, M. Sun, and R. J. Sclabassi, “Decomposition of magnetoencephalographic data into¨ components corresponding to deep and superficial sources,” IEEE Transactions on Biomedical Engineering, vol. 55, no. 6, pp. 1716–1727, 2008.

J. Nurminen, S. Taulu, and Y. Okada, “Effects of sensor calibration, balancing and parametrization on the signal space separation method,” Physics in Medicine & Biology, vol. 53, no. 7, pp. 1975–1987, 2008.

J. Nurminen, S. Taulu, J. Nenonen, L. Helle, J. Simola, and A. Ahonen, “Improving MEG performance with additional tangential sensors,” IEEE Transactions on Biomedical Engineering, vol. 60, no. 9, pp. 2559–2566, 2013.

I. Kumihashi and K. Sekihara, “Array-gain constraint minimum-norm spatial filter with recursively updated gram matrix for biomagnetic source imaging,” IEEE Transactions on Biomedical Engineering, vol. 57, no. 6, pp. 1358–1365, 2010.

K. Sekihara and S. S. Nagarajan, Electromagnetic Brain Imaging: A Bayesian Perspective, Springer-Verlag, Berlin, Heidelberg, 2015.