A Multi-Node Magnetic Positioning System with a Distributed Data Acquisition Architecture

We present a short-range magnetic positioning system that can track in real-time both the position and attitude (i.e., the orientation of the principal axes of an object in space) of up to six moving nodes. Moving nodes are small solenoids coupled with a capacitor (resonant circuit) and supplied with an oscillating voltage. Active moving nodes are detected by measuring the voltage that they induce on a three-dimensional matrix of passive coils. Data on each receiving coil are acquired simultaneously by a distributed data-acquisition architecture. Then, they are sent to a computer that calculates the position and attitude of each moving node. The entire process is run in real-time: the system can perform 62 position and attitude measurements per second when tracking six nodes simultaneously and up to 124 measurements per second when tracking one node only. Different active nodes are identified using a frequency-division multiple access technique. The position and angular resolution of the system have been experimentally estimated by tracking active nodes along a reference trajectory traced by a robotic arm. The factors limiting the viability of upscaling the system with more than six active nodes are discussed.


Introduction
Magnetic positioning systems (MPS) are widely investigated as a suitable solution for applications needing localization in space of magnetic transmitters aptly mounted on objects or people [1]. The growing field of the Internet of Things (IoT) provides many application scenarios for localization systems [2], such as wireless sensor networks [3,4], mobile robots [5], and location-based services [6]. The main advantage of using localization systems based on a magnetic field is that they do not need line-of-sight conditions to work. For this reason, an extensively investigated application is the tracking of capsule probes inside the human body, i.e., wireless capsule endoscopy [7,8]. Non-line-of-sight conditions also occur in many indoor applications, whenever the localization of a magnetic transmitter inside buildings is required [9][10][11].
Short-range MPSs [12,13] can be used to accurately navigate surgical instruments [13,14], to investigate motor symptoms of diseases such as Parkinson's disease [15,16], to provide a human-machine interface such as in data gloves [17,18]. In the latter case, magnetic transmitters and/or sensors are mounted on the fingers in order to track their movements. This latter application would require tracking multiple magnetic nodes simultaneously, a topic that has been investigated in the last few years.
In [19], multiple magnetic markers are localized using Hall sensors. The magnetic field of the markers can be analytically calculated by modeling the markers as elementary magnetic dipoles. Measuring the magnetic field, marker positions are obtained by numerically inverting the model. The real-time simultaneous tracking of up to seven markers is reported. However, different markers cannot be identified, i.e., it is not possible to tell which marker is on each of the measured positions unless this information is known a priori. In [20], magnetic transmitters are coils supplied with a 119-kHz AC current. Signal separation and identification are obtained by modulating the coils with on-off-keying using orthogonal code sequences, an approach called Code Division Multiple Access (CDMA). A CDMA technique using Chirp Spread Spectrum signals has been recently investigated for an MPS with twenty magnetic transmitters [21]. In a different approach [22], each magnetic node transmits using a different frequency and signals are separated by Fourier analysis. This technique is called Frequency Division Multiple Access (FDMA).
Receiving coils are coupled with capacitors to form resonant inductor-capacitor (LC) circuits. Tuning transmitter frequencies on the resonance band of the receivers has a double advantage: first, the induced voltage is increased without increasing the absorbed power of the transmitter; second, the LC circuit acts as a filter, selecting only frequencies of interest within the resonance band, filtering wide-and narrow-band noise components.
In this work, we present a short-range MPS, including its performance results. Resonant coils are used as receivers. Also, the transmitters are small LC circuits supplied with an alternating current, thus generating a time-varying magnetic field. Multi-node tracking is based on the FDMA technique. The position and attitude of each node are measured (five degrees of freedom (DOFs) for each node; the attitude is defined as the direction of the solenoid axis). A novel distributed data-acquisition (DAQ) architecture has been implemented for fast performance. The system can perform up to 124 measurements per second of a single magnetic node and 62 measurements per second of six nodes. The position and attitude accuracy are a few millimiters and a few degrees, respectively. The applicability of the present system for measuring hand tremor associated with Parkinson's disease has already been proved by tracking a single node [16]. Multi-node tracking can be exploited in data gloves equipped with magnetic sensors, e.g., to assess hand kinematics and functioning [18,26], to drive industrial equipments by gesture recognition [27], or to devise automatic sign-language recognition systems [28].

Measurement Principle
Our MPS is devised to track small solenoids (TX) supplied with an oscillating current. When traversed by an electric current, each TX is equivalent to a magnetic dipole, according with Ampere's equivalence principle. Since an alternating current is supplied, each TX produces a varying magnetic field. This field induces on passive solenoids (RX) an electromotive force proportional to the time derivative of the flux of the magnetic field that crosses each RX coil surface, in accordance with Faraday's principle.

