*Article* **Processing Chain for Localization of Magnetoelectric Sensors in Real Time**

**Christin Bald and Gerhard Schmidt \***

Digital Signal Processing and System Theory, Institute of Electrical Engineering and Information Technology, Faculty of Engineering, Kiel University, 24143 Kiel, Germany; cbal@tf.uni-kiel.de **\*** Correspondence: gus@tf.uni-kiel.de; Tel.: +49-431-880-6125

**Abstract:** The knowledge of the exact position and orientation of a sensor with respect to a source (distribution) is essential for the correct solution of inverse problems. Especially when measuring with magnetic field sensors, the positions and orientations of the sensors are not always fixed during measurements. In this study, we present a processing chain for the localization of magnetic field sensors in real time. This includes preprocessing steps, such as equalizing and matched filtering, an iterative localization approach, and postprocessing steps for smoothing the localization outcomes over time. We show the efficiency of this localization pipeline using an exchange bias magnetoelectric sensor. For the proof of principle, the potential of the proposed algorithm performing the localization in the two-dimensional space is investigated. Nevertheless, the algorithm can be easily extended to the three-dimensional space. Using the proposed pipeline, we achieve average localization errors between 1.12 cm and 6.90 cm in a localization area of size 50 cm × 50 cm.

**Keywords:** localization; magnetoelectric sensors; real time; pose estimation

**Citation:** Bald, C.; Schmidt, G. Processing Chain for Localization of Magnetoelectric Sensors in Real Time. *Sensors* **2021**, *21*, 5675. https:// doi.org/10.3390/s21165675

Academic Editor: Daniel Ramos

Received: 30 July 2021 Accepted: 20 August 2021 Published: 23 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

For the correct solution of inverse problems, such as source reconstruction of biomedical sources, it is essential to know the exact position and orientation of the measuring sensors with respect to the source besides measuring the biomedical signals. Especially in magnetic measurements the sensors do not necessarily have a fixed position and orientation. Thus, a determination of the position and orientation at the beginning of a measurement is not sufficient. Much more desirable is a continuous estimation of the sensor's position and orientation simultaneously with the measurement [1,2].

Magnetic tracking systems are used in many applications, e.g., in indoor positioning systems [3,4] or to locate medical devices inside the body [5,6]. Moreover, magnetic localization approaches are used to determine the position of the subject relative to the sensor array in biomagnetic measurements. A procedure for determining the subject relative to the measuring sensor array, either once at the beginning or also simultaneously with the measurement, was presented in [7]. In [1,2] a method for determining the positions of the individual sensors in a flexible on-scalp MEG system relative to the subject was investigated, which can also be applied during measurement.

Until now, mainly SQUIDs (Super Conducting Quantum Interference Devices) [8] and recently OPMs (Optically Pumped Magnetometers) [9,10] are used for the measurement of biomagnetic signals. Unfortunately, these sensors require a magnetically well shielded environment and are therefore inconvenient in operation. Magnetoelectric sensors, on the other hand, do not require any shielding, no expensive cooling system, and are also very small in size, which makes them ideal for array applications. The sensors are composed of a magnetostrictive and a piezoelectric layer and use the resonant structure of a cantilever [11]. Detection limits in the sub-nT regime have been reached recently [12–15] using modulation techniques such as the ΔE-effect [16] for the detection of low-frequency magnetic fields. Consequently, magnetoelectric sensors could be a promising alternative for biomagnetic field measurements in the near future [17].

Since the feasibility of measurements with simultaneous localization using magnetoelectric sensors has been shown in a previous work [15], this work will focus on the processing chain for an enhanced localization of magnetoelectric sensors in real time. However, by adapting the coil excitation signals the method can be used for arbitrary magnetic field sensors. This contribution will step through the general processing overview shown in Figure 1.

**Figure 1.** General overview of a medical system operating with magnetic sensors. The measurements are performed simultaneously with localizing the sensors. After transforming the signals into the digital domain, the signals are processed and analyzed. Since the analysis of the measured magnetic signals is not in the focus of this contribution, the corresponding box is depicted in gray.

In Section 2, the magnetoelectric sensor used in this paper will be presented and characterized. After explaining the so-called forward problem in Section 3, the real-time localization approach will be presented in Section 4. The presented localization processing chain will be verified by measurements presented in Section 5. The paper closes with a conclusion and an outlook in Section 6.

#### **2. Magnetoelectric Sensor**

For the proof of principle an exchange bias ΔE-effect sensor as depicted in Figure 2a will be used in this contribution. A sensor of the same type has already been used in a previous work [15]. The sensor is based on a polysilicon cantilever with a size of 1 mm width, 3 mm length, and 50 μm thickness. The cantilever is covered by a 4 μm thick magnetic multilayer (20 × Ta/Cu/MnIr/FeCoSiB) and a 2 μm thick piezoelectric layer (AlN). Further details about the fabrication process of the sensor can be found in [15]. The magnetic multilayer consists of ferromagnetic and antiferromagnetic layers, which ensure the self-biasing of the sensor and thus lead to a shift of the magnetization curve of the sensor [18]. Hence, the sensor can be operated without applying an external bias field, which is especially favorable regarding array applications. The sensor is connected to a low-noise JFET charge amplifier [19] and placed on a printed circuit board. The whole sensor system is encapsulated in a 2.1 mm thick brass cylinder for electrical shielding.

