2.1. Design of Two-Dimensional Magnetic Tracking System
The two-dimensional magnetic tracking system consists of a single-axis magnetic sensor and two sets of magnetic field gradient coils, as shown in
Figure 1. A single-axis magnetic sensor, which is an induction coil or a GMR spin-valve sensor, is used as the target in a two-axis gradient field. The excitation waveform for synchronous detection is provided by a function generator, model AFG-2225 from GW-Instek. The in-phase and quadrature excitation signals are amplified by a two-channel power amplifier to drive the gradient coils. The output signal of the sensor is analyzed by a lock-in amplifier, Model 7270 DSP from the Signal Recovery, to resolve the in-phase and quadrature components. The analog outputs of the lock-in amplifier are connected to a data acquisition (DAQ) device, model USB-6216 from the National Instruments. The graphical user interface (GUI) software for signal recording and real-time processing is coded in LabVIEW. To investigate the spatial response of the system, the gradient coils are mounted on an X-Y translation stage while the sensor is fixed above the center of the detection plane with a fixed gap, as shown in
Figure 1a. Two step-motors controlled by a PLC driver are used for driving the translation stage to move the gradient coil relative to the sensor.
The gradient magnetic field for position tracking is generated by the
x- and
y-gradient coils, as shown in
Figure 2. The two gradient coils have the same size and coil turns. Each gradient coil consists of four elements, including two main coils and two shim coils. The element coils are formed by wrapping 200 turns of enameled copper wires around rectangular frames made of non-magnetic acrylic plastic material. The shim coils are 300 mm × 50 mm × 3 mm in dimension, and the main coils are 300 mm × 115 mm × 3 mm in dimension. The four-element coils are mounted on an acrylic substrate plate. The shim coils for the
x-gradient field are fixed on the top of the plate, while the main coils are at the bottom. Moreover, the shim coils and main coils of the
y-gradient coil are reversely arranged to the
x-gradient coil. The four-element coils for each axis are connected in series, and the winding directions are arranged to generate the constant gradient field at the center of the tracking stage. The two gradient coils are mounted 20 mm apart and aligned orthogonally to each other.
The in-phase and quadrature-phase sinusoidal excitation currents Ii and Iq are injected into the x- and y-gradient coils, respectively. They are at the same frequency with the phases differed by 90°. The magnitude and polarity of gradient fields Bi and Bq are controlled by setting the amplitude and phase of the excitation current. The excitation signals from the function generator are amplified by an in-house-built power amplifier circuit to drive the coils with 1-kHz sine wave currents.
When the sensor moves on the detection plane above the gradient coils, its sensing direction is perpendicular to the plane. For the x-gradient coil, the distances from the detection plane to the shim and main coils are 66 mm and 71.5 mm, respectively.
For the
y-gradient coil, the distances are 50.5 mm and 45 mm for the shim and main coils, respectively. The optimal distances between the element coils are determined by numerical simulation to maximize the linear range of field distribution.
Bi and
Bq, the
z-axis field components for the
x- and
y-gradient coils, are calculated by calculating the contributions of the element coils, as follows:
where the subscript “
s” represents the shim coil and “
m” indicates the main coil. The
z-axis field component induced by each element coil is calculated by calculating the contribution of four current segments. For example, the field component for the +
x shim coil is:
where
Bz,s(
x1),
Bz,s(
x2),
Bz,s(
y1), and
Bz,s(
y2) are the contributions of AB, CD, BC and DA current segments, respectively, as shown in
Figure 3. For each segment, the
z-axis field component was calculated by the equation derived from the Biot–Savart law:
in which
and
The contribution of the other element coils, i.e.,
Bz,s(
x2),
Bz,s(
y1), and
Bz,s(
y2), were calculated in the similar way. The obtained field distributions for the
x- and
y-gradient coils on the detection plane are shown in
Figure 4. It was found that
Bi and
Bq increase monotonically along the axis of the gradient. The field
Bi of
x-coil is almost linear in the range of
x = ±0.1 m, while the field
Bq of
y-coil exhibits more nonlinearity. The shorter distance from the detection plane to the
y-coil causes deterioration in the uniformity of the gradient of
Bq but enhances the mean-field gradient. When a linearity-correction algorithm is applied, the usable dynamic range of position tracking for the
y-coil becomes larger than that for the
x-coil.
The field profiles generated by the two orthogonal gradient coils were analyzed numerically to achieve the maximum linear range, which corresponded to the detection area of the tracking system. The simulation of field distributions was performed using the ANSYS software with the excitation current of 4.28 mA. The simulated field distributions shown in
Figure 5.
Figure 5a,b illustrate the vector field distribution of
x- and
y-coil, respectively. From
Figure 4 and
Figure 5c,d, it can be observed that the calculated and simulated field distributions of
Bi and
Bq are similar in terms of the curve shape and amplitude.
To verify the simulation results, the field distribution generated by the combined gradient coils driven by a 1-kHz current was measured by an induction coil, as shown in
Figure 6. The working range on the
x–
y plane is 200 mm × 200 mm in dimension and the sampling points are spaced by 20 mm for both the
x- and
y-directions. The field distribution on the detection plane was obtained by recording the in-phase and quadrature-phase outputs from the sensor output, i.e.,
Vi and
Vq, respectively, using a dual-phase lock-in amplifier. It was found that the field distribution of
x-coil was more linear than that of the
y-coil, in agreement with the simulation results. The simulation results showed that, in the working range, the maximum amplitude was about 4.2 µT for the
Bq (
x-gradient) field and more than 6.3 µT for the
Bq (
y-gradient) field with the excitation current of 4.28 mA. The measured field distribution is shown in
Figure 6, where the maximum field amplitude in the working range is about 4.18 µT for
Bi and 6.8 µT for
Bq. The minor discrepancy in
Bq values between the simulated and experimental results may be attributed to the gain errors in the amplifying circuits for generating the driving currents in the
x- and
y-coils. The field distribution was also verified with the coils driven by a direct current using a commercial fluxgate sensor. Since the variation in the size, shape, and alignment of the coil may also induce the error in field magnitude, the comprehensive method for making the correction is experimental calibration.
2.2. Working Principle
The position of the target is obtained by detecting the gradient fields
Bi and
Bq generated by the two gradient magnetic coils, of which the coordinates are correlated to the field distribution. The sensor output, which contains the mixed signals from the
x- and
y-gradient coils, is analyzed by a digital dual-phase lock-in amplifier with a total gain of 50. The use of the phase-sensitive detection technique avoids interference from environmental disturbance, and hence optimizes the signal-to-noise ratios (SNR) [
21,
22]. The process for retrieving the position signals with dual-phase synchronous detection is illustrated in
Figure 7. The waveform of sensor output is multiplied by the in-phase and quadrature-phase reference signals synchronized to the excitation current with a phase difference φ. When φ is tuned to the correct value, the dual-phase outputs
Vi and
Vq correspond to the
x- and
y-coordinates of the target magnetic sensor, respectively. The response time of the system is determined by the low-pass filter with a time constant of 10 ms. The dual-phase outputs of the lock-in amplifier,
Vi and
Vq, are recorded by a 16-bit DAQ controlled by a GUI software. The linearity correction algorithm is performed in real time to retrieve the position of the sensor on the detection plane.
Let
Bi and
Bq be the in-phase and quadrature-phase amplitudes of the
z-axis component generated by the planar gradient coils, as shown in
Figure 2. In the linear range, they have constant gradients along the
x- and
y- axes, respectively. In the ideal case, their intensities are linear functions of the spatial coordinates
x and
y. When the origin of an
x-
y plane is set to be at the center of the gradient coils,
Bi and
Bq are proportional to
x and
y coordinates, respectively:
where
gx =
dBi/d
x and
gy =
dBq/d
y are the constant gradients of
Bi and
Bq. In our system, the waveforms of
Bi and
Bq are sinusoidal at the same frequency (1 kHz), but the phases differ by 90°. When a magnetic field sensor detects the
z-axis component
Bz at a field point on the same
x-y plane, it can be represented by
Bi and
Bq, as follows:
where
B0 = (
Bi2 +
Bq2)
1/2 is the amplitude and
ϕ = atan(
Bq/
Bi) is the phase of the
z-axis component. In terms of
B0 and
ϕ,
Bi =
B0 cosϕ and
Bq =
B0 sinϕ. The spatial coordinates
x and
y can be resolved from the sensor outputs using a dual-phase synchronous detection with a phase reference to the excitation currents in the coils. Given the field-to-voltage transfer coefficient
k for the field sensor, the orthogonal output voltages,
Vx and
Vy, correspond to the fields as follows:
In general, the output of the detection circuit has a reference phase difference
φ for
Bz. With a non-zero
φ, the corresponding in-phase and quadrature-phase outputs
Vi and
Vq are both functions of
φ:
where
ϕ = arctan(
Vy/
Vx). With the phase differences
φ and φ−90°, the in-phase and quadrature-phase outputs are associated with
Vx and
Vy, as follows:
This can be expressed in matrix representation:
Therefore, the coordinates
Vx and
Vy on the voltage plane are:
According to (18), the change in reference phase
φ corresponds to a rotation on the
Vx-
Vy voltage plane, which results in a rotation in the spatial coordinates
x and
y about the
z-axis. By adjusting the phase angle φ, one may redefine the orientation of the coordinate axes, providing that the gradients are perfectly constant. Using Equations (12) and (13) and (18), the spatial coordinates
x and
y can be calculated from the output voltage as follows:
For the gradients |gx| = |gy|, the aspect ratio between the x- and y- axes is 1, which means that the coordinates are not distorted under a rotation induced by changing the reference phase angle φ.
When the excitation currents were applied to the gradient coils, no crosstalk between the coils was observed since the currents were delivered by independent power amplifiers. As the phases of applied currents in the gradient coils were orthogonal to each other, the distortion in position tracking was dominated by the non-uniform field gradient, which can be corrected by the proposed algorithm.
2.3. The Linearity Correction Method
To correct the remaining non-uniformity of gradient field distribution, the linearity correction algorithm was applied to suppress the influence of distorted magnetic fields to improve the accuracy of position tracking. As an initial guess of the sensor’s position, the measured magnetic field was directly converted into the position based on the sensor outputs of in-phase
Vi and quadrature-phase
Vq using Equation (19). The linearity correction algorithm was performed by applying the cubic spline interpolation method to obtain smooth contour curves of the implicit function
Bz(
x,
y) = constant for both in-phase and quadrature-phase components. To implement the linearity correction algorithm, the databases of
Vi and
Vq were collected by stepping the sensor on the tracking stage with 121 sampling points. The sampling points were spaced by 20 mm for both the
x- and
y- directions. The
x- and
y- databases were expanded by using a cubic spline interpolation for the two-dimensional field distribution. Hence, the size of the interpolated databases was extended to 441 data points for
x- and
y- gradient magnetic fields. The interpolated field distributions of
x- and
y- coils consist of twenty-one natural cubic spline functions for
y- and
x- directions, respectively, as shown in
Figure 6. The interpolated curves of
x- and
y- gradient field distribution have a smoother profile, similar to the theoretically predicted distributions. Therefore, the gradient field distribution can be precisely mapped with the limited number of sampled field points.
The cubic spline functions for the
x- and
y-gradient field distributions are
Bi(
x, ym) and
Bq(
xm, y), respectively. The in-phase interpolating function
Bi(
x, ym) consists of the different cubic functions
Bin(
x) in twenty equally spaced intervals from
x1 to
x21:
where
Bin(
x)=
ain +
binx +
cinx2 +
dinx3 (
din ≠ 0) is a natural cubic function with
n = 1, 2, 3,..., 21. The function
Bi(
x, ym) passes through
Bi(
xn) =
Bin, where
n = 1, 3, 5,…, 21 are the measured data points and
n = 2, 4, 6,…, 20 are the interpolated data points. The coordinate
ym corresponds to the
m’th grid line,
m = 1, 2,..., 21, of the
x-gradient field distribution. Similarly, the quadrature-phase cubic spline interpolating function
Bq(
xm,
y) along
y-direction is:
where
Bqn(
y) =
aqn +
bqny +
cqny2 +
dqny3 is a natural cubic function with
n = 1, 2, 3,..., 21. The function
Bq(
xm,
y) passes through
Bq(
xn) =
Bqn, where
n = 1, 3, 5,…, 21 are the measured data points and
n = 2, 4, 6,…, 20 are the interpolated data points. The coordinate
xm corresponds to the
m’th grid line,
m = 1, 2,..., 21, of the
y-gradient field distribution. Suppose that the field components detected by the magnetic sensor at (
x0,
y0) are (
Bi0,
Bq0), and the in-phase and quadrature-phase components of sensor output are
Vi0 and
Vq0. The sensor position (
x0,
y0) can be retrieved by performing linear interpolations to calculate the intersection of two contour lines,
Bi(
x,
y) =
Bi0 and
Bq(
x,
y) =
Bq0 using the extended databases, as shown in
Figure 8.
At the sensor position (x0, y0), the corresponding field components are (Bi0, Bq0). To implement the linearity correction algorithm, the two-dimensional cubic spline interpolation was applied by using the inverse interpolating function x(Bi, ym) and y(Bq, xm). Since Bi and Bq increase monotonically along the axis of the gradient, the field components are compared to the values of the Bi and Bq interpolated databases, respectively, to find the two contour lines Bi(x,y) = Bi0 and Bq(x,y) = Bq0. For instance, to find the contour line Bi(x,y) = Bi0, comparison of Bi0 value to the Bi− interpolated database is performed. For each ym value (m = 1, 2, 3,..., 21), an interval (Bin−1, Bin) is found, satisfying a condition Bin−1 ≤ Bi0 ≤ Bin. The (Bin−1, Bin) intervals corresponding to the (xn−1, xn) interval satisfying a condition xn−1 ≤ x0 ≤ xn (n = 2, 3,...,21) are the landmarks in the database. Therefore, a natural cubic spline function in (xn−1, xn) interval is determined using Equation (20). The coordinates (x0m, ym) corresponding to Bi0 value is given by implementing the cubic spline interpolation with inverse function x(Bi, ym). The points (x0m, ym) forming a contour line of x− interpolated database are Bi(x,y) = Bi0. Similarly, the obtained points (xm, y0m) forming a contour line of y-interpolated database are Bq(x,y) = Bq0.
The position of the sensor is determined by finding an intersection of two contour lines. To find out the intersection point of two contour lines, the four nearest points around the intersection are determined by comparing the distance between the fixed interval points on x- and y- contour lines. The two y-contour points closest to the x-contour line are identified by obtaining two minimum distances, and the two nearest points on the x- and y- contour lines are identified in the same way. The four obtained points form two intersecting segments, of which the intersection point (x0, y0) gives the sensor position.