System Features
The system is based on a matrix of RXs mounted at fixed positions. By measuring the signal induced on each RX solenoid and by inverting a simple mathematical model based on the dipole approximation (i.e., each coil is treated like an elementary magnetic dipole), it is possible to extract the position and attitude of each mobile TX node.
The MPS presented in this work is an enhanced realization of the system presented in [12]. It includes 24 RX sensors: 16 sensors are mounted on the xy plane, while 4 sensors are mounted on both the xz and yz planes. On the horizontal plane, sensors are arranged on a regular 4 × 4 grid, with the distance between sensors being 8 cm. On vertical planes, sensors are positioned on the vertices of a square in which the sides are 8 cm, at the center of both planes. The active volume of the system is a 30 × 30 × 30 cm 3 box, but since a TX cannot be put too close to an RX because the signal magnitude would be too high, implying saturation problems, the effective active volume should be considered a 25 × 25 × 25 cm 3 box.
The supporting structure of the RX matrix is made of wood and has been cut with a numerically controlled milling machine ( Figure 1). RX solenoids are mounted into carved slots. The xy plane provides also a support for the robotic arm used to move the TXs along controlled trajectories. A matrix of slots was also carved on the xy plane to secure the robot on fixed positions and to easily place its end effector on fixed grid points as a reference to calibrate the movement of the robotic arm by performing the procedure described in [12].
The system can track up to six TX nodes simultaneously, while the system in [12] could track one node only. In order to discriminate different nodes, each node is supplied with a different voltage frequency, thus implementing the FDMA technique. The received signal on each RX is the sum of several frequency components that can be separated by Fourier analysis.
Signals on each RX are acquired by analog-to-digital converters (ADC). In [12], a single ADC was multiplexed to sample all channels, using a pseudo-simultaneous strategy: accordingly, the sampling frequency F s for each channel was equal to the sampling frequency of the ADC, divided by the number of channels. For the enhanced MPS described in this work, the idea is to use a distributed DAQ architecture completely parallelized [29] to sample all channels simultaneously. The signal on each channel is acquired by an independent ADC. Signal parameter extraction is performed by several microcontroller units (MCU) working in parallel. All extracted parameters are collected by a single MCU and then transferred to a computer that estimates the position and attitude of each TX. Figure 1. A picture of the system: (a) the whole magnetic positioning system with a robotic arm mounted on it. The reference frame is indicated. (b) A rear view of the yz plane: the round, carved slots for the receivers are visible; only four receivers are mounted on this plane; signals on the receivers are acquired by a microcontroller unit (the red board in the picture); and the microcontroller is mounted on a custom circuit board (the yellow board) described in Section 2.10. The other four receivers with a microcontroller unit are mounted on the xz plane; sixteen receivers are mounted under the xy plane with four microcontrollers. (The diagonal wire was used to wirelessly power the small solenoid (TX) in a different setup not treated in this paper.).

Principle of Operation
The active TX solenoid is treated as an elementary magnetic dipole m tx [22]: wheren tx is the versor orthogonal to the coil surface and defines what we call the attitude of the TX; r tx and S tx are its radius and coil surface, respectively; N tx is the number of windings; and I is the root mean square (RMS) value of the sinusoidal current of frequency f 0 flowing in the coil, for which the peak value is I 0 = I √ 2. The RMS magnetic field induced by a TX at the center of the ith RX receiver is as follows: where d i is the distance between TX and the ith RX, andn d,i is the unit vector associated with d i = r rx,i − r tx . By approximating B as being homogeneous over the whole RX coil surface, the RMS voltage induced on the ith RX is as follows: where N rx is the number of windings of each RX coil, S rx is its surface, andn rx,i is its unit vector assuming a left-handed orientation of the coil. If we collect position and attitude together as a single vector θ = [r txntx ] T , its value can be estimated by minimizing the following cost function: whereṼ rms,i is the measured voltages on each RX. Accordingly,θ = arg min θ F (θ).

