3D long-range through-the-wall magnetoquasistatic coupling and application to indoor position sensing

Abstract

Methods and systems for indoor position sensing are disclosed. The described methods and systems are based on magnetoquasistatic field coupling theory and can be implemented in two- and three-dimensional, long-range, through-the-wall applications, where the transmitting devices are implemented outdoor, the receiving device is implemented indoor, or vice versa. Measurement systems implemented to characterize the disclosed methods are also presented for both two- and three-dimensional applications involving indoor position sensing. Orientation sensing methods and systems are also disclosed.

Claims

1. A long-range, through-the-wall position sensing method comprising: providing one or more transmitting devices, each configured to transmit quasistatic magnetic fields along one or more transmitting axes; providing at least one receiving device configured to receive quasistatic magnetic fields along one or more receiving axes; placing the one or more transmitting devices outside an enclosed space having an outer wall; placing the at least one receiving device inside the enclosed space; transmitting through the outer wall, by the one or more transmitting devices, the quasistatic magnetic fields in one or more frequency bands; detecting, by the at least one receiving device, one or more quasistatic magnetic field couplings between each axis of each of the one or more transmitting axes and each axis of the one or more receiving axes; based on the detected quasistatic magnetic fields, calculating orientation-invariant ranges between the at least one receiving device and each of the one or more transmitting devices, the calculating comprising contributions from complex images of the one or more transmitting devices, thereby accounting for induced eddy-currents in ground due to transmitted quasistatic magnetic fields, and wherein the contributions from complex images are based on a dependency of the detected quasistatic magnetic fields on an inverse of distances of the complex images from the at least one receiving device to a third power.

2. The long-range, through-the-wall position sensing method of claim 1, wherein the one or more transmitting axes are three orthogonal transmitting axes, and the one or more receiver axes are three orthogonal receiving axes.

3. The long-range, through-the-wall position sensing method of claim 1, wherein the calculating the orientation-invariant ranges comprises calculating a sum function based on the square of the detected quasistatic magnetic fields.

4. The long-range, through-the-wall position sensing method of claim 3, further comprising inverting the sum function to calculate the orientation-invariant ranges.

5. The long-range, through-the-wall position sensing method of claim 4, further comprising calculating a ground truth position of the at least one receiving device based on the calculated orientation-invariant ranges.

6. A two-dimensional system comprising at least two transmitting devices placed outside an enclosed space and at least one receiving device placed inside the enclosed space, the system configured to operate according to a method comprising: providing one or more transmitting devices, each configured to transmit quasistatic magnetic fields along one or more transmitting axes; providing at least one receiving device configured to receive quasistatic magnetic fields along one or more receiving axes; placing the one or more transmitting devices outside an enclosed space having an outer wall; placing the at least one receiving device inside the enclosed space; transmitting through the outer wall, by the one or more transmitting devices, the quasistatic magnetic fields in one or more frequency bands; detecting, by the at least one receiving device, one or more quasistatic magnetic field couplings between each axis of each of the one or more transmitting axes and each axis of the one or more receiving axes, and based on the detected quasistatic magnetic fields, calculating orientation-invariant ranges between the at least one receiving device and each of the one or more transmitting devices, the calculating comprising contributions from complex images of the one or more transmitting devices, thereby accounting for induced eddy-currents in ground due to transmitted quasistatic magnetic fields, wherein: the at least two transmitting devices and the at least one receiving device are all at substantially same heights from ground, and the calculating the ground truth position of the at least one receiving device comprises trilateration.

7. The two-dimensional system of claim 6, wherein the at least two transmitting devices are each located at a distance of more than 10m from the outer wall.

8. The two-dimensional system of claim 7, wherein the at least two transmitting devices are positioned in a distance of at least 8m from each other.

9. The two-dimensional system of claim 8, wherein the enclosed space is a house with a size of larger than 8000 ft.sup.2.

10. A three-dimensional system comprising three or more transmitting devices placed outside an enclosed space and at least one receiving device placed inside the enclosed space, the system configured to operate according to method of claim 5.

11. The three-dimensional system of claim 10, wherein the three or more transmitting devices are four transmitting devices positioned in a rectangular configuration, and wherein a plane of the rectangle is parallel to the outer wall.

12. The three-dimensional system of claim 11, wherein the calculating the orientation-invariant ranges comprises calculating a sum function based on the square of the detected quasistatic magnetic fields, and inverting the sum function.

13. The three-dimensional system of claim 12, wherein the inverting the sum function is based on a non-linear algorithm to calculate the orientation-invariant ranges.

14. The three-dimensional system of claim 13, wherein the non-linear algorithm comprises one of a Levenberg-Marquardt algorithm, Gauss-Newton algorithm or gradient descent algorithm.

15. The three-dimensional system of claim 12, wherein the inverting the sum function is based on a linear method, and wherein the orientation-invariant ranges are calculated directly using powers received by the receiving device by each of the transmitting devices based on an analytical formula.

16. The three-dimensional system of claim 15, wherein a ground truth position of the at least one receiving device is directly calculated based on analytical formulae including the orientation-invariant ranges and a dimension of the rectangle.

17. The long-range, through-the-wall position sensing method of claim 1, wherein the one or more frequency bands include frequencies in the range of tens or hundreds of kHz.

18. The long-range, through-the-wall position sensing method of claim 5, wherein the at least one receiving device comprises one or more mobile receiving devices, and wherein the calculating the position is performed for the one or more mobile receiving devices.

19. The long-range, through-the-wall position sensing method of claim 1, wherein the one or more transmitting devices comprise one or more mobile transmitters, and wherein the at least one receiving device comprises two or more fixed receivers, the method further comprising calculating a location of the one or more mobile transmitting devices relative to the two or more fixed receivers.

20. The long-range, through-the-wall position sensing method of claim 1, wherein the one or more transmitting devices comprises one or more mobile transmitters, and wherein the at least one receiving device comprises two or more fixed receivers, the method further comprising calculating a location of the one or more mobile transmitters relative to a reference fixed receiver of the two or more fixed receivers, the reference fixed receiver being movable from a first reference location to a second reference location upon knowledge of a relative position change of the second reference location to the first reference location.