As shown in [15], the localization of the magnetoelectric sensor can be performed simultaneously with a measurement without loss of information or degradation of the signals. The first bending mode was used to localize the sensor, while an artificial heart signal was measured in the second mode using the ΔE-effect. Hence, also in this contribution only frequencies around the first bending mode will be used for the transmission of the localization signals. The amplitude and phase response around the first bending mode of the magnetoelectric sensor used are shown in Figure 2c. The characterization measurements have been performed in a magnetically, electrically, and acoustically shielded environment [11]. The magnitude and phase response of the sensor have been measured applying a magnetic field of *b*ac = 1 μT on the sensor's long axis. The performance of the sensor can be determined as described in [20]. The sensor has a resonance frequency of *<sup>f</sup>*<sup>r</sup> <sup>=</sup> 7.712 kHz and a <sup>−</sup>3 dB bandwidth of *bw*–3dB <sup>=</sup> 10.2 Hz. Since the brass cylinder

acts as a low-pass filter with a cut-off frequency of approximately *f*<sup>c</sup> ≈ 1.5 kHz [15], the sensor's performance will be improved when removing the brass cylinder. Nevertheless, the encapsulation is necessary for electrical shielding and acts as a mechanical protection. Moreover, the limit-of-detection in the first bending mode with brass cylinder is still sufficient, because the coils can simply emit higher field amplitudes. The maximum sensitivity of the sensor is reached when a magnetic field is applied on the sensitive axis of the sensor. However, the sensitive axis of the sensor is not necessarily equal to the long axis of the sensor. There can be a tilt *γ* between these two axes [15,21], which is visualized in Figure 2b.

**Figure 2.** (**a**) Exchange bias ΔE-effect sensor used in this study. The sensor is based on a cantilever of size 3 mm × 1 mm. The cantilever is placed on a printed circuit board and connected to a low-noise JFET charge amplifier [19]. The sensor is encapsulated by a brass cylinder for electrical shielding and mechanical protection. (**b**) Visualization of the relationship between the sensitive and the long axis of the sensor. (**c**) Magnitude and phase response of the sensor in the first bending mode applying a magnetic field of *b*ac = 1 μT. The sensor has a resonance frequency of *f*<sup>r</sup> = 7.712 kHz and a bandwidth of *bw*–3dB = 10.2 Hz.

#### **3. Forward Problem**

For the localization of the magnetoelectric sensor, coils are placed outside the localization area as shown in Figure 3. If the distance between the sensor and the coil is large enough, the magnetic field of the coil *i* at the sensor position*rs*(*t*) can be approximated by the field of a magnetic dipole [22]:

$$\vec{b}\_{i}(t, \vec{r}\_{s}(t)) = \frac{\mu\_{0}}{4\pi} \frac{\mathfrak{d}\vec{r}\_{\rm cs,i}(t) \left(\vec{m}\_{\rm c,i}(t)^{\rm T}\vec{r}\_{\rm cs,i}(t)\right) - \vec{m}\_{\rm c,i}(t) \left\|\vec{r}\_{\rm cs,i}(t)\right\|\_{2}^{2}}{\left\|\vec{r}\_{\rm cs,i}(t)\right\|\_{2}^{5}} \tag{1}$$

Here, *μ*<sup>0</sup> is the permeability of vacuum, *m* c,*i*(*t*) the magnetic dipole moment of the coil *i*, and*<sup>r</sup>*cs,*i*(*t*) = *<sup>r</sup>*s(*t*) <sup>−</sup>*r*c,*<sup>i</sup>* the distance vector between the sensor at position*r*s(*t*) and the coil *i* at position*r*c,*i*. The superscript T denotes the transpose of the vector. The positions of the coils*r*c,*<sup>i</sup>* are fixed during the measurement and therefore time independent. In this study, *N*<sup>c</sup> = 6 coils have been used. The coils have an effective diameter of 2.6 cm and consist of about 350 turns of enameled copper wire with a wire cross section of 0.13 mm2. The impedances of the six coils, separated into magnitude and phase, are shown in Figure 4.

The signal measured by the sensor

$$u\_{\rm in}(t) = h\_{\rm s}(t) \* \left(\vec{d}\_{\rm s}^{\rm T}(t) \sum\_{i=1}^{N\_{\rm c}} \vec{b}\_{i}(t, \vec{r}\_{\rm s}(t))\right) \tag{2}$$

can be described as a voltage at the output of the sensor system. At the location of the sensor a superposition of the magnetic fields of the coils is present. Due to the directional characteristic of the sensor *d*<sup>s</sup>(*t*), only a part of the applied magnetic field is picked up. The conversion of magnetic field into voltage by the sensor system (including the charge amplifier) is described by the impulse response *h*s(*t*). Equation (2) is valid at least for the frequencies around the first bending mode [15]. For simplification, no noise sources are considered here.

**Figure 3.** Real (**a**) and schematic (**b**) measurement setup for the localization of magnetoelectric sensors. The coils are placed outside of the localization area and transmit orthogonal signals, which are measured by the sensor. The localization area (box bounded by white stripes in (**a**)) is of size 50 cm × 50 cm.