System Calibration
An accurate knowledge of receiver positions and orientations is necessary for accurate tracking. A further source of error is the uncertainty on m tx [23] and on other quantities such as the RX surface S rx and the effective number of windings N rx . Moreover, each RX coil is coupled with parallel capacitance, thus forming a resonant circuit (see Section 2.5) with a voltage gain G Q (ω) varying with frequency. The signal on each RX is also amplified before being acquired, using an instrumentation amplifier (INA; see Section 2.10) with gain G ina . Hence, the measured voltage is G Q G ina V rms , and the uncertainties on G Q and G ina have to be taken into account.
A calibration procedure is thus needed in order to set as accurately as possible all these system parameters. Each active TX is moved along a reference trajectory, i.e., over a known values ofθ (k), where k indicates the time step. Then, the position and attitude of each RX, θ rx,i = [r rx,inrx,i ] T , are evaluated by minimizing the cost function: In (5), we have explicitly indicated also the calibration constants C rx,i defined in [12] that account for uncertainties on all other parameters listed above. Basically, since all the other parameters are multiplicative factors in the expression of V rms , the voltage actually measured is C rx,i G Q G ina V rms . By solving the minimization problem, C rx,i is also estimated, thus compensating for all the uncertainties on the system parameters. The calibration procedure is described with full detail in [12]. For both positioning and calibration, the nonlinear minimization problem can be solved using the concentrated cost function method introduced in [30], i.e., decomposing the problem into a linear and a reduced nonlinear problem. Once the system has been calibrated, as it is used to track a trajectory, the estimate ofθ (k) is made more accurate and stable by introducing a Kalman filter [31,32] based on the nearly constant velocity model [32]. The application of the concentrated cost function and Kalman filter to the MPS are illustrated in [12].

Coils and Alternating Voltage Supply
Both TX and RX solenoids are connected to a capacitor in parallel, so as to have resonant LC circuits. In this realization of the MPS, we have N tx = 36, S tx = π5 2 mm 2 , N rx = 252, and S rx = π9.5 2 mm 2 . All solenoids are made with a 33 AWG (American Wire Gauge) enameled copper wire. An alternating voltage is supplied to each TX by a signal generator. RXs are disconnected from any external supply, with the only source of voltage being the electromotive force induced by the oscillating magnetic fields of each TX. Both TX and RX solenoids are shown in Figure 2, together with their respective circuit models.

Figure 2.
A picture of the TX (a) and RX (b) solenoids with their circuit model: TX is supplied with voltage V gen from a waveform generator with an equivalent resistance R eq . On the RX solenoid, a voltage V Φ is induced by the external magnetic flux generated by TXs.

The TX Circuit
TX is supplied with voltage V gen generated by a waveform generator with an equivalent resistance R eq . The impedance of the TX is easily calculated as with the following amplitude: where ω 0 = (LC) − 1 2 is the resonance frequency when R = 0. If R is slightly larger than 0, the shift of the resonance frequency from ω 0 is negligible. Both V tx and V gen are directly measurable at different frequencies to obtain experimental values of |H tx (ω)|.

The RX Circuit
In the circuit model of an RX, there is no external supply. A voltage V Φ is induced across the coil by an external oscillating magnetic flux Φ = Φ 0 sin (ωt + φ). By Faraday's principle, V Φ = − dΦ dt = −ωΦ 0 cos (ωt + φ). Thus, the following applies to the RX circuit: with I = dQ dt . Substituting ansatz Q/C = V C = V C0 sin ωt into the equation, the frequency response is easily obtained: If the capacitor is removed, leaving the circuit open, i.e., C = 0, then |H rx (ω)| = 1 ⇒ V C0 = ωΦ 0 , i.e., the voltage measured across the inductor is equal to the derivative of the magnetic flux, as expected by Faraday's principle. It is possible to measure |H rx (ω)| keeping a TX and an RX at fixed positions for a set of ω values; a direct measure of the voltage across the RX solenoid without a capacitor gives ωΦ 0 , while measuring the same voltage with the capacitor gives V C0 . Their ratio is the experimental measurement of |H rx (ω)|.

System Frequency Response
A LCR (inductor-capacitor-resistance) meter at 200 kHz was used to measure C and L.
The measured values of resistance, capacitance, and inductance of the TXs are R eq = 50 Ω, R = 1 Ω, C = 50 nF, and L = 16.5 µH, and that for RXs are R = 6 Ω, C = 676 pF, and L = 1150 µH. When the voltages are measured using the oscilloscope, the input parasitic capacitance of the instrument (≈ 10 pF in our case) has to be added to C. The overall frequency response of the system can be estimated as H (ω) = H tx (ω) H rx (ω). Theoretical responses calculated using these parameters are shown in Figure 3, together with their experimental values for a TX and an RX. The theoretical curve for TX clearly matches the experimental one, while the theoretical curve for RX (the one indicated with "ideal" in the plot) is appreciably different from experimental values. This is due to the parasitic capacitance of the inductor and to the increase of resistance due to the skin and proximity effects as the frequency is increased. The varying magnetic field generated by the AC current induces eddy currents which cancel the current flow inside the conductor and concentrate the flow only in a layer near the surface; the depth of this layer is called skin depth. The proximity effect is instead due to eddy currents induced on a conductor winding by the varying magnetic field of the nearby windings. By including these effects as explained in the Appendix A, we obtained a much better theoretical curve (solid blue line in Figure 3).