21. A long-range, through-the-wall position and orientation sensing method comprising: providing a transmitting device configured to transmit quasistatic magnetic fields along one or more transmitting axes; providing one or more receiving devices configured to receive quasistatic magnetic fields along one or more receiving axes of each of the one or more receiving devices; placing the transmitting device outside an enclosed space having an outer wall; placing the one or more receiving devices inside the enclosed space; transmitting through the outer wall, by the transmitting device, the quasistatic magnetic fields in one or more frequency bands; detecting, by each of the one or more receiving devices, one or more quasistatic magnetic field couplings between each axis of each of the transmitting device and each axis of the one or more receiving axes of the one or more receiving devices, and based on the detected quasistatic magnetic fields, calculating: ranges between the transmitting device and each of the one or more receiving devices; angular directions to each of the one or more receiving devices; orientation of each of the one or more receiving devices, and contributions from a complex image of the transmitting device, thereby accounting for induced eddy-currents in ground due to transmitted quasistatic magnetic fields, and wherein the contributions from complex images are based on a dependency of the detected quasistatic magnetic fields on an inverse of distances of the complex images from the at least one receiving device to a third power.

22. The long-range, through-the-wall position and orientation sensing method of claim 21, wherein the one or more transmitting axes are three orthogonal transmitting axes, and the one or more receiver axes of each of the one or more receiving devices are three orthogonal receiving axes.

23. The long-range, through-the-wall position and orientation sensing method of claim 21, wherein the calculating the ranges comprises calculating a sum function based on the square of the detected quasistatic magnetic fields.

24. The long-range, through-the-wall position and orientation sensing method of claim 23, further comprising inverting the sum function to calculate the ranges.

25. The long-range, through-the-wall position and orientation sensing method of claim 24, further comprising calculating a position of the receiving device based on the calculated ranges.

26. The long-range, through-the-wall position and orientation sensing method of claim 24, wherein the calculating of the angular directions and the orientations is based on ratio functions defined based on selected numbers of the detected quasistatic magnetic fields.

Description

DESCRIPTION OF THE DRAWINGS

(1) FIG. 1 shows an exemplary measurement arrangement comprising transmitting and receiving devices according to an embodiment of the present disclosure.

(2) FIG. 2 show exemplary transmitting and receiving devices according to an embodiment of the present disclosure.

(3) FIG. 3 shows an exemplary measurement arrangement according to an embodiment of the present disclosure.

(4) FIGS. 4A-4B show exemplary measurement results according to embodiments of the present disclosure.

(5) FIGS. 5A-5C show exemplary measurement results according to embodiments of the present disclosure.

(6) FIG. 6A shows an exemplary measurement arrangement according to an embodiment of the present disclosure.

(7) FIG. 6B shows an exemplary arrangement of transmitting devices according to an embodiment of the present disclosure.

(8) FIG. 6C shows an exemplary measurement arrangement according to an embodiment of the present disclosure.

(9) FIGS. 7A-7E show exemplary measurement results according to embodiments of the present disclosure.

(10) FIG. 8 shows a table summarizing exemplary measurements results according to embodiments of the present disclosure.

DETAILED DESCRIPTION

(11) Throughout the disclosure, the terms “long-range” and “through-the-wall” are used to describe the position sensing applications where the transmitting and receiving devices are set up outside and inside (or the other way around) of an enclosed area (such as a dwelling) and where the transmitting device is located at a distance larger than 10 m from an outer wall of the enclosed area, the outer wall being in between the transmitting and receiving device such that the transmitted fields travel through the wall to reach the receiving device.

(12) FIG. 1 shows an exemplary measurement arrangement (100) in accordance with an embodiment of the present disclosure. Measurement arrangement (100) comprises a transmitting device (Tx) located outside a dwelling (102) and a receiving device (Rx) positioned inside the dwelling. Transmitting device (Tx) is set up at a height (h.sub.m) above the ground and at a horizontal distance (d) from outer wall (101) of the dwelling (102). Receiving device (Rx) is set up at a height (h.sub.n) from the ground and inside the dwelling. Also shown are the orthogonal 3-axis MQS transmitter moments (m.sub.1, m.sub.2, m.sub.3) and the orthogonal 3-axis receiver moments (n.sub.1, n.sub.2, n.sub.3). The complex image theory moments (m.sub.1,i, m.sub.2,i, m.sub.3,i) corresponding to the transmitting device (Tx) are also shown in FIG. 1.

(13) With further reference to FIG. 1, transmitting device (Tx) comprises three orthogonal coils, each generating a unique magneto-quasi-static (MQS) field, and similarly, receiving device (Rx) comprises three orthogonal detection coils. The field sensed inside the dwelling is dominated by the source (shifted in location) and due to an exponential attenuation (or loss) that is only amplitude dependent as shown in [1]. The voltage measured at a small receive coil inside the dwelling is given by V=−jω[{right arrow over (B)}.Math.{circumflex over (n)}]a′ and by combining [3] and [1], the following is obtained:

(14) $\begin{array}{cc}{V}_T\approx -j\frac{{\mathrm{\omega \mu }}_{o}a\prime }{4\pi }\left[e^{-a\mathrm{\kappa \Delta }r}\frac{3\left(\hat{m}.Math.\hat{r}\right)\left(\hat{n}.Math.\hat{r}\right)-\hat{m}.Math.\hat{n}}{r^3}+e^{-a\mathrm{\kappa \Delta }{r}_i}\frac{3\left({\hat{m}}_i.Math.{\hat{r}}_i\right)\left(\hat{n}.Math.{\hat{r}}_i\right)-{\hat{m}}_i.Math.\hat{n}}{r_i^3}\right],& \left(1\right)\end{array}$
where {right arrow over (B)} is the vector magnetic field of an electrically-small coil driven by a current, {circumflex over (m)} and {circumflex over (n)} represents the surface normal unit vectors for transmitting and receiving devices (Tx, Rx) respectively, ω is the radial frequency, the subscript i denotes the complex image contributions for ground effects, α′ is the receive coil's area, |{right arrow over (m)}|=1 Am.sup.2, r and r.sub.i represent the distance from the transmitting device (Tx) and the complex image of the transmitting device (Tx) to the receiving device (Rx) respectively. The exponential function from [1] may be applied to both source and image terms, where ακ=α/δ.sub.b(1−j) is a scaling of the inverse complex skin-depth, where δ.sub.b=√{square root over (2/ωμσ.sub.b)}. The vector image range is approximately {right arrow over (r)}.sub.i=(h.sub.m+h.sub.i){circumflex over (z)}+ρ{circumflex over (ρ)}, where h.sub.i=h.sub.m+δ.sub.g(1−j), ρ represents the cylindrical coordinate, and δ.sub.g=√{square root over (2/ωμσ.sub.g)}. For the 2D case and for long-ranges (ρ>>h.sub.m+h.sub.i), the vertical dipole moments along the z-axis can be neglected due to image cancellation as shown in FIG. 1, and {circumflex over (m)}.sub.2-3.Math.{circumflex over (z)}={circumflex over (m)}.sub.2-3,i.Math.{circumflex over (z)}={circumflex over (n)}.sub.2-3.Math.{circumflex over (z)}=0, and {circumflex over (m)}.sub.2-3={circumflex over (m)}.sub.2-3,i [2].

(15) Continuing with the calculations above, the 2D long-range operation permits the simplification Δr.sub.i≈Δr. By using steps in [2], the total voltage can be then simplified to give:

(16) $\begin{array}{cc}{V}_T\approx C_g{e}^{-a\mathrm{\kappa \Delta }r}\left[3\left(\hat{m}.Math.\hat{\rho }\right)\left(\hat{n}.Math.\hat{\rho }\right)C_{\varphi }-\hat{m}.Math.\hat{n}\right]\times \left[\frac{1}{r^3}+\frac{1}{r_i^3}\right],& \left(2\right)\end{array}$
where C.sub.g is a field-coupling gain, C.sub.ϕ ≈(ρ/r).sup.2.fwdarw.1 at long-ranges [2], and (2) is valid for all dipole configurations. The complex shift (multiplier of internal distance) can be neglected due to internal ranges (Δr, Δr.sub.i) that are small relative to skin-depth (δ.sub.b) [1], so the total range, r and r.sub.i, are given by [2].
Orientation-Invariant Range Dependence

(17) The coupling in equation (2) above shows dependence on orientation, which may complicate the range-dependent field drop-off analysis. To remove this dependence and study the field power spreading/attenuation, the approach in [2] may be used to sum the squared components of all combinations of field couplings in FIG. 1. Since the z-components are negligible for long-ranges, as noted previously, the solution is given by:

(18) $\begin{array}{cc}\begin{array}{c}{F}_r=\underset{m_q,n_k}{.Math.}{V}^2\left(m_q,n_k\right)\approx C_{\mathrm{sg}}{e^{-\mathrm{\alpha \Delta }r}\left[\frac{1}{r^3}+\frac{1}{r_i^3}\right]}^2,\\ =C_{\mathrm{sg}}{e^{-\alpha H\left(r-d\right)}\left[\frac{1}{r^3}+\frac{1}{r_i^3}\right]}^2,\end{array}& \left(3\right)\end{array}$
where the steps in [2] are used, C.sub.sg=5C.sub.g.sup.2 is the system gain, and the complex attenuation constant is defined as α=2ακ. For the 2D case studied here, Δr≈H(ρ−d)≈H(r−d), since ρ≈r, where H(x) represents the Heaviside unit step function of x.
Description of Exemplary Electronics Implementations in Arrangement (100) of FIG. 1

(19) In what follows, exemplary electronic elements to implement the measurement arrangement (100) of FIG. 1 are described, although the person skilled in the art will understand that varieties of other electronic devices may also be envisaged when implementing such measurement arrangements in accordance with the teachings of the present disclosure.

(20) With reference to FIG. 1, measurement arrangement (100) may include two fixed-location tri-axial transmitting devices, a mobile receiving device, and a processing system (not shown). Each transmitting device may include three orthogonal square coils. Signaling for each coil is obtained from a Novatech DDS9m signal generator. These signals are fed into FET (field-effect transistor) gate drivers, powered at, for example, 24V, to drive the coils with high current. The coils in each axis are all resonated with cylindrical high-voltage mica capacitors that can withstand high currents without dissipating excess heat. A first transmitting device is made up of about 30 turns of 18-AWG copper magnet wire wound sparsely around a 32×32×32 cm cube made up of polycarbonate sheets. The first transmitting device is resonated and programmed to transmit MQS fields in the X, Y, and Z axes at 101.0, 92.4, and 70.55 kHz, respectively. A second transmitter may be identically constructed, however using a more portable antenna made up of 50 turns of 18-AWG copper magnet wire wound tightly around a 13.5×13.5×13.5 cm acrylic cube. The second transmitting device may be resonated and programmed to transmit in the X, Y, and Z axes at 76.8, 68.6, and 60.0 kHz, respectively. The receiving device was fabricated with about 2.5×5×0.5 inch in size, excluding the battery system.

(21) FIG. 2 shows a system block diagram (200) representing an exemplary combination of transmitting and receiving devices in accordance with an embodiment of the present disclosure. A set of three magnetically-loaded inductors may be utilized as pickup antennas to maximize portability. The signals received by each of the RX coils are sent through a low-noise pre-amplifier before being multiplexed/switched. The output is band-pass filtered. A series of low noise op-amps and differential amplifiers buffer the signal before it is digitized by a 16-bit ADC (analog-to-digital conversion) at a sampling rate of 1 MSPS. The receiver's DSP (digital signal processing) engine may be designed to scale with the number of transmitting devices. The DSP reads the ADC output of the current channel, mixes with an I/Q NCO (numerically controlled oscillator) at the specified frequency with a sum-and-dump filter followed by a CORDIC (factor combining) stage. After a time specified by the operator, the receiver switches to the next available channel, and re-starts the processing. The switch may be set to have a dwell time of 250 ms per channel and an integration time of 1 ms per sample. The result is a reduction in the output data rate (set by the length of the accumulator), providing magnitude solutions for each of the frequencies specified to the DSP engine. Further processing can be completed by a short FFT (Fast Fourier Transform) to correct for imbalance in the unsynchronized clock domains. As part of this study, the data was directly transferred to MATLAB® (a mathematical tool) via Xilinx's ILA (integrated logic analyzer) module, and processed using the sum-of-squares method after the unknown losses and gains in each axis/link are calibrated or removed.