**Figure 4.** Coil impedances separated into magnitude and phase. In (**a**) the whole spectrum from 100 Hz up to 1 MHz is shown, so that the resonance of the coils can be seen. In (**b**) the frequency range is scaled to the frequency range of the excitation signals.

#### **4. Localization Processing Chain**

For the estimation of the sensor's position and orientation an inverse problem must be solved. The processing chain for solving this inverse problem is shown in Figure 5.

**Figure 5.** Processing chain for the localization of magnetoelectric sensors in real time. The input signal of the sensor is matched filtered with the equalized coil signals. The matched filter outputs at time lag zero are compared with the lead-field matrix entries. A first order Kalman filter smooths the estimated position-orientation-pairs over time to mitigate possible outliers.

For this purpose, the localization area is first divided into a discrete grid containing *N*p different position-orientation-pairs *<sup>P</sup>* = [*p*1,..., *<sup>p</sup><sup>j</sup>* ,..., *<sup>p</sup>N*<sup>p</sup> ], with *<sup>p</sup><sup>j</sup>* = [*r* T p,*j* , *d*- T p,*j* ] <sup>T</sup> consisting of a position vector*r*p,*<sup>j</sup>* (containing x, y, and z components) and an orientation vector *d*p,*<sup>j</sup>* (directivity described by roll *ϕ*, pitch *θ*, and yaw *ψ*). The lead-field matrix

$$A = \begin{bmatrix} a\_{1\times} \dots \omega\_{j\times} \dots \omega\_{N\_{\text{p}}} \end{bmatrix} = \begin{bmatrix} a\_{11} & \dots & a\_{1j} & \dots & a\_{1N\_{\text{p}}} \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ a\_{i1} & \dots & a\_{ij} & \dots & a\_{iN\_{\text{p}}} \\ \vdots & \ddots & \vdots & \ddots & \vdots \\ a\_{N\_{\text{c}}1} & \dots & a\_{N\_{\text{c}}j} & \dots & a\_{N\_{\text{c}}N\_{\text{p}}} \end{bmatrix} \tag{3}$$

describes the forward problem for the defined position-orientation-pairs in *P*. That means, the lead-field matrix entry of row *i* and column *j* is defined as

$$a\_{i\dot{j}} = \vec{d}\_{\text{p},\dot{j}}^{\text{T}} \frac{\text{3 } \vec{r}\_{\text{cp},\dot{j}} \left(\vec{d}\_{\text{c},i}^{\text{T}} \vec{r}\_{\text{cp},\dot{j}}\right) - \vec{d}\_{\text{c},i} ||\vec{r}\_{\text{cp},\dot{j}}||\_2^2}{||\vec{r}\_{\text{cp},\dot{j}}||\_2^5} \tag{4}$$

and thus describes the influence of coil *i* on the sensor, if the sensor would have the position and orientation described by *<sup>p</sup><sup>j</sup>* . The distance vector is defined as*<sup>r</sup>*cp,*ij* =*r*p,*<sup>j</sup>* −*r*c,*<sup>i</sup>* and the orientation vector *d*p,*<sup>j</sup>* can be described by the angles *θ* and *ψ* (*ϕ* is always zero here) using [23]:

$$\vec{d}\_{\mathbb{P}^j} = \begin{pmatrix} \cos(\theta)\cos(\psi) \\ \cos(\theta)\sin(\psi) \\ -\sin(\theta) \end{pmatrix} \tag{5}$$

Equation (4) is a reduced magnetic dipole equation. The prefactor of Equation (1) is neglected and the magnetic dipole moment of the coil is reduced to the orientation of the coil.

#### *4.1. Signal Generation and Equalizer*

To separate the mixed signals received by the sensor, the signals of the coils must be orthogonal. Two signals *xi*(*n*) and *xj*(*n*) are orthogonal for *<sup>n</sup>* ∈ {*L*1, ... , *<sup>L</sup>*2}, if the following condition

$$\sum\_{n=-L\_1}^{L\_2} \mathbf{x}\_i(n) \mathbf{x}\_j^\*(n) \;= \begin{cases} E\_{i\prime} & i = j \\ 0, & i \neq j \end{cases} \tag{6}$$

is fulfilled [24]. Extending this equation to the constant *Ei* being 1, the signals are called orthonormal [24]. This is necessary to extract the individual coil amplitudes from the sensor signal and make them comparable. Different approaches can be used for the generation of orthogonal signals, e.g., using a TDMA (Time Division Multiple Access), an FDMA (Frequency Division Multiple Access) or a CDMA (Code Division Multiple Access) approach [25]. Due to the small bandwidth of the sensor, a TDMA approach is used in this contribution. Thus, the excitation signals

$$\mathbf{x}\_{\rm ex}(n) = \cos(2\pi f\_{\rm r}(n - \kappa\_i)) \,\mathrm{w}(n - \kappa\_i) \quad \text{with } \kappa\_i = (i - 1)L\_{\rm mf} \quad \forall i \in \{1, \ldots, N\_{\rm c}\} \tag{7}$$