Driving the TX
While designing the system, it is important to choose the frequencies of the signals supplied to each TX in order to match the high-response band of H(ω). Because of the device parameter tolerance, the actual responses of the solenoids are different. It is important to verify that the experimental response curves are largely superimposed. The chosen operating frequency band is delimited in Figure 3 by the green vertical lines. The alternating voltage is supplied to each TX through a square wave generated with a Cypress PSoC 5LP microcontroller (CY8CKIT-059 development board, Figure 4). The square wave has been used because it is easy to generate using the digital outputs (GPIO) of the PSoC. Each TX is connected to a different GPIO. The square wave component at fundamental frequency is selected by the narrow bandwidth of H (ω). The PSoC5 operates at low voltages (1.7 to 5.5 V); hence, it is suitable to be used with an autonomous power supply (e.g., a battery) as wearable electronics (e.g., data gloves). The amplitude of the square wave is controlled by regulating the voltage used to power the PSoC5. When several TXs are used instead of a single node, the signal received by each RX is the sum of the voltages induced by all TXs. Since, between each RX and the ADCs, an amplification stage is applied to measure small signals, this can cause saturation of the ADCs, in which the full scale is 3 V, when the TXs are too close to an RX. The solution we adopted is lowering the power voltage as the number of TXs increases. TX operating frequencies and PSoC5 power voltages are reported in Table 1; frequency values have been chosen to be as regularly spaced as possible given the clock divider values available on the PSoC5.  Table 1) is delimited by the vertical green lines.

Functional Scheme
The system is based on two main functional blocks, shown in Figure 5: the distributed DAQ, and the control and optimization software. The distributed DAQ is based on a network of MCUs communicating through a Serial Peripheral Interface (SPI), based on a master-slave architecture. We used a star configuration: the master-out-slave-in (MOSI), master-in-slave-out (MISO), and clock (CLK) lines from the master are split into six lines each, with each branch connected to a slave. There are six slaves. Each slave acquires signals on four RXs and performs the Fourier analysis to estimate voltage amplitudes and phases. The master triggers data acquisition on all slaves simultaneously and then collects all acquired data. There is also a Reset line (RST) to restart all slaves by operating only on the master MCU.The control and optimization software runs on a standard personal computer. This application communicates with the master MCU through a serial connection via the USB port. It triggers the master to start acquisition and then reads all voltage amplitudes and phases stored on the master. These data are used by the optimization routine to estimate the position and attitude of all TXs.
At the hardware level, the MOSI, MISO, and CLK lines of the SPI interface should be split from the master towards all slaves using a star topology. A short circuit simply made with cables and connectors could be the easiest solution, but we verified that it was not possible to connect more than two slaves. The fan-out of SPI ports on the Texas Instrument (TI) MCUs is not suitable to drive all the slaves, resulting in corrupted and practically unusable SPI signals. To solve this problem, we devised a custom splitter board using suitable bus drivers.
The splitter board ( Figure 6) has an input connector for power supplies and master SPI ports and six output connectors for the slaves. Power supply input pins are simply short-circuited with the corresponding output pins.
We mounted three-state buffers on the board. We chose TI SN74ABT125 bus drivers since their maximum output current is 32 mA (one of the highest values we found on the market) and their switching time is very short (≈5 ns satisfies the requirements if one considers that our SPI CLK frequency is 3.125 MHz). The output enabled by these buffers is driven by CS lines of the SPI. The MOSI line from the master is short-circuited to the inputs of six buffers, while each CS line is used to enable only one buffer at a time. Each buffered output is connected to a different slave. The same scheme has been implemented for the CLK line from the master. In this scheme, each slave is driven exactly by one buffer; hence, the fan-out is not a problem anymore. The scheme is reversed for the MISO lines. The six MISO lines from each of the slaves are connected to the inputs of six buffers, while their outputs are all short-circuited and connected to the MISO pin of the master. In Figure 6, all lines driven by a buffer are indicated in red. The 5 V pin on the input connector is used to power the buffers.