(22) In what follows, some exemplary measurement results highlighting the performance of the disclosed methods and devices will be presented.

(23) Measurement Configuration

(24) FIG. 3 shows an exemplary measurement arrangement (300) according to further embodiments of the present disclosure. Two transmitting device (Tx1, Tx2) are used as part of the measurements in a larger than 8000 square feet house in an urban environment in West Los Angeles, Calif. In the exemplary embodiment of FIG. 3, transmitting device is located 11.1 m from the outer wall of the house, and the two transmitting devices (Tx1, Tx2) are located 8 m apart.

(25) With continued reference to FIG. 3, static (no motion) receive measurement locations throughout the interior of the house (302) are used. A total of 117 interior locations were surveyed relative to an external predefined origin collocated with transmitting device (Tx1). A floor plan was co-registered to the measurement grid and origin. FIG. 3 shows all surveyed locations (303) inside the house, along with overlay of the floor plan, and the locations of the two transmitting devices (Tx1, Tx2). At each surveyed location (303) in the interior, the receiving device is positioned almost horizontally and left to be static prior to data collection. The transmitting and receiving devices are both located at a height of about 2 m above ground. Data is retrieved-by and stored-in the laptop connected to the receiving device.

(26) Field Power Model Analysis

(27) In order to show the effectiveness of the disclosed model, data using the summed power of equation (3) for all points in the building was collected (at surveyed points (303) of FIG. 3), and graphed as a function of the relative ground truth range (range of each surveyed point from the transmitters). The result is shown in FIG. 4A for transmitting device (Tx1) of FIG. 3 and FIG. 4B for transmitting device (Tx2) of FIG. 3.

(28) The C.sub.sg system gain constant in equation (3) may be found based on a laboratory calibration of all system gains. The α (which is used to derive σ.sub.b of the house) and the δ.sub.g (which is used to derive σ.sub.g of the ground) constants are found using a numerical least-squares minimization of equation (3) to the data collected. The results show that the model agrees well with the data collected in the house. The constants derived give a ground and building conductivity of σ.sub.g≈9-16 mS/m and σ.sub.b<4 nS/m, respectively. The through-the-wall model with complex image theory (CIT) to account for ground effects as given by equation (3) shows an average improvement of 0.2-1.4 dB (mean) over the free-space theory from transmitting devices (Tx1, Tx2) of FIG. 3.

(29) 2D Position Sensing Results and Analysis

(30) Referring back to FIG. 3, by numerically inverting measured power between both transmitting devices (Tx1, Tx2) and the mobile receiving device for the respective ranges using equation (3), the 2D position of the device can be found by solving the 2D trilateration problem.

(31) FIG. 5A is a reproduction of FIG. 3 further showing, this time, the solved or inverted position at each point with the floor plan underlayed on the figure. Both the actual surveyed points (303) and the corresponding location is estimated by the model (equation (3)) are shown. The transmitters locations are also noted in the figure. The estimation error for each point is represented by a line connecting each surveyed point (303) to its corresponding estimation (303′).

(32) FIG. 5B shows the cumulative error distribution of all points shown in FIG. 5A. The mean and RMS (root-mean-squared) error were found to be 0.53 m and 0.76 m for all y-locations and 0.75 m and 1.08 m for all x-locations. About 90% of measurements show error below 0.98 m and 1.4 m for y-locations and x-locations, respectively. However, the worst-case peak error is found to be 3.5 m and 5.4 m for y-locations and x-locations, respectively

(33) With further reference to FIG. 5A, in order study repeatability, the mobile receiver is located at 25 different locations inside the house (sub-set of surveyed locations (303)) and 20 repeat measurements were conducted at each location, which are all inverted and studied for repeatability or biases.

(34) FIG. 5C shows the result of the repeatability measurements. It is observed that there are some small amounts of variability. Nevertheless, in most cases, the variability is small or insignificant relative to the average error found in FIG. 5C. Based on this analysis and field observations, it is found that the largest sources of error is due to measurement locations near large metallic household appliances in the kitchen (refrigerator) and laundry (washing/drying machines) areas, that support strong eddy-currents and strong secondary scattered magnetic fields.

(35) 3D Orientation-Invariant Range Dependence

(36) FIG. 6A shows a measurement arrangement (600A) according to embodiments of the present disclosure. The measurement arrangement (600A) is overall similar to the one described with regards to FIG. 1, except that there is a 2×2 array of transmitting devices (Tx1, . . . , Tx4) are used outdoor for the arrangement (600A) of FIG. 6 to enable a more accurate 3D position sensing. Transmitting devices (Tx1, . . . , Tx4) are set up at different heights and the goal here is to analyze the results in three dimensions (3D) taking into account the height above the ground. In what follows, a model for this scenario is derived and some measurement results highlighting the performance of such model is described in detail.