are cosine signals at the resonance of the magnetoelectric sensor [15]. The signals are weighted with a Hann window *<sup>w</sup>*(*n*) [26] of length *<sup>L</sup>*sig, so that a smoothed in- and outfading of the signals is ensured. Additionally, the condition *L*mf ≥ *L*sig must be fulfilled. If *L*mf > *L*sig, there is a pausing time between two consecutive coil signals. This is important when considering the impulse response of the sensor. The excitation signals are repeated every *<sup>L</sup>*<sup>r</sup> = *<sup>N</sup>*c*L*mf samples

$$\mathbf{x}\_{\text{out},i}(n) = \mathbf{x}\_{\text{ex},i}(n - \lambda L\_{\mathbf{f}}) \quad \text{with } \lambda \in \mathbb{Z} \tag{8}$$

and transmitted by the coils after D/A conversion. It should be noted that *<sup>L</sup>*<sup>r</sup> = *<sup>L</sup>*mf if an FDMA or a CDMA approach is chosen, because the signals are transmitted simultaneously by all coils.

As can be seen from Equation (1) the magnetic field of a coil is proportional to the driving current. Since the output of the D/A converter is proportional to a voltage, the excitation signals *<sup>x</sup>*ex(*n*)=[*x*ex,1(*n*),..., *<sup>x</sup>*ex,*N*<sup>c</sup> (*n*)]<sup>T</sup> are linearly deformed in amplitude and phase. This can be described by the impulse response of the coil *<sup>h</sup>*c,*i*(*n*), denoting the relationship between the voltage and the current of the coil *i*. Additionally, the signals are modified by the impulse response of the magnetoelectric sensor *h*s(*n*), as can be seen from Equation (2). Thus, the signals *x*ex(*n*) must be equalized either prior to the deformation due to the coil and sensor impulse response or before they are forwarded to the matched filter. In this contribution, the matched filter impulse responses are adapted to match with the modified transmitted signals, so that they are again comparable to the coil signals measured by the sensor. Each coil excitation signal is adjusted individually by the equalizer

$$h\_{\text{eq},i}(n) = \oint\_{\text{c},i} \hat{h}\_{\text{c},i}(n) \* \hat{h}\_{\text{s}}(n) \,. \tag{9}$$

with ˆ *<sup>h</sup>*c,*i*(*n*) and <sup>ˆ</sup> *h*s(*n*) denoting the approximated impulse responses of the coil *i* and the magnetoelectric sensor, respectively. The factor *g*ˆc,*<sup>i</sup>* describes the influence of other components in the measurement setup. This includes for example different gains of the individual coil amplifier channels and different conversion factors of the coils from current to magnetic field. This is, e.g., due to variances in the number of windings. These values are approximated by a constant for the considered frequency range. The values for the six coils and amplifier channels used in this study are given in Table 1.

**Table 1.** Parameter of the coils and amplifier channels used in this study. The conversion factors of the coils describing the relationship between current and magnetic field are described by the column *conversion*. The gains of the coil amplifier channels are normalized to the maximum value (channel 6). Both values are determined at the resonance frequency of the sensor.


The equalized signals are calculated according to

$$\mathfrak{x}\_{\mathsf{eq},i}(n) = \mathfrak{g}\_{\mathsf{eq},i} \underbrace{\mathfrak{x}\_{\mathsf{eq},i}(n) \* h\_{\mathsf{eq},i}(n)}\_{\mathfrak{x}\_{\mathsf{eq},i}(n)} \tag{10}$$

and forwarded to the matched filter. The weighting factor *<sup>g</sup>*eq,*<sup>i</sup>* = <sup>1</sup> *Ei* ensures that each equalized signal has the same correlation output value one between the signals *<sup>x</sup>*eq,*i*(*n*) and *<sup>x</sup>*˜eq,*i*(*n*) at lag zero. The factor *Ei* is the auto-correlation output of *<sup>x</sup>*˜eq,*i*(*n*) at lag zero. In Figure <sup>6</sup> the cross correlation of the signals *<sup>x</sup>*eq,*i*(*n*) and *<sup>x</sup>*˜eq,*i*(*n*) at time lag zero is visualized.

**Figure 6.** Cross correlation between individual equalized coil excitation signals.

It is obvious that adjacent coils signals have cross correlation values different than zero. This is due to the decay behavior of the sensor impulse response. Nevertheless, the shown values are still sufficient for separating the coil signals.

#### *4.2. Matched Filter*

As described in Equation (2), the input signal of the sensor is a superposition of the magnetic coil signals. Additionally, noise sources superimpose the signal. To obtain a high signal-to-noise ratio (SNR) the coil amplitudes can be increased, which leads to a high energy consumption. Alternatively, a matched filter [27] can be used, which increases the SNR and thus can perform well also with lower energy consumption. This additionally makes the algorithm more robust against distortions. Hence, to obtain the amplitudes of the coil signals measured by the magnetoelectric sensor, the sensor input signal *<sup>x</sup>*in(*n*) is matched filtered with the equalized coil excitation signals

$$\mathbf{x}\_{\rm mf,i}(n) = \mathbf{x}\_{\rm in}(n) \* \mathbf{x}\_{\rm eq,i}(L\_{\mathbf{r}} - n) \,. \tag{11}$$