Control and Optimization Software
The control software has two distinct modes: free hand motion and robot motion; when in the first mode, all TXs are moved freely by an operator and the software estimates their positions. When operating in the second mode, the software also controls a robotic arm to move the TXs along predefined trajectories; at each acquisition, the actual positions are read from the robot internal MCU, thus allowing to compare the positions estimated by the MPS with the nominal or groundtruth ones. This can be done only after the robot has been calibrated using the reference grid, as explained in [12].
The used robotic arm is a Dobot Magician by Shenzhen Yuejiang Technology. Being made of aluminium and acrylonitrile butadiene styrene (ABS) polymer, it is not ferromagnetic, and in [12], we verified the it does not produce any measurable distortion of magnetic fields. The Dobot is controlled through a dedicated application programming interface (API) [33].
The control software is entirely developed in C++. A simplified flow chart of the software is shown in Figure 7. The robot is controlled by a dedicated thread running independently in order to move the robot and acquire data simultaneously. A predefined set of trajectory points is loaded, the robot thread is started, and then the program enters into its main loop. Through the serial port, the master MCU is triggered (which in turn triggers all slaves), and the application waits until data acquisition has finished and then reads all the voltage amplitudes and phases from the master. The optimization problem is solved by up to six concurrent threads (one for each TX). As in [12], we verified that six concurrent threads are faster than solving six optimization problems sequentially. The next data acquisition on the distributed MCU network is triggered immediately after the optimization threads have been started to perform concurrently both computation and acquisition for faster performances. Once the trajectory is completed (or the acquisition is stopped manually when running in the free hand mode), the main loop terminates. We tested the application on a computer equipped with a six-core Intel i7-8750H CPU at 2.20 GHz, running Windows 10.

FDMA and Band-Pass Sampling
As explained in Section 2.5, each TX solenoid is supplied with a different voltage signal frequency. The operating band is centered aroung f c = 182 kHz, and the bandwidth is B = 10 kHz. The separation between two near TX frequencies is about 2 kHz. In order to discriminate two near frequency components, the frequency resolution of the acquisition system has to be a few hundred hertz or better. The frequency resolution of the discrete Fourier transform (DFT) is ∆ f = f S N S , where f S is the sampling frequency and N S is the number of samples. Acquiring at Nyquist rate ( f S ≥ 2 × 182 kHz) would imply the acquisition of 2048 or more samples. Each MCU acquires signals on four RXs simultaneously, and 4 × 2048 samples exceed its memory capacity. In order to reduce N S , f S has to also be reduced below the Nyquist frequency. Since the bandwidth of TX signals is narrow, it is possible to apply the band-pass sampling technique [34]. If X ( f ) is the Fourier transform of a signal x (t), then the Fourier transform of the signal sampled with sampling rate f S = 1/T S is i.e., X S ( f ) is the periodic repetition, with period f S , of X ( f ). In Figure 8, the original spectrum for six TXs (blue) at arbitrary positions and a couple of undersampled signal spectra, with f S = 270 kHz (red) and f S = 150 kHz (green), are shown. It is important to choose f S such that spectral repetitions do not overlap. The nonoverlapping conditions are as follows: where m is an arbitrary integer indicating how many spectral repetitions there are in the interval [− f c , f c ], i.e., within the original spectrum. The optimal sampling frequency is the mean f S = ( f S + f S ) 2. In Figure 8, we used m = 1 (red) and m = 2 (green). The DFT returns only the part of the spectrum within the shadowed area (i.e., a single spectral repetion); hence, when using an odd m, the spectrum is flipped and care must be taken in order to correctly associate each peak with its correspondig TX. We chose f S = 270 kHz (m = 1) and N S = 1024 ⇒ ∆ f = 264 Hz. The memory capacity of MCUs supports N S = 1024. A noticeable difference of the present MPS compared to the DAQ system devised in [29] is that no fitting of the acquired signal is performed. Only the amplitude and phase of the signal are needed in the computation, and they are directly obtained from the spectrum calculated with the Fast Fourier Transform algorithm. Flat-top windowing is applied to the acquired data, optimized to give minimal scalloping loss [35].

