SYSTEM AND METHOD FOR CALIBRATING MAGNETIC SENSORS IN REAL AND FINITE TIME

Abstract

The present invention relates to a calibration method of magnetic sensors, for removing from the measurements the so-called bias and obtaining the actual measurement of the magnetic field. This method, in addition to measurements of such magnetic sensors, uses measurements of angular rotation sensors available in various types of commercial devices including smartphones. The present invention also relates to a corresponding system for determining the instantaneous real time orientation and/or position of a mobile device with respect to a magnetic field.

Claims

1. A method for real-time calibration of the magnetic field detected by at least one magnetic field sensor, to the magnetic field sensor being connected in an integral manner at least an angular rotation sensor, the method comprising the following steps performed by an electronic processor locally or remotely connected to said at least one magnetic field sensor and said at least one angular rotation sensor, the method comprising: A. acquiring, along time, a vector signal of the magnetic field by said at least a magnetic field sensor, said magnetic field vector signal being expressed in a sensor reference system, and comprising a magnetic bias; B. acquiring a rotation vector signal by said at least an angular rotation sensor, along time; C. sampling said vector signal of magnetic field and said rotation signal obtaining respectively a sampled magnetic field vector signal,
m(t)=[m.sub.x(t),m.sub.x(t),m.sub.x(t)] and a sampled rotation signal
(t)=[.sub.x(),.sub.y(t),.sub.z(t)] at time t; D. constructing the following orientation matrix: $S_{}\left(t\right)=\left[\begin{array}{ccc}0& -{}_z\left(t\right)& {}_y\left(t\right)\\ {}_z\left(t\right)& 0& -{}_x\left(t\right)\\ -{}_y\left(t\right)& {}_x\left(t\right)& 0\end{array}\right]$ wherein the following steps are performed by said electronic processor: E. numerically solving the following algebraic equation:
D{circumflex over (b)}(t)=Y where $D=\left[\begin{array}{c}{D}_1\\ \mathrm{.Math.}\\ D_n\end{array}\right].Math..Math.\mathrm{and}.Math..Math.Y=\left[\begin{array}{c}{Y}_1\\ \mathrm{.Math.}\\ Y_n\end{array}\right]$ in which {circumflex over (b)}(t) is an estimate of the bias b and D.sub.q is a function of S.sub.(t) and Y.sub.q is a function of m(t) and S.sub.(t), {circumflex over (b)}(t) being independent of the time derivative of the magnetic field components of the three-dimensional vector m(t), on each respective time interval q=1, . . . , n defined for [(q1)]T,qT], wherein n is a positive integer and T a predefined time interval in which measurements by sensors of steps A and B are to be performed, or the sampling time of step C; and F. calculating a calibrated magnetic field m.sub.c(t) as m.sub.c(t)=m(t){circumflex over (b)} expressed in said sensor reference system.

2. The method according to claim 1, wherein: $.Math..Math.T=T$ $n=\frac{T_w}{T}$ in which T is a specified period of time and T.sub.w is a pre-defined observation window, sufficiently wide as to ensure that the matrix D is full rank, and wherein:
D.sub.q=.sub.0.sup.T.sub.K()S.sub.(+(q1)T)dY.sub.q=.sub.0.sup.T.sub.K()S.sub.(+q1)T)m(+(q1)T)d.sub.0.sup.T{dot over ()}.sub.K()m(+(q1)T)d wherein .sub.K(t)C.sup.K is a function called modulating function defined in the space of K times differentiable functions, with K positive integer, and is defined on a finite time interval [0,T] and satisfies the following boundary conditions: ${{\frac{d^i.Math.{}_K\left(t\right)}{{\mathrm{dt}}^i}.Math.}_{t=0}=\frac{d^i.Math.{}_K\left(t\right)}{{\mathrm{dt}}^i}.Math.}_{t=T}=0,i=0,1,\mathrm{.Math.}.Math.,K-1$ and wherein the algebraic equation:
D{circumflex over (b)}=Y is solved numerically in the time window T.sub.w, and wherein the same algebraic equation is then solved for a subsequent time interval from T to T+T.sub.W and so on for successive windows of width T.sub.w.