The matched filter output can be evaluated every *L*<sup>r</sup> samples, since all coil signals have then been completely transmitted once

$$m\_i(k) \, = \, \text{x}\_{\text{mf},i}(kL\_\mathbf{r}) \,. \tag{12}$$

with the value *mi*(*k*) corresponding to the amplitude of the coil signal *<sup>i</sup>* at the sensor. These amplitudes are summarized in the vector *<sup>m</sup>*(*k*) = [*m*1(*k*), *<sup>m</sup>*2(*k*),..., *mN*<sup>c</sup> (*k*)] <sup>T</sup> and forwarded to the localization algorithm.

#### *4.3. Localization Algorithm*

For the estimation of the position and orientation of the sensor, the matched filter output vector *m*(*k*) is compared with the columns of the lead-field matrix *A*. As stated in Equation (4), each column *a<sup>j</sup>* describes the coil amplitudes that would be measured by the sensor (after being filtered by the matched filter), if the sensor would occupy the defined position and orientation pair described by *<sup>p</sup><sup>j</sup>* . To be more robust against gain uncertainties and to ensure comparability between the measured coil amplitudes and the lead-field matrix columns, the vectors are normalized to the respective absolute maximum value beforehand. The values for the cost function *<sup>c</sup>*(*k*)=[*c*1(*k*), ... , *cj*(*k*), ... , *cN*<sup>p</sup> (*k*)]<sup>T</sup> are calculated by [15]:

$$c\_{j}(k) = \sum\_{i=1}^{N\_{\mathcal{C}}} \left| \frac{a\_{ij}}{\max\left\{ |a\_{j}| \right\}} - \frac{m\_{i}(k)}{\max\left\{ |m(k)| \right\}} \right|^{2} \forall j \in \{1, \ldots, N\_{\mathcal{P}}\} \tag{13}$$

The estimated sensor position and orientation is then given by the position-orientation pair of the forward problem with the minimum cost function value

$$\mathfrak{p}(k) = \mathfrak{p}\_{l(k)} \text{ with } l(k) = \underset{j}{\text{argmin}} \ c\_{j}(k) \tag{14}$$

and can also only be determined every *L*<sup>r</sup> samples. It is obvious that localization errors occur due to the discretization of the localization area. If the sensor is not directly located on a defined grid point, the localization error will be at least the distance between the closest grid point and the sensor's location. Unfortunately, even higher localization errors can occur in some forward model configurations, due to the shape of the cost function *c*(*k*). Further information can be found in Appendix A. To overcome this problem a higher resolution is required. However, this will increase the computational complexity dramatically and hence can endanger the real-time capability of the system. By increasing the resolution iteratively, the localization can be performed with high accuracy and a moderate increase in computational complexity.

The flow chart for the iterative localization process is shown in Figure 7. The first iteration is calculated as described before. Instead of considering only one estimated position-orientation-pair as stated in Equation (14), the *N*<sup>b</sup> best position-orientation-pairs are taken into account. The minima and maxima of the included position-orientation-pairs plus a fraction (one step size) of the previous grid form the boundaries of the new grid. The new grid is again divided into *N*p position-orientation-pairs and the forward problem as described in Equation (4) is calculated. From then on, the steps are repeated as described in the upper part of this section. Grid refinement stops when either the resolution between

two adjacent grid points is less than a specified resolution or the maximum number of iterations *N*it has been reached. The estimated position-orientation pair is given in the last iteration as described by Equation (14).

**Figure 7.** Flow chart of the iterative localization approach for one time step *k*. As long as neither the maximum number of iterations nor the desired resolution is reached, the algorithm keeps refining the localization grid.

#### *4.4. Postprocessing*

To mitigate possible outliers in the localization results, a linear Kalman filter for smoothing the estimated localization outcomes is used. The measurement equation of the system can be described by

$$
\mathfrak{p}(k) = \mathrm{Hs}(k) + \mathfrak{n}\_{\mathrm{m}}(k), \tag{15}
$$

with the state variables *s*(*k*), the observation model *H* transforming the states into the measurement variables, and supposing white Gaussian measurement noise *n*m(*k*) [28]. Assuming a linear system, the state variables are updated via the equation

$$\mathbf{s}(k) = \mathbf{F}\mathbf{s}(k-1) + \mathfrak{n}\_{\mathbb{P}}(k) \,. \tag{16}$$

The matrix *F* is the state transition matrix and *n*p(*k*) is white noise with zero mean [28]. Due to the noise processes the measurement variables are subject to errors. The Kalman filter attempts to predict the states and thus reduces the influence of the noises stated in Equations (15) and (16). The calculations are performed according to the descriptions in [28]. Based on the *<sup>N</sup>*<sup>m</sup> <sup>=</sup> 6 measured variables summarized in *<sup>p</sup>*ˆ(*k*), there are *<sup>N</sup>*<sup>s</sup> <sup>=</sup> <sup>3</sup> · *<sup>N</sup>*<sup>m</sup> <sup>=</sup> <sup>18</sup> states available, when additionally considering the velocity and acceleration of the measured variables. The initialization of the variables and covariance matrices are described in Appendix B. The enhanced localization output is denoted as *p*ˆ enh(*k*). It is worth noting that the Kalman filter outputs should not lie outside the localization area and are thus restricted to the boundaries of the localization grid.