Microcontroller Programming
Following the analysis in [29], we built the distributed DAQ architecture using Texas Instruments (TI) Delfino TMS320F28379D [36] microcontrollers, mounted on the TI LaunchXL development board. TI Delfino has a 200 MHz dual-core CPU and four 12-bit ADCs with a maximum sampling rate of 3.5 MSa/s. The MCU can be programmed using C++ and a dedicated compiler provided by TI. There are also two control law accelerators (CLA), i.e., additional co-processors that can execute simple operations in parallel with the main CPU. Other features relevant for our purposes are the direct memory access module (DMA), two SPI interfaces, and two serial communication interfaces (SCI; known as Universal Asynchronous Receiver-Transmitter (UART)). We used one SCI to connect the master MCU with the personal computer. The TI Delfino includes a floating point unit (FPU) co-processor for each core, supporting 32-bit single-precision floating point operations, that we used with the specific mathematical library provided by TI implementing a device-optimized FFT algorithm [37]; an interprocessor communication (IPC) module for communications between the two cores, that we used to control parallel operations; and a pulse width modulator (PWM), that we used to generate a square wave applied to trigger the ADCs at fixed sampling frequency. Figure 9 shows a simplified flow chart of both master and slaves operations. Once started, the master sends a set of configuration parameters (sampling frequency, number of TXs, TXs frequencies, and the bandpass parameter m) through SPI to all slaves to avoid reprogramming each slave whenever one needs to use different configurations. When the init phase has completed, the master enters a wait state until it receives an external trigger from the computer through the SCI. Immediately after, the master triggers all slaves simultaneously. Each slave uses its four independent ADCs to acquire signals on each RX. A DFT-based spectrum estimation is obtained for each RX. The maxima of its magnitude are found, and amplitudes and phases are stored in a vector. A cyclic redundancy check (CRC) [38] is applied, and the checksum is appended to the vector. After 6 ms (the fixed duration of signal acquisition and processing), the master reads through SPI the resulting vector from each slave and recalculates the checksum to verify that data have not been corrupted during transfer. If the checksum is correct, a time reference in µs, to be used in the Kalman Filter, is read from an internal timer of the master MCU. All data are sent to the computer through the SCI, and finally, the master returns in its initial state, waiting for an external trigger.

SPI Communication
SPI communication is controlled by the master. Slaves are programmed as finite-state machines; a global variable is used to define the state. On a slave, SPI communication is handled by a CLA to allow communication without interrupting the operations of the main CPU. When a word is received through the SPI, the CLA checks the state variable, executes the operations associated with the current state, and in case, changes the state.
In the initial state, the slave assumes the received words as configuration parameters and stores them in proper variables. As the last parameter has been received, the slave switches to the ready state and it can accept a trigger from the master. As signal acquisition has to be started simultaneously on all RXs, the master sends a trigger signal to all slaves at the same time, lowering all CS lines simultaneously. When triggered, a slave switches to the busy state and starts signal acquisition and processing. As all data are ready to be transferred, the slave switches to the data_ready state. In this phase, for any dummy word sent by the master, the slaves writes a data word on the MISO line until all data are transferred; then, the slave switches back to its ready state.

Signal Acquisition and Processing
As a slave is triggered, the PWM, by writing on the proper control register, is internally connected to the ADC's interrupt line to trigger signal acquisition at fixed sampling frequency f S . As a sample is acquired and converted, its value is written on the ADC result buffer and an interrupt signal is sent to the DMA that reads this buffer and writes the value on a vector cell, automatically increasing a pointer to the next cell for the subsequent sample. Once N S samples have been acquired, the PWM trigger is disconnected, the DMA pointer is reset, and the acquisition phase finishes.
There are four data vectors, one for each RX. The following sequence of operations is executed on each vector: the flat-top window function is applied to the data; the FFT is calculated; and the maxima of the magnitude spectrum are found: since TX frequencies are known a priori, maxima are directly found around their nominal values; hence peak finding is very fast. Since the CPU is dual-core, two sequences of operations are executed in parallel, halving the execution time. According to the FPU library documentation, for N S = 1024, the execution time of the FFT algorithm is 152 µs.

Slave Boards
The slave circuit is shown in Figure 10. Each slave MCU is connected to four RX solenoids. The signal of each RX has to be amplified to match the ADC input range. We implemented all these parts on a custom circuit board, one for each slave. The input connector of each board is wired to an output connector of the splitter board. A voltage regulator is mounted on the slave board, powered through the 6.5 V pin, and used to power the MCU with 3.3 V. The value 6.5 V from an external supply was chosen to compensate a slight voltage drop along the wires and the drop-out of the regulator. Each RX is connected to the inputs of an AD8421 instrumentation amplifier (INA) for which the nominal gain is set to 8. There are four INAs powered through the ±12 V connectors. Since the output signal of an RX is alternating, while the ADCs are monopolar, an offset has to be added to the signal to make it oscillate between 0 and 3 V, that is, the ADC voltage range. The offset is added through the V re f pin of the INA. It is obtained with a resistive divider connected to the output of the voltage regulator and connected to V re f through an op-amp buffer.

Testing the Magnetic Positioning System
To test the performance of the MPS, we set up six calibration and measurement stages, from a single TX to six TXs. For short terminology, we call one-node (or one-TX) configurations (or systems) an MPS tracking a single TX node and two-node, three-node, etc. configurations, an MPS tracking respectively two, three, etc. TX nodes simultaneously. Since, in different configurations, the TXs are supplied with different voltages and there can be mutual induction effects too, parameters resulting from the calibration in a particular configuration cannot be used with other configurations. A calibration procedure has to be performed for each configuration. The calibrated MPS is then used to track a reference trajectory in order to estimate the position and attitude errors. In Section 3.1, we describe the calibration and measurement operations, and in Section 3.2, the results will be presented.