3. The method according to claim 2, wherein: $\frac{d^i.Math.{}_K^s\left(t\right)}{{\mathrm{dt}}^i}=\{\begin{array}{cc}\underset{j=0}{\overset{K}{.Math.}}.Math.{\left(-1\right)}^j.Math.\left(\begin{array}{c}K\\ j\end{array}\right).Math.g_{\mathrm{ji}}\left(t-j.Math.\stackrel{\_}{T}\right),& i=0,\mathrm{.Math.}.Math.,K-1\\ \underset{j=0}{\overset{K}{.Math.}}.Math.{\left(-1\right)}^j.Math.\left(\begin{array}{c}K\\ j\end{array}\right).Math.\left(t-j.Math.\stackrel{\_}{T}\right),& i=K\end{array}.Math..Math.\mathrm{where}.Math..Math.\stackrel{\_}{T}=\frac{T}{K}.Math..Math.\mathrm{and}.Math..Math.g_{\mathrm{ij}}\left(t-j.Math.\stackrel{\_}{T}\right)=\{\begin{array}{cc}\frac{1}{\left(K-i-1\right)!}.Math.{\left(t-j.Math.\stackrel{\_}{T}\right)}^{K-i-1},& t\left[j.Math.\stackrel{\_}{T},T\right]\\ 0,i=K,& \mathrm{otherwise}\end{array}$

4. The method according to claim 1, wherein: $D_q={\left[\left(\frac{1}{s+}\left[K\left(t,t\right).Math.S_{}\left(t\right)\right]\right)\right]}_q$ $Y_q={\left[K\left(t,t\right).Math.m\left(t\right)-\left[V_{\frac{}{}.Math.K\left(t,\right)}.Math.m\right].Math.\left(t\right)+\left[V_K.Math.S_{}.Math.m\right].Math.\left(t\right)\right]}_q$ wherein Y.sub.q and D.sub.q axe discretizations on each interval q=1, . . . , n of the matrices shown in the argument $\frac{1}{s+}$ [.Math.] is a filter or the function K (t,t)S.sub.(t) in the argument, s is the Laplace variable, p is a non-null positive real number, is a time variable, V.sub.K is the Volterra operator, and $V_{\frac{}{}.Math.K\left(t,\right)}$ is the Volterra operator applied to the derivative of K(t, ), K(t,) being defined as:
K(t,)=e.sup.(t)(1e.sup.)

5. The method according to claim 4, wherein D.sub.q(t) is a so-called function persistently exciting in .sup.33, wherein there are two constants r>0, T>0 such that for every t0 the following occurs:
.sub.t-T.sup.tD.sub.q.sup.T()D.sub.q()drI>0 and the bias estimation is given by the solution of the following adaptive differential equation:
{circumflex over ({dot over (b)})}(t)=M.sub..sup.1(t)(sign(E(t))|{square root over (E(t)|)}+{dot over ()}.sub.(t){dot over (M)}.sub.(t){circumflex over (b)}(t)) wherein $.Math.M_f\left(t\right)={\left[{\left(\frac{1}{s+}\left[K\left(t,t\right).Math.S_{}\left(t\right)\right]\right)}^T.Math.\left(\frac{1}{s+}\left[K\left(t,t\right).Math.S_{}\left(t\right)\right]\right)\right]}^f$ ${}_f\left(t\right)={\left[{\left(\frac{1}{s+}\left[K\left(t,t\right).Math.S_{}\left(t\right)\right]\right)}^T.Math.\left(K\left(t,t\right).Math.m\left(t\right)-\left[V_{\frac{}{}.Math.K\left(t,\right)}.Math.m\right].Math.\left(t\right)+\left[V_K.Math.S_{}.Math.m\right].Math.\left(t\right)\right)\right]}^f$ $.Math.E\left(t\right)=M_f\left(t\right).Math.\hat{b}\left(t\right)-{}_f\left(t\right)$ wherein the apex indicates a filtering operation of the argument values in the square brackets of M.sub.(t) and .sub.(t).

6. The method according to claim 5, wherein the filtering operation is of the low-pass type.

7. A method for the real time determination of the instantaneous orientation of a mobile device with respect to an axis of a reference system of a magnetic field, in particular the geomagnetic field, and/or for the determination of the position of said device with respect to said reference system, the mobile device comprising or being constituted by at least a magnetic field sensor and at least an angular rotation sensor, in which the magnetic field is determined according to the method of claim 1, and wherein a further step G is performed, by calculating the orientation of said device with respect to said magnetic field on the basis of the calibrated vector m.sub.c(t), on the basis of said reference system, and of a sensor of inclination of said sensor with respect to the plane perpendicular to said axis as measured by at least a sensor field in the case where the inclination is null.

8. A computer program comprising code means configured to perform the steps A to F according to claim 1.

9. A computer program comprising code means configured to perform step G of claim 7.

10. A system for real time determination of the instantaneous orientation and/or the position of a mobile device with respect to a magnetic field defined in a reference system, comprising: at least a magnetic field sensor; at least an angular rotation sensor; at least a field sensor only in the case of the position determination; connected integrally to said device, and in local or remote communication with an electronic processor, wherein said computer program according to claim 9 is installed on said electronic processor.

11. The system according to claim 10, wherein said mobile device is a smartphone or tablet or PDA equipped internally with at least a magnetic field sensor and at least an angular rotation sensor, and said electronic processor is its CPU.

12. The system according to claim 10, wherein said magnetic field is the geomagnetic field.

Description

DETAILED DESCRIPTION OF EMBODIMENTS OF THE INVENTION

[0040] Further features and advantages of the invention will become apparent from the description of preferred embodiments, which are not exclusive or limiting, of a method for the real-time and finite-time calibration of magnetic sensors according to the invention, illustrated only by way of indication and not by way of limitation in the accompanying drawings, in which:

[0041] FIGS. 1A and 1B show the locus of the magnetic field measurements acquired respectively by using a uncalibrated magnetometer and a calibrated magnetometer.

[0042] FIGS. 2A, 2B, 3C illustrate the two reference systems. In particular, FIG. 2A shows a known reference system, external to the sensor, and known as the NED (North, East, Down) reference system. FIG. 2B shows a reference system fixed to the device which carries out the magnetic and the angular rotation measurements. FIG. 2C shows said overlapped reference systems highlighting the angles that can be determined between said systems.

[0043] FIG. 3 shows the axes of the reference system fixed to a device capable of measuring both the magnetic field vector and the orientation vector.

[0044] FIG. 4 illustrates a representative flow diagram of a method of estimation according to the invention;

[0045] FIGS. 5 and 6 show two flow diagrams relative to two embodiments of a method for calibration of magnetic sensors according to the invention; and

[0046] FIG. 7 shows an embodiment of the system according to the invention.

MATHEMATICAL DERIVATION

[0047] According to the invention, it is first noted that the previous EQ. 6 can be generically represented in the form:

((t))b=(m(t),(t),{dot over (m)}(t))EQ. 7

wherein m(t) are the measurements obtained by the magnetometer, (t) the measurements obtained by the tool measuring the angular rotations, {dot over (m)}(t) is the derivative of the measurement obtained by the magnetometer, ((t)) is a generic 33 matrix whose elements are a function of the angular rotations measurements and (m(t),(t),{dot over (m)}(t)) is a vector function 31 of the magnetic field measurements, the measurements of angular rotations and the derivative of the magnetic field measurements.

[0048] With particular reference to the above figures, the numeral 400 globally indicates a real-time and finite time calibration method of a magnetometer that, according to the present invention, has a particular characteristic in that it comprises: [0049] A first step (410) of processing during which said differential equation formulation of the problem (EQ. 6) ((t))b=(m(t), (t),{dot over (m)}(t)) is transformed into a new algebraic equation that cannot be solved in closed form ((t))b=(m(t), (t),{dot over (h)}(t)), wherein the original derivation operation applied to the measured data of the magnetic field {dot over (m)}(t) of that known differential equation is no longer present on the measured signal m(t) but it has been transferred to another term {dot over (h)}(t) of that new algebraic equation that cannot be solved in closed form. Such new algebraic equation cannot be solved in closed form since the function ((t)) is not invertible. [0050] A second step (420) of reformulation of said new algebraic equation that cannot be solved in closed form, obtained in said step of processing, in a new algebraic equation reformulated and solvable in finite-time ((t))b=(m(t), (t), {dot over (h)}(t)) and which allows to calculate the solution in terms of estimation of the optimal bias in finite time. Said reformulated algebraic equation can be solved in closed form since the function ((t)) is constructed so as to be invertible. [0051] A third step (430) for estimating the bias in finite time which provides: [0052] A first data acquisition operation (431) during which the measurements of the magnetic field vector and the angular rotation vector are acquired. [0053] A second operation of construction (432) of the n-dimensional matrices related to the angular rotation and magnetic field measurements obtained in said data acquisition operation. [0054] A third operation of numeric resolution (433) of said new algebraic equation reformulated and identified in said step of reformulation for the obtaining of the optimal numerical estimate of the calibration bias in finite time using said n-dimensional matrix of angular rotation measurements and said acquisitions of the earth's magnetic field vector.

[0055] In other words, said first processing step that implements a calibration method according to the present invention produces the effect of transforming the known differential equation Eq. 7 into a new algebraic equation which has the following form:

((t))b=(m(t),(t),{dot over (h)}(t))EQ. 8

where it is apparent that the derivative operation has been eliminated from the term m(t) and it has been introduced onto another generic term h(t) that does not depend on either the magnetic field measurements or the orientation measurements.

[0056] For example, in accordance with the present invention, two preferred embodiments are described which respectively obtain said result of eliminating the derivatives from the magnetic field measurement by applying two different procedures which make said function {dot over (h)}(t) dependent on results respectively obtained by applying the theory of Volterra operator and that of the modulating functions.

[0057] In a first preferred embodiment (500), said step of processing (410) provides to transform the differential equation described by EQ. 6 into a new algebraic equation using the Volterra operator. The prior art relating to the application of Volterra operator for the transformation of differential equations is described for example by T. Burton Volterra Integral and Differential Equations, Elsevier, 2005, ISBN 978-0-444-51786-9.

[0058] Given a function .sub.loc.sup.2 (.sub.0), where .sub.loc.sup.2 is the space of square integrable functions on every compact subset of .sub.0, and .sub.0 is the set of real numbers greater than or equal to zero, its image through the operator (integral linear) of Volterra V.sub.K induced by a Hilbert-Schmidt kernel function (the space of these functions is indicated with ) K(.Math.,.Math.): .fwdarw. (wherein is the set of real numbers) is indicated by [V.sub.K](.Math.), and it is defined by the product:

[V.sub.K](t).sub.0.sup.tK(t,)()d,t.sub.0EQ. 9

[0059] With reference to said first preferred embodiment, the chosen Hilbert-Schmidt Kernel function K(.Math.,.Math.) must respect the following properties for a given i2: [0060] K(.Math.,.Math.) and admits i-th derivatives with respect to the second argument along g the whole positive real axis;

[00002]${\frac{{}^j}{{}^j}.Math.K\left(t,\right).Math.}_{=0}=0,t{}_0,j\left\{0,\mathrm{.Math.}.Math.,i-1\right\}.$

[0061] Said first preferred embodiment (500) applies the theory of Volterra operator and the properties of said chosen kernel function to said known differential equation of formulation of the problem (EQ. 6) to transform it into a new algebraic equation that cannot be solved in closed form, wherein the derivative of the signal m(t) does not appear any longer.

[0062] Advantageously, said preferred embodiment uses the following Kernel function:

K(t,)=e.sup.(t-)(1e)EQ. 10

[0063] Said Kernel function eliminates said derivative of the signal m(t) and at the same time allows that some of the terms of said algebraic equation are easily calculated in a way known in itself, those terms being the analytical representation of a dynamic system of the first order (low pass filter)

[0064] In said first preferred embodiment (500), said step of processing (410) provides (510): [0065] a first operation (511) that provides to implement the Volterra operator (EQ. 9) and said function Kernel (EQ. 10) To said known differential equation (EQ. 6) OF The formulation of the problem, getting a new intermediate equation:

[V.sub.KS.sub.b](t)=[V.sub.KS.sub.m](t)+[V.sub.K{dot over (m)}](t)EQ. 11 [0066] a second operation of exploitation of the properties of said Volterra operator (512) through which said intermediate equation (EQ. 11) can be rewritten so as to obtain said new algebraic equation that cannot be solved in closed form, wherein the derivative of the term m(t) no longer appears, which was present in that original known equation defining the problem (EQ. 6), but said branching operation has been transferred onto the Kernel function:

[00003]$\begin{array}{cc}\left[V_K.Math.S_{}\right].Math.\left(t\right).Math.b=\left[V_K.Math.S_{}.Math.m\right].Math.\left(t\right)+K\left(t,t\right).Math.m\left(t\right)-\left[V_{\frac{}{}.Math.K\left(t,\right)}.Math.m\right].Math.\left(t\right)& \mathrm{EQ}..Math.12\end{array}$

[0067] Advantageously, said second operation (512) of exploitation of the properties of the Volterra operator allows to rewrite some of the terms of EQ. 12 by the following relationships holding for the above low-pass filter of the first order:

[00004]$\begin{array}{cc}\left[V_K.Math.S_{}\right].Math.\left(t\right)=\frac{1}{s+}\left[K\left(t,t\right).Math.S_{}\left(t\right)\right].Math.\left[V_K.Math.S_{}.Math.m\right].Math.\left(t\right)=\frac{1}{s+}\left[K\left(t,t\right).Math.S_{}\left(t\right).Math.m\left(t\right)\right].Math.\left[V_{\frac{}{}.Math.K\left(t,\right)}.Math.m\right].Math.\left(t\right)=\frac{1}{s+}\left[\frac{}{}.Math.K\left(t,t\right).Math.m\left(t\right)\right]& \mathrm{Eq}..Math.13\end{array}$

[0068] In said substitutions in (EQ. 13), the term s represents the Laplace variable, the notation P(s)[u(.Math.)] indicates the output of a system in the time domain P(s) having as input the signal u(.Math.). The execution of these replacements (EQ. 13) into said intermediate equation (EQ. 12) obtained in said first step of processing (511) provides, as a result, said new algebraic equation that cannot be solved in closed form:

H(t)b=(t)EQ. 14

wherein we have indicated with the term (t) the following quantity

[00005]$\left(t\right)=K\left(t,t\right).Math.m\left(t\right)-\left[V_{\frac{}{}.Math.K\left(t,\right)}.Math.m\right].Math.\left(t\right)+\left[V_K.Math.S_{}.Math.m\right].Math.\left(t\right)$$\mathrm{and}$$H\left(t\right)=\frac{1}{s+}\left[K\left(t,t\right).Math.S_{}\left(t\right)\right]$

which can be calculated in real time and can be assumed as known since it does not depend on the unknown value of the bias b. Said new algebraic equation cannot be solved in closed form since the matrix H(t) is not reversible.

[0069] Advantageously, in said first preferred embodiment (500), it is assumed that the orientation measurements obtained in said acquisition step (431) are in agreement with the hypothesis of persistent excitability described in the following. Based on said assumption, then the matrix H(t) of said new algebraic equation (EQ. 14) is persistently exciting in .sup.33 in the sense that there are two constants r>0, T>0 such that for every t0 the following occurs:

.sub.t-T.sup.tH.sup.T()H()drI>0EQ. 15

[0070] Note that as the said matrix H(t).sup.33 is by definition anti-symmetric and, by the theorem of Jacobi, is not full rank and therefore is not invertible. If we assume that the system is under the condition of persistent excitability, then a filtered version of the matrix H.sup.T()H() (obtained by passing through a filterfor example a low pass filterthe elements of the matrix) has full rank and is therefore invertible.

[0071] In said first preferred embodiment (500), in said step of reformulation (420) of said new algebraic equation that cannot be solved in closed form (EQ. 14), it is provided to re-write that new algebraic equation assuming that said assumption of persistent excitability is valid, by proceeding (520): [0072] To a first reformulation operation (521) which provides to multiply both members of said new algebraic equation (EQ. 14) by the matrix H.sup.T, getting H.sup.T (t)H(t)b=H.sup.T (t)(t); [0073] To a second reformulation of operation (522), which provides to replace in the result of said first operation the values (t)=H.sup.T(t)(t) and M(t)=H.sup.T(t)H(t), obtaining M(t)b=(t). If the value of the bias estimate b coincides with the exact unknown value of the bias, the previous can be written as M(t)b(t)=0. Note that the terms (t) and M(t) are calculated at each iteration. [0074] A third optional reformulation operation (523) which provides: [0075] A definition step which defines two filtered versions M.sub.(t) and .sub. (t) of M(t) and of (t):

[00006]$\{\begin{array}{c}{\stackrel{.}{M}}_f\left(t\right)=-.Math..Math.M_f\left(t\right)+M\left(t\right)\\ {\stackrel{.}{}}_f\left(t\right)-.Math..Math.{}_f\left(t\right)+\left(t\right)\end{array}$ [0076] wherein .sub.>0, M.sub.(0)=0.sub.3x3 and (0)=0.sub.3 with 0.sub.3x3 and 0.sub.3 respectively matrices of zeroes and vector of zeroes; and [0077] A step of replacement that replaces said filtered versions, M.sub.(t) and .sub.(t) in place of the nominal quantities M(t) and (t) in the intermediate result obtained in said second resolution operation, obtaining M.sub.(t)b.sub.(t)=0. [0078] To a fourth optional resolution operation (524). In said fourth resolution operation, it is observed that if {circumflex over (b)}(t) is an estimate of the sought bias b, then the equation M.sub.(t){circumflex over (b)}(t).sub.(t)=0 is not satisfied, and instead it has a residual E(t)=M.sub.(t){circumflex over (b)}(t).sub.(t) wherein (t)={dot over (M)}.sub.(t){circumflex over (b)}(t)+M.sub.(t){circumflex over ({dot over (b)})}(t){dot over ()}.sub.(t). Said fourth operation of resolution in finite time to obtain that |E(t)|.fwdarw.0, chooses in a particularly advantageous manner that:

{circumflex over ({dot over (b)})}(t)=M.sub..sup.1(t)(sign(E(t)){square root over (|E(t)|)}+{dot over ()}.sub.(t){dot over (M)}.sub.(t){circumflex over (b)}(t)) [0079] which was obtained using the known principle of the sliding mode (C. Edwards, SK Spurgeon, Sliding Mode Control, Taylor & Francis 1998). [0080] In this case, if said persistent excitability assumption is valid, then the matrix M.sub. has full rank and is reversible, therefore said equation is now solved and achieves the optimal estimation {circumflex over (b)}(t) in finite time of the sought bias vector b.

[0081] Other solutions could be envisaged to acquire a sufficient number of samples of the signals of interest (magnetometer and gyroscope) so as to set a least squares problem to obtain an estimate of b.

[0082] In said first preferred embodiment (500), said estimation step bias (430) provides that (530): [0083] the first data acquisition operation (431) uses a magnetometer and a triaxial gyroscope (without any need of an accelerometer as it is done in the prior art) for the magnetic field and orientation measurement (531). [0084] the second operation of construction of the magnetic field vector m(t)=[m.sub.x, m.sub.y, m.sub.z] and the matrix constructed from the orientation data (432) from respectively (532) the magnetic field measurements and the angular rotation measurements. If =[.sub.x, .sub.y, .sub.z] is a 3D acquisition of the orientation measurement, then the matrix is constructed as:

[00007]$S_{}\left(t\right)=\left[\begin{array}{ccc}0& -{}_z\left(t\right)& {}_y\left(t\right)\\ {}_z\left(t\right)& 0& -{}_x\left(t\right)\\ -{}_y\left(t\right)& {}_x\left(t\right)& 0\end{array}\right]$ [0085] A third operation of numeric resolution (433) of said new reformulated algebraic equation identified in said step of reformulation (520), that performs (533) and ensures the optimal estimation in finite time of the bias vector b.

[0086] This guarantee of finite-time resolution pertains to the theory described for example by Mohammad Pourmahmood Aghababa in Finite-time synchronization of two different chaotic systems with unknown parameters via sliding mode technique, Applied Mathematical Modelling, Volume 35, Issue 6, June 2011, Pages 3080-3091.

[0087] In a second advantageous embodiment (600), in accordance with the present invention (400), said step of processing (410) provides to transform the differential equation described by EQ. 6 into a new algebraic equation that cannot be solved in closed form using the modulating function theory (610). The prior art relevant to the use of the modulating functions theory is described for example by H. A. Preising and D. W. T. Rippin in Theory and application of the modulating function method-I. Review and theory of the method and theory of the spline-type modulating functions published in Computers & Chemical Engineering, Volume 17, Issue 1, January 1993, Pages 1-16. This known technique is used to estimate the unknown parameters of a differential equation.

[0088] A function .sub.K(t)C.sup.K (differentiable K times), defined over a finite time interval [0, T] and which meets the following boundary conditions:

[00008]$\begin{array}{cc}{{\frac{d^i.Math.{}_K\left(t\right)}{{\mathrm{dt}}^i}.Math.}_{t=0}=\frac{d^i.Math.{}_K\left(t\right)}{{\mathrm{dt}}^i}.Math.}_{t=T}=0,i=0,1,\mathrm{.Math.}.Math.,K-1& \mathrm{EQ}.Math..Math.16\end{array}$

is known as the modulating function.

[0089] A generic function (t).sup.1, where .sup.1 is the space of absolutely integrable functions defined on [0, T] is modulated by calculating the inner product of said generic function called modulating function .sub.K(t):

(t),.sub.K(t)=.sub.0.sup.T(t).sub.K(t)dtEQ. 17

[0090] The constraints at the boundaries, present inside EQ 16, substantially ensure that the boundary conditions of the function (t) are irrelevant after the modulation. Moreover, thanks to these constraints, one can transfer the derivation operation from the function (t) onto the modulating function .sub.K (t), in such a way that there is no more need to approximate the time derivative of a signal affected by noise:

[00009]$\begin{array}{cc}\frac{d^i.Math.f\left(t\right)}{{\mathrm{dt}}^i},{}_K\left(t\right)={\left(-1\right)}^i.Math.f\left(t\right),\frac{d^i.Math.{}_K}{{\mathrm{dt}}^i},i=0,\mathrm{.Math.}.Math.,K-1& \mathrm{EQ}..Math.18\end{array}$

[0091] Said second advantageous embodiment uses a spline function as a modulating function .sub.K(t). The n-th derivative of a spline function of order K and characteristic time

[00010]$\stackrel{\_}{T}=\frac{T}{K}$

is expressed by:

[00011]$\begin{array}{cc}\frac{d^i.Math.{}_K^s\left(t\right)}{{\mathrm{dt}}^i}=\{\begin{array}{cc}\underset{j=0}{\overset{K}{.Math.}}.Math.{\left(-1\right)}^j.Math.\left(\begin{array}{c}K\\ j\end{array}\right).Math.g_{\mathrm{ji}}\left(t-j.Math..Math.\stackrel{\_}{T}\right),& i=0,\mathrm{.Math.}.Math.,K-1\\ \underset{j=0}{\overset{K}{.Math.}}.Math.{\left(-1\right)}^j.Math.\left(\begin{array}{c}K\\ j\end{array}\right).Math.\left(t-j.Math..Math.\stackrel{\_}{T}\right),& i=K\end{array}& \mathrm{EQ}.Math..Math.19\end{array}$

[0092] wherein the term:

[00012]$g_{\mathrm{ij}}\left(t-j.Math..Math.\stackrel{\_}{T}\right)=\{\begin{array}{cc}\frac{1}{\left(K-i-1\right)!}.Math.{\left(t-j.Math..Math.\stackrel{\_}{T}\right)}^{K-i-1},& t\left[j.Math..Math.\stackrel{\_}{T},T\right]\\ 0,.Math.i=K,& \mathrm{otherwise}\end{array}$

[0093] We denote by T.sub.w the observation time of the measurements from the device (300) relevant to the magnetic field vector and to the vector of angular rotation and we denote with N the sub-intervals of length

[00013]$T=\frac{T_w}{N}$

seconds.

[0094] Advantageously, using said theory of modulating functions and said modulating function of the spline type, in said step of processing (410) of said second advantageous embodiment (600), in accordance with the present invention (400), one proceeds (610) to apply, for each time interval q=1, N defined for t[(q1)]T, qT], said inner product operator (EQ. 17) between that modulating function (EQ 19) and the differential equation of the calibration problem formulation (EQ. 6). Through this operation the original differential equation known (EQ. 6) is transformed into a new algebraic equation where the derivative of the term m(t) has been transferred to the term .sub.K(t) of the modulating function:

[.sub.0.sup.T.sub.K()S.sub.(+(q1)T)d]b=.sub.0.sup.T.sub.K()(S.sub.(+(q1)T)m(+(q1)T)+{dot over (m)}(+q1)T))dEQ 20

[0095] If we use the EQ. 18 in the previous EQ 20, we get:

.sub.0.sup.T.sub.K()S.sub.(+(q1)T)d]b=.sub.0.sup.T.sub.K()S.sub.(+(q1)T)m(+(q1)T)d+.sub.0.sup.T{dot over ()}.sub.K()m(+q1)T)dEQ 21

which can be written as a matrix equation where the derivative of the term m(t) has been transferred to the term .sub.K(t) of the modulating function:

G.sub.qb=.sub.qEQ. 22

wherein:

G.sub.q=.sub.0.sup.T.sub.K()S.sub.(+(q1)T)dEQ. 23

.sub.q=.sub.0.sup.T.sub.K()S.sub.(+(q1)T)m(+(q1)T)d.sub.0.sup.T{dot over ()}K()m(+(q1)T)d

[0096] Note that the EQ. 22 in accordance with said step of processing (410) of the proposed method (400) cannot be solved in closed form.

[0097] In said second advantageous embodiment (600), in accordance with the present invention (400), said step of reformulation of the algebraic equation (420) provides (620) to calculate an estimate {circumflex over (b)}(t) of bias b starting from said algebraic equation (EQ. 22) as obtained in said step of processing (610) using a least squares solution of the equation:

{circumflex over (b)}=G.sup.+ZEQ. 24

[0098] Said equation solvable by least squares is obtained by defining matrices

[00014]$G=\left[\begin{array}{c}{G}_1\\ \mathrm{.Math.}\\ G_N\end{array}\right].Math..Math.\mathrm{and}.Math..Math.Z=\left[\begin{array}{c}{}_1\\ \mathrm{.Math.}\\ {}_N\end{array}\right]$

which allow to rewrite EQ. 22 in the form Gb=Z where b is the bias of which the term {circumflex over (b)} of EQ. 24 is the optimum estimate.

[0099] In accordance with the present invention (400), the calibration method in question provides that the width of the observation window, T.sub.w, is chosen so as to ensure that the matrix G be of maximum rank.

[0100] Consequently, in said second advantageous embodiment (600), in accordance with the present invention (400), said bias estimation step (430) provides that (630): [0101] the first data acquisition operation (431) uses a triaxial magnetometer and a gyroscope for the magnetic field and orientation measurement (631) by acquiring said measures on a time window of T.sub.w seconds. [0102] the second operation of construction of the orientation array and the magnetic field array (432) provides for the relevant construction (632) starting from the measures relating to said time slot obtained in said first data acquisition operation by concatenating said magnetic field measurements inside a matrix

[00015]$Z=\left[\begin{array}{c}{}_1\\ \mathrm{.Math.}\\ {}_N\end{array}\right]$

and said orientation measures within a matrix

[00016]$G=\left[\begin{array}{c}{G}_1\\ \mathrm{.Math.}\\ G_N\end{array}\right]$ [0103] A third operation of numerical solution (433) of said algebraic equation reformulated as identified in said reformulation step (620) that performs (633) the optimal estimation in finite time of the bias vector b using a least squares technique on the data. [0104] A fourth step in which the same algebraic equation is then solved for a subsequent time interval from T to T+T.sub.w and so on for successive windows of magnitude T.sub.w.

A Calibration Method According to the Invention

[0105] Downstream of the above defrivation, the real-time calibration method of a magnetic field measurement according to the invention can be expressed in general terms as follows.

[0106] It uses at least one magnetic field sensor connected integrally to at least one angular rotation sensor and the corresponding measures to process a correction of the measurement bias. To this end, an electronic processor connected locally or remotely to said at least one magnetic field sensor and to said at least one angular rotation sensor performs the following steps: [0107] A. acquiring a vector signal of the magnetic field by said at least one magnetic field sensor along time, said magnetic field vector signal being expressed in a sensor reference system, and comprising a magnetic bias; [0108] B. acquiring a rotation vector signal by said at least an angular rotation sensor, along time; [0109] C. sampling said vector signal of magnetic field and said rotation signal obtaining respectively a sampled magnetic field vector signal, m(t)=[m.sub.x(t),m.sub.x(t),m.sub.x(t)] and a sampled rotation signal (t)=[.sub.x (t), .sub.y (t), .sub.z (t)] at the time t; [0110] D. constructing the following orientation matrix:

[00017]$S_{}\left(t\right)=\left[\begin{array}{ccc}0& -{}_z\left(t\right)& {}_y\left(t\right)\\ {}_z\left(t\right)& 0& -{}_x\left(t\right)\\ -{}_y\left(t\right)& {}_x\left(t\right)& 0\end{array}\right]$

[0111] Subsequently, on the basis of these data, the electronic processor estimates the bias according to the following steps: [0112] E. numerically solving the following algebraic equation:

D{circumflex over (b)}(t)=Y [0113] where

[00018]$D=\left[\begin{array}{c}{D}_1\\ \mathrm{.Math.}\\ D_n\end{array}\right].Math..Math.\mathrm{and}.Math..Math.Y=\left[\begin{array}{c}{Y}_1\\ \mathrm{.Math.}\\ Y_n\end{array}\right]$ [0114] in which {circumflex over (b)}(t) is an estimate of bias b and D.sub.q is determined on the basis of S.sub.(t) and Y.sub.q is determined on the basis of m(t) and S.sub.(t), with {circumflex over (b)}(t) being independent of the time derivative of the magnetic field components of the three-dimensional vector m(t), on each respective time interval q=1, . . . . , n defined for t[(q1)]T,qT], in which n is a positive integer and T a predefined time interval in which measurements by the sensors of steps A and B, or the sampling time of step C; and [0115] F. calculating a calibrated magnetic field m.sub.c(t) as m.sub.c(t)=m(t){circumflex over (b)} as expressed in said sensor reference system.