(37) For 3D applications, in order to derive an orientation-invariant range dependent formula, some conversions and simplifications to Equation (2) may be applied. For simplicity, it is assumed [3].sup.2={1,2,3}×{1,2,3} and c=−jωμ.sub.0a/4π.Math.e.sup.−ακΔr. With reference to the embodiments of FIGS. 1 and 3, it is assumed that the distance between the transmitting and receiving devices is far enough to allow e.sup.−ακΔr≈e.sup.−ακΔri. Let {circumflex over (m)}.sub.1, {circumflex over (m)}.sub.2, {circumflex over (m)}.sub.3 be the unit vectors in the direction of the transmitter's magnetic moments. Similarly, {circumflex over (n)}.sub.1, {circumflex over (n)}.sub.2, {circumflex over (n)}.sub.3 are the unit vectors in the direction of the receiver's magnetic moments. Observe both {circumflex over (n)}.sub.1, {circumflex over (n)}.sub.2, {circumflex over (n)}.sub.3 and {circumflex over (m)}.sub.1, {circumflex over (m)}.sub.2, {circumflex over (m)}.sub.3 form rotated permutations of the standard basis. Such a basis preserves the magnitude and angle between all vectors of R.sup.3, i.e. we have v.Math.u=({circumflex over (n)}.sub.1.Math.v, {circumflex over (n)}.sub.2.Math.v, {circumflex over (n)}.sub.3.Math.v).Math.({circumflex over (n)}.sub.1.Math.u, {circumflex over (n)}.sub.2.Math.u, {circumflex over (n)}.sub.3.Math.u) for all v,u∈R.sup.3. Based on such observation, the following can be derived:

(38) $\begin{array}{cc}{V}_{m.fwdarw.n}^2=\frac{c^2}{r^6}\underset{\forall \left(i,j\right)\in {\left[3\right]}^2}{.Math.}{\left(3\left({\hat{m}}_i.Math.\hat{r}\right)\left({\hat{n}}_j.Math.\hat{r}\right)-\left({\hat{m}}_i.Math.{\hat{n}}_j\right)\right)}^2=\frac{6c^2}{r^6}& \left(4\right)\end{array}$
Similar calculations show V.sub.m.sub.i.sub..fwdarw.n.sup.2=6c.sup.2/r.sub.i.sup.6. Power measured by the receiver is given by:

(39) $\begin{array}{cc}P\stackrel{.}{=\underset{\forall \left(i,j\right)\in {\left[3\right]}^2}{.Math.}}{\left(V_{m_i.fwdarw.n_j}+V_{m_i^i.fwdarw.n_j}\right)}^2=\frac{6c^2}{r^6}+\frac{6c^2}{r_i^6}+2\underset{\forall \left(i,j\right)\in {\left[3\right]}^2}{.Math.}{V}_{m_i.fwdarw.n_j}{V}_{m_i^i.fwdarw.n_j}& \left(5\right)\end{array}$

(40) Let θ be the elevation angle of {right arrow over (r)}, i.e. the angle between {right arrow over (r)} and the xy-plane, similarly let θ.sub.i be the elevation angle of {right arrow over (r)}.sub.i. Geometry shows that r cos θ=r.sub.i cos θ.sub.i. For further simplification a temporary notation is introduced as f(v)=(1, 1, −1).Math.v. Using this notation, it can be observed that {right arrow over (r)}.Math.f({right arrow over (r)})=cos(2θ) and {right arrow over (r)}.sub.i.Math.f({right arrow over (r)}.sub.i)=cos(2θ.sub.i). Because both {circumflex over (r)} and {circumflex over (r)}.sub.i point towards the receiver, we can deduce ({right arrow over (r)}.Math.{right arrow over (r)}.sub.i)({right arrow over (r)}.Math.f({right arrow over (r)}.sub.i))=cos(θ.sub.i+θ)cos(θ.sub.i−θ). Hence,

(41) $\begin{array}{cc}\stackrel{\prime }{\underset{\forall \left(i,j\right)\in {\left[3\right]}^2}{.Math.}}{V}_{m_i.fwdarw.n_j}{V}_{m_i^i.fwdarw.n_j}=c^2\frac{3\left({\cos }^2\theta +{\cos }^2{\theta }_i-1\right)+1}{r^3{r}_i^3}\frac{P}{c^2}=\frac{6}{r^6}+\frac{6}{r_i^6}+\frac{6\left(1+r^2/r_i^2\right){\cos }^2\theta -4}{r^3{r}_i^3}.& \left(6\right)\end{array}$
Description of Measurement Configuration

(42) An experiment was conducted at the same >8000 square foot residential building in West Los Angeles (as described with regards to FIG. 3), but on the back side of the house. The same transmitting device driver electronics used in the embodiment described previously with regards to FIG. 3 are used here to drive high-current signals into the antennas. However, the antenna coils were built by wounding 12 AWG Magnet wire 30-40 turns around a 2″×12″×12″ cube made out of ¼″-thick acrylic sheets. With these new antennas, the transmitting devices can generate MQS fields with a maximum range of 75 m when powered at 30V. Each node transmitted three unique MQS fields from its orthogonal coils, with frequencies ranging from 55.9 to 171.6 kHz. The mobile receiving device antenna was also modified to limit the self-coupling between the three receiving device coils. This was achieved by assembling 10 planar ferrite cores together and wounding about 150 turns of 32 AWG magnet wire around each side. By integrating the three coils together in a cube-like fashion, the coils were completely orthogonal to each other, greatly reducing the mutual inductance between the loops (up to 12 dB improvement).

(43) Continuing with the 3D experiment, four tri-axial transmitter nodes were fixed in a vertical rectangular configuration, as shown in FIG. 6B, on a trailer positioned about 5 m away from the back side of the house. Two transmitters were raised up 7.5 m above the ground using the trailer's antenna masts. The remaining two transmitters were placed directly underneath on the trailer beam at about 1 m above the ground. The horizontal baseline distance between the transmitters on the masts was roughly 9 m, while the distance between the transmitters on the beams is 8 m. According to embodiments of the present disclosure, the rectangular configuration of the transmitters may be in a plane parallel to the outside wall of the residential building.