Calibration and Measurement
In order to operate with all TXs mounted on the robotic arm, we 3D-printed the two supports shown in Figure 11 to be used as end effectors of the robot, fixing all TXs in the designated slots. The first support was used for calibration, since it keeps all TXs at the same attitude along the z axis to avoid the calibration accuracy varying among the TXs because of different attitudes. Defining the attitude in terms of elevation and azimuth angle, the Dobot arm, by construction, keeps the elevation fixed while moving, and only the azimuth changes along the trajectory. However, when the attitude is directed along the z axis (elevation 90 • ), it does not depend on the azimuth because of the cylindrical symmetry of this configuration.
For the measurement stage instead, we used the second support (attitudes are reported in the caption of Figure 11). As the robot moves, the elevations are kept constant at 90 • (n 2 ), 0 • (n 5 ), and 45 • (n 1 ,n 3 ,n 4 , andn 6 ), while the azimuth changes along the trajectory (except of course forn 2 ). With this support, we can estimate the position and angle accuracy of the MPS for three different elevations. Figure 11. TX node supports for calibration (left) and measurement (right): the robot reference frame is reported. In this reference frame, the attitudes of the second support, using direction cosines notation, Following the analysis of [12], the calibration trajectory was chosen to span as large a part as possible of the active volume of the MPS. The calibration trajectory is shown in Figure 12a. Of course, each TX on the support travels along a slightly different trajectory, and this figure shows only one out of six trajectories. The system is calibrated by minimizing the cost function (5), as explained in Section 2.3 and in [12]. All TXs are active while the signal is being acquired, as explained in Sections 2.8 and 2.9.2. There is a minimization problem for each TX, each one independent from the others, i.e., different sets ofθ rx,i andĈ rx,i values are obtained for each TX. After the calibration procedure was performed for all configurations, we proceeded with the measurement stages, mounting all TXs on the second support shown in Figure 11 and using the robot to move it along the trajectory of Figure 12b. Since the support has slots with different attitudes, we performed several measurements, testing each of the TXs on each and every slot to estimate the accuracy of each TX for different attitudes and positions. As an example, for the one-node configuration, we tracked a single TX six times, one for each slot; for the six-node configuration, we tracked six TXs six times, the first time mounting TX1-TX6 in slots 1-6; then in slots 2-6 and 1; then in slots 3-6, 1, and 2; and so on. The results are discussed in the next section.

Experimental Results
The first interesting result is that, although they are operated using different frequencies, the accuracy is practically the same for each TX, i.e., varying the frequency within the band defined in Section 2.8, with f c = 182 kHz and B = 10 kHz, does not produce any appreciable difference in the accuracy of the system. All the results presented in the following are thus an average over all the TXs of each configuration (see Table 1). The measurement rates for different configurations are reported in Table 2. We define the position accuracy σ r of the system as the Euclidean distance between the position measured with the MPS and the real position provided by the robot, while the angular accuracy σ θ is defined as the angle between the measured attituden and the real attitude n, i.e., σ θ = cos −1 n·n |n||n| . Measuring the accuracies for each trajectory point, we get distributions for both σ r and σ θ . The percentiles of σ r and σ θ are reported in Tables 3 and 4 for all configurations. The last row of both tables contains the percentiles calculated on the union of all distributions with respect to each configuration. The cumulative distribution functions (CDF) of σ r and σ θ for each configuration and for the union distribution are reported in Figures 13 and 14. The CDFs of different configurations are discernible, but the percentiles do not differ by more than ≈0.5 mm and ≈0.5 • , the accuracy of the system thus being practically the same for any configuration.   Since the arrangement of the RXs is not symmetrical (16 RXs on the xy plane and 4 RXs on both the xz and yz planes), we examined if the accuracy has some dependence from the attitude elevation angle. In fact, it has the best position accuracy obtained for the vertical attitude (90 • elevation) and the worst position accuracy for the horizontal attitude (0 • elevation); opposite results are obtained for the angle accuracy, the best values being obtained for the horizontal attitude and the worst being obtained for the vertical attitude. The percentiles for elevations 90 • , 45 • , and 0 • are reported in Tables 5 and 6. The corresponding CDFs are reported in Figures 15 and 16.   In order to understand the effect of mutual induction among the TXs and the effect of lowering the voltage supply as the number of TXs increases (see Section 2.5), we calculated the signal-to-noise ratio (SNR) for the 1-node configuration and the 6-nodes configuration. We define the signal as the ideal voltage V rms (k) calculated using (3) at each point of the reference trajectory traced by the robotic arm, while the noise is defined as the difference between the ideal and the measured voltagesṼ rms (k): In Figure 17, we report the two distributions, and in Table 7, we report their mean and standard deviation. The SNR is not so much affected by the voltage reduction and the mutual induction between nodes; this is consistent with the position and attitude errors that do not increase too much as the system is upscaled.