[0116] At this point, one can apply various special methods for the solution of the equation: D{circumflex over (b)}(t)=Y, since the general method already shown solves the problem of estimating a time derivative of the magnetic field.

[0117] For example by defining:

[00019]$.Math..Math.T=T$$n=\frac{T_w}{T}$

in which T is a specified period of time and T.sub.w is a pre-defined observation window sufficiently wide as to ensure that the matrix D is full rank, one can use the expressions above derived for the solution by the modulating functions, in which the algebraic equation: D{circumflex over (b)}(t)=Y is solved numerically in the time window T.sub.w, and wherein the same algebraic equation is then solved for a subsequent time interval from T to T+T.sub.w and so on for successive windows of width T.sub.w.

[0118] Another example of solution is given by the application of above illustrated kernel functions and Volterra operator with and subsequently to Eq. 14 (in this case D=M.sub.(t), Y=.sub.(t) and n=1.

[0119] More specific solutions are possible once the basic concept of the general formulation above is known.

Method for the Determination of Orientation and/or Position of a Sensor According to the Invention.

[0120] From erstwhile discussion, it is therefore clear that the correct magnetic field is a vector in the sensor reference system. For example, if the sensor is contained in a smartphone, the measured magnetic field will be expressed in the coordinates of the smartphone reference system. If one wants to then calculate the orientation of the smartphone (sensor) with respect to the magnetic field reference system, one will need to know the inclination of the smartphone with respect to an axis of the reference system. This, in some cases, may involve the use of a field sensor (e.g. gravitational field and therefore the use of an accelerometer), which measures the inclination of the smartphone with respect to said axis (in cases where it is not null by definition). This field sensor is also important for the determination of the position of the sensor in the magnetic field reference system.

System for the Determination of Orientation and/or Position of a Sensor According to the Invention

[0121] In FIG. 7 a diagram of a particular embodiment not limiting the invention is provided. A smartphone 300 indicates on the screen (optional) its orientation with respect to the local magnetic field. To do this, the magnetic sensor 330 provides the data to the CPU 320, which calibrates them by the data of the gyroscope 340, according to the method illustrated above. It is understood that this can also be done by a cloud server in communication with the smartphone, while the two sensors 330, 340 must be integral to the smartphone.

[0122] According to the method just described, an accelerometer or sensor field (eg. gravitational), not shown in the figure, can be included as well.

[0123] The CPU has code means capable of running a computer program that performs the steps of the methods above described starting from the sensor data.

[0124] The smartphone can be in general a mobile or a tablet device, a PDA or another.

Innovations and Advantages

[0125] In extreme synthesis, compared to conventional methods described, the method subject-matter of the present invention presents the following differences: [0126] It allows the determination of the optimal estimate of bias without the need for the magnetometer to undergo major changes in its orientation; [0127] It guarantees the convergence of the calculation method, and then determining the bias value; [0128] It guarantees the finite time of execution of the method of calculation.

[0129] From these differences the following advantages ensue: [0130] guarantee of obtaining of a resulting value for the estimation of the sought bias; [0131] guarantee on the execution time of the calculation which is performed in a timeframe that can be pre-determined and set in advance; [0132] non-dependence on orientations to which the magnetometer must be subjected, in order for the estimate of bias to be obtainable.

[0133] The invention thus conceived is susceptible to numerous modifications and variants, all falling within the scope of the appended claims.

[0134] Moreover, all the details may be replaced with other technically equivalent elements.

[0135] In practice, the means and materials employed, as well as the configurations, shapes and the dimensions, may be varied depending on the contingent requirements and state of the art.

[0136] Where the constructional characteristics and techniques mentioned in the following claims are followed by reference signs or numbers, such signs and reference numbers have been applied with the sole purpose of increasing the intelligibility of the claims and, consequently, they do not constitute in any way limitation of the interpretation of each element identified, purely by way of example, by such signs and reference numbers.