#### **5. Measurements and Results**

For the proof of principle, the measurements were performed in a two-dimensional space, i.e., only considering the x and y components of the sensor position and only considering the angle yaw *ψ* for the orientation of the sensor. The z component of the sensor's position, as well as the orientation angles roll *ϕ* and pitch *θ* are assumed to be zero. Nevertheless, the proposed method can easily be performed in the three-dimensional space without any restrictions. More coils should be used for this, positioned in the threedimensional space. Additionally, a smaller initial grid resolution might be used to keep the computational complexity low. The measurements were performed with a real-time system developed at the chair of Digital Signal Processing in Kiel [29]. A picture of the graphical user interface of the tool is shown in Figure 8.

**Figure 8.** Graphical user interface of the real-time system used for localizing the magnetic sensors. The estimated position and orientation of the sensor is shown graphically in the 3D view and as text in the lower left corner. The number of iterations *N*it, the number of considered position-orientation pairs for refining the localization grid *N*<sup>b</sup> and the desired resolution in position and orientation are adjustable during runtime.

> The parameter used for the measurements (as defined in the sections above) are listed in Table 2. The localization area is limited to values between 0 cm and 50 cm in x and y direction and between −90◦ and 90◦ for the yaw angle.


**Table 2.** Parameter of the coil signals used in this study.

The waiting time between two consecutive coil signals must be rather high due to the decay of the impulse response of the sensor. Consequently, a total time of *<sup>L</sup>*<sup>r</sup> *<sup>f</sup>*<sup>s</sup> <sup>=</sup> 896 ms, with *f*<sup>s</sup> being the sampling rate, is needed to completely transmit the signals and thus to generate one localization result. This time can be shortened tremendously when using a sensor with a higher bandwidth. The sensor is placed in fixed positions and with fixed orientations to determine the accuracy of the algorithm. The tested position-orientation pairs *<sup>p</sup>*s,*<sup>j</sup>* were chosen randomly and are shown in Figure 9. The arrow directions represent the sensor's long axis.

**Figure 9.** Position-orientation-pairs of the sensor used for testing the localization algorithm.

The tilt between the sensor's sensitive and long axis has been approximately determined by manually rotating the sensor in a Helmholtz coil. The tilt was measured outside a shielded environment and results in *γ* ≈ −45◦. However, the tilt of the sensor depends on various factors, such as the strength of the bias field [30] (e.g., the earth's magnetic field) and is therefore only an approximation.

The results of the localization for the position-orientation pair *<sup>p</sup>s*,3 over time are shown as an example in Figure 10. Here, *x*s,*j*, *y*s,*j*, and *ψ*s,*<sup>j</sup>* denote the x and y component and the yaw angle of the sensor at position *<sup>p</sup>s*,*<sup>j</sup>* , respectively. The localization output without the Kalman filter is described by *<sup>x</sup>*ˆs,*j*(*k*), *<sup>y</sup>*ˆs,*j*(*k*), and *<sup>ψ</sup>*ˆs,*j*(*k*) and after smoothing by the Kalman filter by *x*ˆenh s,*<sup>j</sup>* (*k*), *<sup>y</sup>*ˆenh s,*<sup>j</sup>* (*k*), and *<sup>ψ</sup>*ˆenh s,*<sup>j</sup>* (*k*).

**Figure 10.** Real and estimated position and orientation of *<sup>p</sup>*s,3 over time. The variances in the localization result are due to the presence of noise and cross talk in the measurement hardware.

There are some variances of the localization result over time. This is mainly due to the presence of noise, which leads to slightly varying amplitudes at the sensor and thus to small variations in the localization outcomes. The offset error might be due to cross talk between the coil amplifier channels and coupling of the magnetic coil signals into the cables and electronics. Additionally, it can be seen that the Kalman filter smooths the estimation output over time, so that outliers in the localization results do not have such a high impact.

To quantify the accuracy of the algorithm, a localization error is defined according to

$$\bar{\varepsilon}\_{\mathbf{r},j} = \frac{1}{N\_{\text{meas}}} \sum\_{k=0}^{N\_{\text{meas}}-1} \sqrt{\left(\hat{\mathfrak{x}}\_{\mathbf{s},j}^{\text{enh}}(k) - \mathfrak{x}\_{\mathbf{s},j}\right)^2 + \left(\hat{\mathfrak{y}}\_{\mathbf{s},j}^{\text{enh}}(k) - \mathfrak{y}\_{\mathbf{s},j}\right)^2} \tag{17}$$

for the position estimation and

$$\mathcal{E}\_{\rm d,j} = \frac{1}{N\_{\rm meas}} \sum\_{k=0}^{N\_{\rm meas}-1} \sqrt{\left(\hat{\Psi}\_{\rm s,j}^{\rm csh}(k) - \psi\_{\rm s,j}\right)^2} \tag{18}$$

for the orientation estimation. The number of localization outcomes is set to *N*meas = 50, according to a measurement time of 44.8 s. The accuracy of the localization results for all tested position-orientation-pairs is shown in Figure 11.