(44) FIG. 6C shows measurement arrangement (600C) for the 3D experiment in accordance with embodiments of the present disclosure. Each transmitting device position was surveyed using a Leica Flexline TS06 Laser Total Station to precisely measure their relative 3D locations from the Leica's reference point. On each of the three floors of the house, a visible reference point was marked and measured to also determine the 3D position relative to the Leica. From these floor reference points, 270 interior locations (603) spaced equally by 1 m were carefully measured and marked on all three floors of the house. With these reference measurements collected by the Leica, the transmitting devices locations and marked locations inside the house were all translated to the same coordinate system to confirm ground truths. FIG. 6C shows the ground truth locations of the four outdoor transmitter nodes on the trailer and all locations (603) on the 3.sup.rd floor of the house where data was collected. At each of the 270 interior locations (603), the receiving device sampled and processed the field measurements of all four transmitting devices. The data was immediately streamed to the connected laptop and saved using MATLAB®.

(45) Linear and Non-Linear Inversion Approaches

(46) Due to the complexities of the full equation that characterizes the MQS fields in 3-dimensions, the range inversion can only be solved nonlinearly when high accuracy is needed via careful consideration of precise locations of the reference nodes. The Levenberg-Marquardt Algorithm may be used for this purpose. A gain constant C.sub.sg and a ground depth constant α for each transmitter and rewrote Equation (2) in terms of x, y, and z may be introduced.

(47) $\mathrm{Let}\rho =\sqrt{\Delta x^2+\Delta y^2},\theta ={\tan }^{-1}\left(\frac{\Delta z-h}{\rho }\right),r=\sqrt{{\rho }^2+{\left(\Delta z-h\right)}^2},\mathrm{and}{r}_i=\sqrt{{\rho }^2+{\left(\Delta z+h+\alpha \right)}^2}.$
Then the following can obtained:

(48) $\begin{array}{cc}P=C_{\mathrm{sg}}\left[\frac{1}{r^6}+\frac{1}{r_i^6}+\frac{\left(1+r^2/r_i^2\right){\cos }^2\theta -2/3}{r^3{r}_i^3}\right]& \left(7\right)\end{array}$

(49) With reference to equation (7), in accordance with the teachings of the present disclosure, using the measured powers from all four transmitter nodes and their stationary locations as the inputs, the Levenberg-Marquardt Algorithm is able to directly optimize for and output a precise position of the receiver.

(50) To obtain a linear version of the problem we simplify the problem by allowing r.sub.i≈√{square root over (r.sup.2+α(α+2h).sup.2)} and therefore the following can be obtained:

(51) $\begin{array}{cc}P\approx C_{\mathrm{sg}}\left[\frac{1}{r^3}+\frac{1}{r_i^3}\right],& \left(8\right)\end{array}$

(52) Equation (8) is an under-approximation, but it is simplified to a function of just one variable r. Using four transmitters in a rectangular configuration, it is possible to obtain coordinates of the receiver. It is simple to write equations for r.sub.1, r.sub.2, r.sub.3, and r.sub.4 in terms of x, y, and z. Using only simple algebraic manipulation, one can derive the equations below:

(53) 0$\begin{array}{cc}z=\frac{1}{4h}\left(r_1^2+r_2^2-r_3^2-r_4^2+2h^2\right),x=\frac{1}{2D}\left(D^2+r_1^2-r_2^2\right),y=.Math.\sqrt{r_1^2-x^2-z^2}.Math.& \left(9\right)\end{array}$
Here, D is the distance between transmitter 1 and transmitter 2. Data was gathered to test the accuracy of the Levenberg-Marquardt Algorithm and the linear approximation. In what follows, the results are described in detail.
Field Power and Range Dependence Analysis

(54) With reference to FIG. 6C, at each surveyed location (603), the receiving device measured the three unique fields generated per transmitting device in each of its own orthogonal coils, resulting in nine total power components saved for each respective transmitting device. This 3×3 matrix of field power measurements is summed together and plotted as a function of the actual distance from the transmitting node as shown in FIG. 7A. Using a weighted least-squares optimization solver on the collected data, the constants C.sub.sg and α for in equation (7) are determined for each transmitter and are listed below:
C1=1, α1=3
C2=1, α2=6.1578
C3=1, α3=11.8503
C4=1, α4=14.4662

(55) Each solid curve in FIG. 7A represents the range inversion from the powers measured using the optimized constants on each respective floor. Measured values vs distance and corresponding to each solid curve are also shown. Both values obtained from the model and the measured values correspond to transmitting device (Tx1). The full equation accounts for the height from the ground and is able to closely characterize the fields using the same constants irrespective of which floor it is on, as shown by the plotted results in FIG. 7A.

(56) FIG. 7B shows the errors in range inversion at various distances inside the house. It can be noticed that the measurements on the first floor, have the highest range error compared to the second and third floors. This was consistent across all four transmitting devices, which may be due to metal appliances and other sources of environmental perturbations located on the first floor.

(57) 3D Position Sensing Results and Analysis

(58) With the simplified linear approximation (equations 8-9), the range from each transmitting device can be directly calculated. Due to the approximations made to remove the height dependence, the gain constants need to be re-optimized to account for the new generalization model. From the data collected on all three floors, the optimized constants for the simplified model (using the same α ground constants) are listed below:
C1=0.44787, α1=3
C2=0.50837, α2=6.16
C3=0.55996, α3=11.85
C4=0.60395, α4=14.47