Conclusions
A magnetic positioning system capable of tracking the position and orientation of multiple coils in real time was presented and characterized. A median position accuracy of 3.1 mm and a median orientation accuracy of 2.4 degrees were observed with respect to reference measurements obtained using a robotic arm. The system can perform 62 measurements per second when tracking six nodes and up to 124 measurements per second when tracking one node only. The results show that the measurement accuracy of the system is very similar among different configurations. This is a clear hint that the system could be upscaled to more than six nodes. The most important limiting factor would be the saturation problem. Indeed, it is not possible to lower the transmitter supply voltage anymore because of the sensitivity limitations of the most distant RXs. This would result in a reduction of the effective number of receivers, compromising measurement accuracy. There are some possible solutions to overcome the saturation problem for future developments. One is to implement an automatic gain control on the INAs to dynamically reduce the gain as the signal amplitude approaches the saturation threshold. Another possibility is to increase the number of receivers, thus assembling them on a finer grid and then tuning different sets of TXs and RXs on different resonance bands. Lastly, ADCs with a greater range and a finer resolution could be used, e.g., 16-bit ADCs in place of 12-bit ADCs. This solution could not be implemented with the off-the-shelf components we used, since the 12-bit ADCs with fixed range are integrated in the microcontroller IC, but it is perfectly feasible in principle.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Inductor and Resistance Models for RX Solenoids
Since each RX coil has 252 closely packed windings of enameled wire arranged in 6 layers, it has an intrinsic parasitic capacitance C p that has to be taken into account. The parasitic capacitance is modeled adding a capacitor in parallel to the R − L series of Figure 2 (bottom). An equivalent impedance Z eq is obtained as Z −1 eq = (R + iωL) −1 + iωC p , that can be written using an equivalent resistance R eq and inductance L eq as Z eq = R eq + iωL eq : Using a bench LCR meter, we measured both L eq and R eq at 1, 10, 100, and 200 kHz, obtaining the values reported in Figures A1 and A2. The low-frequency values have to be used for L and R. According to (A1), experimental L eq can be fitted to obtain C p = 65 pF ( Figure A1, magenta). Using this value of C p to calculate R eq according to (A2), we obtain the yellow curve of Figure A2, which clearly does not match the experimental values. Figure A1. Measured inductance of an RX coil (blue) and its model given by Equation (A1) (magenta). Figure A2. Blue represents the measured resistance of an RX coil. Yellow represents R eq given by Equation (A2). Red represents R eq given by the Dowell model accounting for the skin and proximity effects, using Equation (A3). Magenta represents the Dowell model with corrected porosity factor.
In order to correctly model R eq in the AC regime, one has to take into account the skin effect and the proximity effect. In the first case, the varying magnetic field generated by the AC current induces eddy currents which cancel the current flow inside the conductor and concentrate the flow only in a layer near the surface; the depth of this layer, called skin depth, is given by δ ω = 2ρ µω, where ρ is the resistivity of the conductor and µ is its magnetic permeability; the second effect is instead due to eddy currents induced on a conductor winding by the varying magnetic field of the nearby windings. We can define the ratio F R = R AC /R DC of the AC resistance including skin and proximity effects over the DC resistance. An analytical model for F R was derived by Dowell [39]: where A = π 4 3 4 d δ ω √ η; d is the wire diameter; η = d p is the so-called porosity factor, accounting for the space between the windings; p is the distance between the centers of the adjacent winding conductors; and N l is the number of windings layers. The first term of the sum on the right-hand side of (A3) accounts for the skin effect, while the second term accounts for the proximity effect. Applying the Dowell model, i.e., substituting R with F R R in (A2), we obtained the red curve in Figure A2, which clearly overestimates the resistance of the inductor at higher frequencies. This behavior is due to the proximity term and in particular to the definition of the porosity factor. This is a known limit of the Dowell model, and usually, the porosity factor has to be fitted from finite elements simulations for different coil geometries, as in [40,41]. In the proximity term of (A3), we just treated η as a parameter that we determined by fitting the experimental resistance of Figure A2, thus obtaining a good model of R eq (magenta curve in Figure A2). We verified that the equivalent inductance L eq is not appreciably influenced by skin and proximity effects. Substituting R and L with R eq (ω) and L eq (ω) into (8), we obtained a theoretical response curve very similar to the measured one (blue curve in Figure 3).