The localization error is lying between 1.12 cm and 6.90 cm for the position estimation and results in a mean error of about 3.44 cm. The error for the orientation estimation is between 3.02◦ and 16.76◦. This results in a mean error of about 11.23◦. When considering fixed positions and orientations of the sensor, the localization output can be averaged over time and compared to the real sensor position/orientation afterwards. When doing so, the localization accuracy can be slightly improved and results in values between 0.46 cm and 6.52 cm for the position estimation and 1.54◦ and 15.35◦ for the orientation estimation. The average error reduces to 3.14 cm and 10.54◦, respectively.

The high errors for the estimation of the sensor's position can result from the noise and the cross talk in the measurement system. Due to the different distances between the coils and the sensor the SNR is dependent on the sensor position/orientation. For example,

looking at position *<sup>p</sup>s*,9 an average SNR of about 9.5 dB is obtained at the sensor. Higher coil currents would lead to an improved SNR. Nevertheless, the goal was to localize with a minimum amount of energy. To avoid cross talk, the cables as well as the amplifier channels must be shielded. Additionally, the remaining cross talk can already be considered when setting up the forward problem or with an appropriate initial calibration. However, since the focus of this work is on the real-time localization pipeline and the calibration will be very extensive, it is not the subject of this work, but will be taken into account in our future work. The high error variance of the orientation estimation can partly be due to the change of the sensor's sensitive axis with respect to the bias field. Even a rotation in the earth's magnetic field can tilt the sensitive axis [30]. This problem only occurs with the sensors presented here and not with other types of magnetic field sensors. Furthermore, possible calibration errors do not only influence the position estimation but also the orientation estimation.

**Figure 11.** Mean localization errors for all tested position-orientation pairs. The index *j* depicts the respective position-orientation pair *<sup>p</sup>*s,*<sup>j</sup>* of the sensor as defined in Figure 9.

#### **6. Conclusions and Outlook**

An algorithmic pipeline for localizing magnetic sensors in real time was presented in this contribution. Besides the localization algorithm itself, pre- and postprocessing steps for an enhanced estimation of the sensor's position and orientation have been described. The potential of the proposed algorithm was emphasized by measurements with a magnetoelectric exchange bias ΔE-effect sensor. Nevertheless, the proposed method can be applied to any type of magnetic field sensor. Only the coil excitation signals must be adapted to the properties (frequency range, dynamic range, etc.) of the magnetic sensor used. Using the magnetoelectric sensor, a mean localization error of 3.44 cm has been reached. For the proof of principle the localization of the sensor has been limited to the two-dimensional space. Nevertheless, the localization can be easily extended to the three-dimensional space.

The achieved results are comparable with other magnetic position estimation approaches. In [4]a3 × 3 m grid was used for localizing, achieving an accuracy of less than 10 cm. Localizing with a 3D sensor in a grid size of 8 × 7 cm an accuracy of 2.6 mm could be reached in [31]. In [2] a localization accuracy of ≤ 2 mm has been achieved, using coils for the localization of the sensors in a flexible MEG system. That shows that our localization method performs well, but can still be improved. However, the localization method investigated in this contribution can be performed in real time and in parallel to magnetic measurements without any degradation. Additionally, the robustness in noisy environments is increased by the usage of a matched filter and the smoothing of the localization outcomes via the linear Kalman filter. Moreover, the usage of a magnetoelectric sensor can be advantageous with respect to later medical applications due to its small size and the low production and operation costs.

Until now, the biggest limiting factor is the hardware used. Moreover, only a simplified model of the magnetoelectric sensor has been used, where the sensor is reduced to a point model. A more detailed model of the sensor, which also considers the dimensions of the sensor as well as a bias field dependent tilting of the sensor's sensitive axis, could improve the localization accuracy. The results can be further improved using multiple sensors included in an array with fixed distances and orientations as described in [1]. To reduce the transmitting time of the coils and thus to increase the rate of localization outcomes, FDMA or CDMA approaches would be beneficial. This would require an adaption of the sensor hardware.

**Author Contributions:** Conceptualization, C.B. and G.S.; methodology, C.B. and G.S.; software, C.B.; validation, G.S.; formal analysis, C.B.; investigation, C.B.; writing—original draft preparation, C.B.; writing—review and editing, C.B. and G.S.; visualization, C.B.; supervision, G.S.; project administration, G.S.; funding acquisition, G.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was funded by the German Research Foundation (Deutsche Forschungsgemeinschaft, DFG) via the collaborative research center CRC 1261.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors would like to thank Alexander Teplyuk for building the coil amplifier and Christine Kirchhof, Phillip Durdaut, and Jens Reermann for building up the sensor.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **Abbreviations**

The following abbreviations are used in this manuscript:


#### **Appendix A. Localization Errors**

As already mentioned in Section 4.3, localization errors can occur due to the discretization of the localization area. It is obvious that the minimal localization error is defined between the real sensor position/orientation and the closest grid point and is only zero, when the sensor is lying directly on a defined grid point. Besides these obvious errors even higher localization errors can occur, if the grid resolution is not high enough. This highly depends on the coil arrangement, i.e., on the setup of the forward problem. In Figure A1 the cost function of a sensor located at point *A* (only considering a 2D plane and assuming a fixed orientation) is shown assuming an exemplary (not well set up) forward problem.