(59) Inserting these constants into Equation 8, power measurements can be directly inverted into ranges from each transmitter. Knowing the rectangular configuration of the transmitters, simple trilateration and trigonometric techniques are used to solve for the relative X, Y, Z positions of the receiving device, with respect to an external reference point (in this case transmitting device (Tx1). Using the predetermined height difference, h, between transmitting devices (Tx1, Tx3) and the calculated range from each transmitter (r.sub.1, r.sub.3), the Z-position is calculated using the law of cosines. The same process is completed with transmitting devices (Tx2, Tx4) and averaged with the previous height result. The resulting height sensing equation simplifies to:

(60) $\begin{array}{cc}z=\frac{r_1^2+r_2^2-r_3^2-r_4^2+2h^2}{4h}& \left(10\right)\end{array}$
With the solved Z-position and ranges from the horizontal transmitting devices (Tx1, Tx2), the 2-D X and Y positions can be solved with these equations:

(61) $\begin{array}{cc}x=\frac{D^2+r_1^2-r_2^2}{2D},y=.Math.\sqrt{.Math.r_1^2-x^2-z^2.Math.}.Math.& \left(11\right)\end{array}$

(62) FIG. 7C shows the height calculated from using both the Simplified and Nonlinear techniques. From the plot it can be noticed that the errors on the first floor consistently has the highest errors, which can also be verified from the power measurements of transmitting device (Tx1) in FIG. 7B. There are instances on the first floor where the Z error can be significant enough to produce ambiguous floor sensing. For example, the maximum Z-position error reaches 2.895 m which could overlap the 2.sup.nd floor threshold and produce an incorrect floor result. Since the 2D position tracking depends on the height calculated, these height sensing errors can propagate into the position inversion results.

(63) FIGS. 7D-7E show the inverted 2D position results between for simplified linear and nonlinear solutions for the first and second floor respectively. The actual ground-truth position of the measurement collected is also shown. FIG. 8 shows table (800) summarizing the average errors in each dimension for each of the techniques (simplified and non-linear) used. Analyzing the overall errors of both approaches, it can be generalized that the non-linear has improved X and Y errors but at the expense of adding Z error. The simplified model works well in directly solving for 3D position using a closed-form solution with comparable results as the nonlinear solver.

(64) According to various embodiments of the present disclosure: One or more receiving devices and one or more transmitting devices may be used in any of the systems and measurement arrangements described throughout the present disclosure. The one or more receiving devices may comprise one or more mobile and/or one or more fixed receiving devices. The one or more transmitting devices may comprise one or more mobile and/or one or more fixed transmitting devices. The position of the one or more transmitting devices may be calculated with reference to a fixed receiver of the one or more receiving devices. The fixed receiver used as the reference may be relocated from one reference location to another as long as the relative positions of the reference locations with respect to one another are known. The position of the one or more transmitting devices may be calculated with reference to the location of two or more receiving devices.

(65) While the technique utilizing either linear or non-linear solution exemplified previously applies true-range multilateration with two (for two dimensional) or more (for three-dimensional), it is noted that: 1) angular/orientation invariant range dependence is particularly required for this and obtained by summing all orthogonal components of transmit to receive couplings, and that 2) two-dimensional and three-dimensional solutions can also be found by solving the multi-axial orthogonal components (referred herein as MSNA to denote multi-axial single-node) through partial ratio of partial sums for orientation or direction angles, and summed components for range, therefore providing up to three-dimensional position and orientation with as few as a single reference orthogonal node outside the building.

(66) The benefit of the multi-node true-range multilateration technique is the simplicity in algorithms and calibrations, and very low coupling between orientation errors and position errors, however, the MSNA technique has the benefit of single node implementation, which is also simpler for system deployment and lower cost at the expense of more complex algorithms. A physics-based limitation in MSNA is that it is not possible to exactly decouple the complex (and if strong) secondary fields from both the ground and building. As a result, there are three exact solution regimes within which the solutions for MSNA can be found to be most accurate: 1. Distance from outer wall to mobile device inside building is approximately less than one bulk-building skin depth (Δr{tilde over (<)}δ.sub.b), where δ.sub.b=√{square root over (2/ωμσ.sub.b)} is the building bulk skin depth with σ.sub.b is the bulk building conductivity. Secondary currents or fields from building can be neglected. 2. Distance from base station or outer reference node to mobile device is approximately greater than distance from base station or outer reference node to ground complex image (r{tilde over (>)}2h+δ.sub.g) where δ.sub.g=√{square root over (2/ωμσ.sub.g)} is the ground bulk skin depth with σ.sub.g is the bulk ground conductivity. Secondary currents or fields from ground can be neglected. 3. Distance from outer wall to mobile device is approximately less than one bulk-building skin depth and distance from base station or outer reference node to mobile device is approximately greater than distance from base station or outer reference node to ground complex image. In this case, both Δr{tilde over (<)}δ.sub.b and r{tilde over (>)}2h+δ.sub.g is satisfied, however, system will have limited range, as frequency of operation is not optimized for maximum possible range. If we assume that Δr≈r−d, where d is the distance inside the building, then solving the set of two requirements imply 2h+δ.sub.g{tilde over (<)}r{tilde over (<)}δ.sub.b+d, which due to earth ground conductivity in the 0.1-10 mS/m imply a very limited range of operation at best within the building. As a result, this regime is unlikely to be useful in most building or house environments.

(67) In most cases noted above (1-3), the approximated MSNA solutions can be found by solving the previously presented coupling equations, given again for the voltage coupling between a single transmit and single receive element:

(68) $\begin{array}{cc}{V}_T\approx C_g{e}^{-a\mathrm{\kappa \Delta }r}\left[3\left(\hat{m}.Math.\hat{\rho }\right)\left(\hat{n}.Math.\hat{\rho }\right)C_{\varphi }-\hat{m}.Math.\hat{n}\right]\times \left[\frac{1}{r^3}+\frac{1}{r_i^3}\right]& \left(12\right)\end{array}$
where the approach for solution involves obtaining a representation of all coupling from orthogonal transmit elements to receive elements. In total, 3×3=9 unique couplings can be found. To simplify the analysis for 2D (two-dimensional), the only need is to develop the 2×2=4 unique solutions that are sensitive to 2D position and orientations.

(69) Applicant will focus on the two-dimensional (2D), x, y, ϕ, 3 degrees of freedom (DOF) problem for simplicity and proof of concept, which requires two orthogonal dipoles at the transmitter (TX) and receiver (RX) (2×2=4 field couplings). To permit linear sensitivity to position and orientation, the field couplings are rewritten in transmit-centric and receive-centric forms (see Ref. [1-2]) as follows:
F.sub.m1=|V|.sub.(m1.fwdarw.n1).sup.2+|V|.sub.(m1.fwdarw.n2).sup.2
F.sub.m2=|V|.sub.(m2.fwdarw.n1).sup.2+|V|.sub.(m2.fwdarw.n2).sup.2
F.sub.n1=|V|.sub.(n2.fwdarw.m1).sup.2+|V|.sub.(n2.fwdarw.m2).sup.2
F.sub.n2=|V|.sub.(n2.fwdarw.m1).sup.2+|V|.sub.(n2.fwdarw.m2).sup.2  (13)
where F.sub.m1 and F.sub.m2 are the transmit-centric equations for the orthogonal transmitter TX 1 and 2 (m1,m2), which describe the total field power measured at all receivers RX(s) due to each TX, and F.sub.n1 and F.sub.n2 are the receive-centric equations for orthogonal receivers RX 1 and 2 (n1,n2), which describe the total power of all transmitters TX(s) measured at each receiver RX (see Ref. [1-2]).

(70) To obtain ranging that is orientation invariant, the contributions from both transmit centric equations are summed, as seen from Ref. [1]. For the present problem, this gives:

(71) $\begin{array}{cc}{F}_r=\left\{\begin{array}{c}{F}_{m1}\\ +F_{m2}\end{array}\right\}\approx 5{C_g^2\left[\frac{e^{-2a\mathrm{\kappa \Delta }r}}{r^3+r_i^3}\right]}^2=5{C_g^2\left[\frac{e^{-2a\kappa H\left(r-d\right)}}{r^3+r_i^3}\right]}^2& \left(14\right)\end{array}$
where H is the Heaviside unit step function. Measurements of F.sub.r obtained from summing all measured powers from the 2×2 orthogonal coupling is then described by a range dependence as noted, and can be inverted for effective range to device, r, when suitable assumptions of the various constants are known (C.sub.g, α, κ, δ.sub.g). This can be solved by numerical fitting, least square minimization or optimization, or similar for range to the device, r.

(72) The angular direction and orientation is substantially simplified, as the various constants due to range and building or ground image is removed by ratio functions needed for orientation or direction as seen from Ref. [1-2] as example.

(73) The angular direction to the receive device is given by (see for example Ref. [1]):

(74) $\begin{array}{cc}{F}_{\varphi m}=\frac{5-3\cos \left(2{\varphi }_{m1}\right)}{5+3\cos \left(2{\varphi }_{m1}\right)},{\varphi }_m={\varphi }_{m1}=\frac{1}{2}{\cos }^{-1}\left[\frac{5-5\times { }^S{F}_{\varphi m}}{3\times { }^S{F}_{\varphi m}+3}\right]& \left(15\right)\end{array}$

(75) The orientation of the receive device is given by (see for example Ref. [1]):

(76) $\begin{array}{cc}{F}_{\varphi n}=\frac{5-3\cos \left(2{\varphi }_{n1}\right)}{5+3\cos \left(2{\varphi }_{n1}\right)},{\varphi }_n={\varphi }_{n1}=\frac{1}{2}{\cos }^{-1}\left[\frac{5-5\times { }^S{F}_{\varphi n}}{3\times { }^S{F}_{\varphi n}+3}\right]& \left(16\right)\end{array}$

(77) It is noted that the 2D description above can be extended to 3D without loss of generality, and that further additional adjustments to this description can be made by simple manipulations of constants to provide suitable improvements in the fitting of this approximate simplified solution to the measurements obtained in the field. For the 2D example given above, 2D device position and 2D device relative orientation in Cartesian coordinate is obtained by simple transformations:
x=r cos ϕ.sub.n
y=r sin ϕ.sub.n
ϕ=ϕ.sub.m−ϕ.sub.n

(78) The methods and systems described in the present disclosure may be implemented in hardware, software, firmware or any combination thereof. Features described as blocks, modules or components may be implemented together (e.g., in a logic device such as an integrated logic device) or separately (e.g., as separate connected logic devices). The software portion of the methods of the present disclosure may comprise a computer-readable medium which comprises instructions that, when executed, perform, at least in part, the described methods. The computer readable medium may comprise, for example, a random access memory (RAM) and/or a read-only memory (ROM). The instructions may be executed by a processor (e.g., a digital signal processor (DSP), an application specific integrated circuit (ASIC), a field programmable logic array (FPGA), a graphic processing unit (GPU) or a general purpose GPU).

(79) A number of embodiments of the disclosure have been described. Nevertheless, it will be understood that various modifications may be made without departing from the spirit and scope of the present disclosure. Accordingly, other embodiments are within the scope of the following claims.

(80) The examples set forth above are provided to those of ordinary skill in the art as a complete disclosure and description of how to make and use the embodiments of the disclosure, and are not intended to limit the scope of what the inventor/inventors regard as their disclosure.

(81) Modifications of the above-described modes for carrying out the methods and systems herein disclosed that are obvious to persons of skill in the art are intended to be within the scope of the following claims. All patents and publications mentioned in the specification are indicative of the levels of skill of those skilled in the art to which the disclosure pertains. All references cited in this disclosure are incorporated by reference to the same extent as if each reference had been incorporated by reference in its entirety individually.

(82) It is to be understood that the disclosure is not limited to particular methods or systems, which can, of course, vary. It is also to be understood that the terminology used herein is for the purpose of describing particular embodiments only, and is not intended to be limiting. As used in this specification and the appended claims, the singular forms “a,” “an,” and “the” include plural referents unless the content clearly dictates otherwise. The term “plurality” includes two or more referents unless the content clearly dictates otherwise. Unless defined otherwise, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which the disclosure pertains.

(83) The references in the present application, shown in the reference list below, are incorporated herein by reference in their entirety.

REFERENCES

(84) [1] D. D. Arumugam. Through-the-wall magnetoquasistatic ranging. IEEE Antennas and Wireless Propagation Letters, 16:1439-1442, 2017. [2] D. Arumugam. Single-anchor 2d magnetoquasistatic position sensing for short to long ranges above ground. Antennas and Wireless Propagation Letters, IEEE, PP(99):1-1, 2015. [3] D. D. Arumugam. Decoupled range and orientation sensing in long-range magnetoquasistatic positioning. Antennas and Wireless Propagation Letters, IEEE, 14:654-657, 2015.