Now assuming a localization grid with a coarse resolution, which is visualized by the black lines. The points, where the lines are crossed are the potential grid positions. The minimal localization error would be the distance between point *A* and point *B*. However, the position-orientation pair with the smallest cost function is found at point *C*. This is due to the shape of the cost function and the small grid resolution, whereby the grid points do not lie on the minimum of the cost function. This results in a large localization error. Thus, a higher grid resolution is needed, which leads to a higher computational complexity and hence can endanger the real-time capability of the system. To overcome this, an iterative localization is implemented, which refines the grid resolution progressively.

**Figure A1.** Exemplary cost function. The sensor is located at point *A*. Due to the relatively coarse grid (black lines), there will be a localization error of at least the distance between the point *B* and point *A*. However, due to the shape of the cost function, the minimum that is crossing the grid lines—and thus the localization outcome—is at point *C*.

#### **Appendix B. Initialization of the Kalman Matrices**

As already stated in Section 4.4 the calculations of the Kalman filter are performed as described in [28]. The initialization of the matrices is filled as described for a discrete Wiener process acceleration model [28].

• State transition matrix

$$F = \begin{bmatrix} I\_{N\_m} & \Delta t I\_{N\_m} & \frac{1}{2} \Delta t^2 I\_{N\_m} \\ \mathbf{0}\_{N\_m \times N\_m} & I\_{N\_m} & \Delta t I\_{N\_m} \\ \mathbf{0}\_{N\_m \times N\_m} & \mathbf{0}\_{N\_m \times N\_m} & I\_{N\_m} \end{bmatrix} \tag{A1}$$

• Measurement matrix

$$H = \begin{bmatrix} I\_{N\_{\mathrm{m}}}, \mathbf{0}\_{N\_{\mathrm{m}} \times 2N\_{\mathrm{m}}} \end{bmatrix} \tag{A2}$$

• Covariance matrix of the process noise

$$\mathbf{Q} = \mathbb{E}\{\boldsymbol{\mathfrak{w}}\_{\mathbb{P}}\boldsymbol{\mathfrak{w}}\_{\mathbb{P}}^{\mathrm{T}}\} = \begin{bmatrix} \frac{\boldsymbol{\Lambda}^{4}}{4} \boldsymbol{I}\_{N\_{\mathrm{m}}} & \frac{\boldsymbol{\Lambda}t^{3}}{2} \boldsymbol{I}\_{N\_{\mathrm{m}}} & \frac{\boldsymbol{\Lambda}t^{2}}{2} \boldsymbol{I}\_{N\_{\mathrm{m}}}\\ \frac{\boldsymbol{\Lambda}t^{3}}{2} \boldsymbol{I}\_{N\_{\mathrm{m}}} & \frac{\boldsymbol{\Lambda}t^{2}}{2} \boldsymbol{I}\_{N\_{\mathrm{m}}} & \boldsymbol{\Lambda}t \boldsymbol{I}\_{N\_{\mathrm{m}}}\\ \frac{\boldsymbol{\Lambda}t^{2}}{2} \boldsymbol{I}\_{N\_{\mathrm{m}}} & \boldsymbol{\Lambda}t \boldsymbol{I}\_{N\_{\mathrm{m}}} & \boldsymbol{I}\_{N\_{\mathrm{m}}} \end{bmatrix} \boldsymbol{\mathcal{O}}\_{p}^{2} \tag{A.3}$$

• Covariance matrix of the measurement noise

$$\mathcal{R} = \mathbb{E}\{\mathfrak{n}\_{\rm m}\mathfrak{n}\_{\rm m}^{\rm T}\} = I\_{N\_{\rm m}}\mathfrak{d}\_{\rm m}^{2} \tag{A4}$$

• State covariance matrix

$$\mathbf{S}(0) = \mathbf{I}\_{3N\_m} \boldsymbol{\vartheta}\_{\mathbf{s}}^2 \tag{A5}$$

Here, *<sup>I</sup>M*<sup>1</sup> denotes a unit matrix of size *<sup>M</sup>*<sup>1</sup> <sup>×</sup> *<sup>M</sup>*<sup>1</sup> and **<sup>0</sup>***M*1×*M*<sup>2</sup> a matrix filled with zeros of size *<sup>M</sup>*<sup>1</sup> <sup>×</sup> *<sup>M</sup>*2. <sup>E</sup>{...} is the expected value operator and <sup>Δ</sup>*<sup>t</sup>* <sup>=</sup> *<sup>L</sup>*<sup>r</sup> *f*s , with *f*s being the sampling rate. Considering the high measurement noise and assuming a slowly moving or fixed sensor (i.e., low process noise), the estimated variances are set to *σ*ˆ <sup>2</sup> *<sup>p</sup>* <sup>=</sup> 0.001, *<sup>σ</sup>*<sup>ˆ</sup> <sup>2</sup> *<sup>m</sup>* <sup>=</sup> 10, and *σ*ˆ <sup>2</sup> *<sup>s</sup>* <sup>=</sup> 500. These values were chosen exemplary.

#### **References**

