Method for estimating the direction of motion of an individual

Abstract

The present invention relates to a method for estimating the value of the angle of the direction of motion or of the walk of a person who carries a device capable of measuring the basic quantities related to his/her own inertia such as a common smartphone. The device is capable of measuring for example acceleration and angular rotation and to determine its orientation relative to Magnetic North. The estimation method in question is particularly suitable for estimating the direction angle of the individual's motion regardless of how the individual carries the device able to measure the main inertial quantities and consequently this method solves the problem of the determination of the relative attitude between the reference system of the carried device and the user reference system that moves within a generic absolute reference system.

Claims

1. Method for the estimation of an angle θ representing axis and direction of motion of an individual in a geomagnetic absolute Cartesian reference system North-East-Down (X.sub.aY.sub.aZ.sub.a), to whom a Cartesian reference system (X.sub.iY.sub.iZ.sub.i) of an individual carrying a device equipped with one or more inertia sensors is associated, the device being associated with a device Cartesian reference system (X.sub.dY.sub.dZ.sub.d) and comprising at least one accelerometer and at least one magnetometer, the method being characterized in that it comprises the execution of the following steps: A. splitting a straight angle into a set AS of n discrete angles θ.sub.j with j=1, . . . , n, where n is a positive integer; B. measuring, at a given time instant t in a time window, in the device Cartesian reference system (X.sub.dY.sub.dZ.sub.d), the device acceleration vector aced via said at least one accelerometer, and the local magnetic field vector m=(m.sub.x, m.sub.y, m.sub.z) by said magnetometer; C. determining an attitude A.sub.d of said device based at least one the local magnetic vector m; D. calculating, within the absolute Cartesian reference system (X.sub.aY.sub.aZ.sub.a), the linear acceleration vector acc.sub.1=(acc.sub.lx, acc.sub.ly, acc.sub.lz) to which said device is subjected during the motion of said individual, via the following substeps: D1. rotating the device acceleration vector aced onto the horizontal plane X.sub.aY.sub.a of said absolute reference system using the attitude A.sub.d, obtaining a new vector of acceleration acc.sub.h; D2. subtracting a gravitional component g=[0, 0, g.sub.za≅9.81 m/s.sup.2, from the new vector of acceleration acc.sub.h of D1 step, obtaining a linear acceleration vector acc.sub.l; E. Building a linear accelerations matrix AQ as follows: E1. calculating, for each angle θ.sub.jεAS, the following function of the linear acceleration vector acc.sub.l:
acc.sub.θ.sub.j=acc.sub.lx*cos(θ.sub.j)+acc.sub.ly*sin(θ.sub.j) E2. appending a new column c containing at least said values acc.sub.74 .sub.j, j=1, . . . , n, and the component acc.sub.lz of said vector acc.sub.l within a matrix AQ:
c=[acc.sub.θ.sub.1, . . . ,acc.sub.θ.sub.n,acc.sub.lz].sup.T wherein steps B to E are repeated for successive time instants in said time window, and wherein the following further steps are carried out: F. Calculating the variances of the {σ.sup.2(acc.sub.θ.sub.1), . . . ,σ.sup.2(acc.sub.1z)} of the quantities stored in the rows of said accelerations matrix AQ; G. estimating a discrete angle θ.sub.dir of the direction axis of motion of said individual as: ${\theta }_{\mathrm{dir}}=\arg \underset{{\theta }_j\in \mathrm{AS}}{\max }{\sigma }^2\left({\mathrm{acc}}_{{\theta }_j}\right)$ H. determining said angle θ representing the axis and the direction of motion of an individual, according to the following sub-steps: H1. Calculating acc.sub.74 .sub.dir a phase shift obtained by calculating a correlation index according to a deviation between reference values formed respectively by a first sequence and a second sequence, wherein both the first and second sequences are in said time window and the first sequence is obtained as
acc.sub.θ.sub.dir=acc.sub.lx*cos(θ.sub.dir)+acc.sub.ly*sin(θ.sub.dir) and the second sequence is given by
acc.sub.lz; wherein in said H1 step the phase shift is obtained by calculating a correlation index according to a deviation between the reference values formed by the row acciz in the matrix AQ and the second sequence formed by the values acc.sub.74 .sub.dir obtained by extracting the row associated with θ.sub.dir in the matrix AQ; H1A. starting from the zero value of said deviation, identifying: a first maximum value of the correlation index corresponding to a first deviation value, from among negative values of said deviation, and a second maximum value of the correlation index corresponding to a second deviation value, from among positive values of said deviation; and H1B. determining the angle θ as follows: If the absolute value of said first deviation value is less than or equal to the absolute value of said second deviation value, then θ=θ.sub.dir; If the absolute value of said first deviation value is greater than the absolute value of said second deviation value, then θ=θ.sub.dir+π.

2. Method according to claim 1, wherein the attitude A.sub.d of step C is represented by a quaternion q.sub.d and in that the D1 step comprises the following substeps: calculating a vector axis by carrying out a cross product of the gravitational vector g=[0,0,g.sub.Za] as expressed in said absolute Cartesian coordinate system (X.sub.aY.sub.aZ.sub.a) and said quaternion q.sub.d, wherein a new gravitational vector g.sub.s is obtained by transforming the gravitational vector g from the absolute reference system (X.sub.aY.sub.aZ.sub.a) to the device reference system (X.sub.dY.sub.dZ.sub.d) using the following equations:
g.sub.s=q.sub.d.Math.g.Math.q.sub.d.sup.−1
axis=g×g.sub.s where the term q.sub.d.sup.−1 is the conjugate of the quaternion q.sub.d, the operator represents the product between quaternions, and in which the gravitational vector g has been converted from a three-dimensional vector to a four-dimensional vector by inserting a zero component; calculating an angle ϕ between the two vectors g and g.sub.s; building a further quaternion q.sub.n from the values of axis and ϕ; rotating the vector acc.sub.d onto the horizontal plane X.sub.aY.sub.a of said absolute reference system by said further quaternion q.sub.n: ${\mathrm{acc}}_h=\left[\begin{array}{c}{\mathrm{acc}}_{h}x\\ {\mathrm{acc}}_{h}y\\ {\mathrm{acc}}_{h}z\end{array}\right]=q_n^{-1}.Math.{\mathrm{acc}}_d.Math.q_n$

3. Method according to claim 1, wherein before step E2, the following sub-steps are executed: determining a number of columns of matrix AQ or an extent of said time window; and if the number of columns of matrix AQ or the extent of said time window is greater than a predetermined respective value, deleting the oldest column of said matrix AQ.

4. Method according to claim 1, wherein the calculation of the F step is carried out for each column appended in the E2 step.

5. Method according to claim 1, wherein a correlation index is calculated using the Pearson correlation index.

6. Method according to claim 1, wherein in step B a digital filtering is applied to the signals measured by said at least one accelerometer.

7. System for the estimation of an angle θ representing the axis and the direction of motion of an individual in an absolute Cartesian geomagnetic reference system North-East-Down (X.sub.aY.sub.aZ.sub.a), to whom an individual's Cartesian reference system (X.sub.i Y.sub.i Z.sub.i) is associated, the system comprising: a device with one or more inertia sensors, the device being associated with a device Cartesian reference system (X.sub.dY.sub.dZ.sub.d) and comprising at least one accelerometer and at least one magnetometer, an electronic processor in data communication with said device; a memory storing instructions configured to, when executed, cause said electronic processor to perform steps A-H of the method according to claim 1.

8. System according to claim 7, wherein said electronic processor is integrated in said device.

9. System according to claim 7, wherein said device is a smart phone or tablet or PDA equipped internally of at least one magnetic field sensor and at least an accelerometer, and said electronic processor is its CPU.

Description

DETAILED DESCRIPTION OF EMBODIMENTS OF THE INVENTION

(1) Further characteristics and advantages of the invention will become more clear from the description of preferred, but not exclusive embodiments of a method for estimating the direction axis and direction of the walk of an individual according to the invention, illustrated only by way of illustration but not by way of limitation in the accompanying drawings, in which:

(2) FIG. 1 shows the reference systems used to facilitate the explanation of the method. FIGS. 1a and 1b show an absolute reference system known as NED (North, East, Down). FIG. 1c shows the reference system fixed to the individual wherein the X axis is aligned with the direction of the walk. FIG. 1d shows the reference frame of the device designed to measure the inertia.

(3) FIG. 2 identifies in detail the yaw angle ψ which substantially expresses the orientation of the device about the gravitational axis (coincident with the Z-axis), evaluated with respect to the X.sub.a axis of the absolute reference system, which by definition is aligned to North.

(4) FIG. 3 shows a total of the three reference systems of the previous figures. In FIG. 3a, the reference system of the individual who carries a device for measuring the inertia and the reference system of the same device are shown; in this case the device is a smartphone. FIG. 3b shows the absolute reference system.

(5) FIG. 4 shows a detail on the horizontal plane perpendicular to the Z axis of the three reference systems by identifying the angle of yaw ψ and the angle of the direction of motion θ estimated by this method.

(6) FIG. 5 shows one of the 180 degree angles

(7) $\left[-\frac{\pi }{2};\frac{\pi }{2}\right]$
extractable from a 360 degree angle, the straight angle having been divided into n=6 sectors of equal amplitude equal to

(8) $\frac{\pi }{6}$
radians.

(9) FIG. 6 shows the two absolute and inertia measuring device reference systems with evidence of the representation of the gravity vector g=(g.sub.x, g.sub.y, g.sub.z) as measured by an accelerometer located inside said device and with the indication of the two angles of roll (r) and pitch (p) representative of the inclination with respect to the axes X and Y.

(10) FIG. 7 shows the 180 degree angle already shown in FIG. 5 with evidence and indication of the directions.

(11) FIG. 8 shows a flow chart of the invention method for estimating the direction of the walk.

(12) FIG. 9a shows a diagram of the rough accelerations measured on the three axes of an accelerometer during a walk, as expressed in the reference system of the orientation measuring device. FIG. 9b shows the linear accelerations, expressed in the absolute reference system, as obtained starting from the raw acceleration and applying a low pass filter, rotating the orientation and eliminating the gravitational component.

(13) FIGS. 10a and 11a show a detail of the signals contained within a time window during the execution of the method. FIGS. 10b and 11b show the cross-correlation graph relating to said signals; and

(14) FIG. 12 shows a snapshot of the error between the real axis of motion and the estimation ones obtained by projecting on the horizontal plane using the attitude A.sub.d or the quaternion q (hereafter described), respectively.

(15) The present invention describes an analytical method capable of real-time estimating the axis and direction of the walk of an individual. In particular, this is done by determining the angle θ∈[0; 2π) shown in FIG. 4. The use of said angle, together with an estimation of the attitude of the inertia measurement device, obtainable by using any of the conventional methods such as that described by Magdwick et al., allows to reconstruct in an univocal and absolute way the axis and direction of motion of the individual's movement on a horizontal plane regardless of how the individual carries the inertia measuring device, i.e. irrespective of the relative orientation between the individual's reference system and the device reference system.

Definitions

(16) It is here defined as “attitude” of a device a minimal set of variables that allow to reconstruct, univocally, the orientation of the reference system integral to the device with respect to a fixed reference system.

(17) Such minimal set can be constituted, according to the known art, in an equivalent manner by: 1) a set of three angles (for example roll, pitch, yaw), 2) a 3×3 rotation matrix, 3) a quaternion (see eg. (F. L. Markley, J. L. Crassidis (2014), Fundamentals of Spacecraft Attitude Determination and Control. Springer).

(18) Using the measurement of the magnetic vector m=(m.sub.x, m.sub.y, m.sub.z) with respect to the North (in the X.sub.aY.sub.aZ.sub.a NED system illustrated in FIG. 3b) and the measurements of roll and pitch angles, it is possible to calculate the angle of yaw ψ, using known methods (e.g. Eq. 3 in the following).

(19) It is here defined as “linear acceleration” of a material point a three-dimensional vector acting on the plane X.sub.aY.sub.a of the inertial reference system (NED) that represents the component of the individual's motion acceleration perpendicular to gravity.

(20) In the following, as absolute reference system, reference will always be made to the NED system (FIG. 3b).

(21) Method for the Calculation of the Travel Direction

(22) Discretization of 180 Degree Angle

(23) The present invention determines the desired angle θ by means of two basic steps: a first computational step determines an estimate of the direction of motion, which provides an angle θ.sub.dir representative of the angle between the reference axis of the inertial system of the device and the axis along which the individual motion takes place. This angle between the axis of motion and the reference axis may coincide or be opposite to the angle of the axis of motion, depending on the direction of motion. A second computational step, in which such ambiguity on direction is dissolved and, starting from θ.sub.dir, the sought angle θ is obtained.

(24) With reference also to FIG. 8, in the first computational step, one chooses: a desired estimation accuracy Δθ; and a generic 180 degree angle (for example

(25) $\left[-\frac{\pi }{2};\frac{\pi }{2}\right]$
extracted from the 360 degree angle. It is noted here that from a 360 degree angle endless 180 degree angles can be extracted, depending on the origin of the angle and the direction of the angle. The arbitrariness of the choice does not diminish the generality of what it is said below.

(26) Subsequently, the 180 degree angle is discretized according to the desired accuracy, i.e. subdivided into equal parts, each having amplitude Δθ. The discretization is performed to reduce the computational cost and, with increasing desired accuracy, i.e. with decreasing Δθ, the computational load associated with calculation of the motion direction angle increases.

(27) In a preferred embodiment, the first step 410 of initialization provides to divide the 180 degree angle

(28) $\left[-\frac{\pi }{2};\frac{\pi }{2}\right)$
into n sectors of equal size. The chosen number of sectors n determines the level of accuracy obtained in the final estimate of the direction of motion (specifically the chosen accuracy will be of

(29) $\Delta \theta =\frac{\pi }{n}).$

(30) For example, if one chooses n=6 (FIG. 5), one can extract from the 180 degree angle

(31) $\left[-\frac{\pi }{2};\frac{\pi }{2}\right)$
n allowable values for the angle θ.sub.dir, which correspond to discrete angles set

(32) $AS=\left\{-\frac{\pi }{2};-\frac{\pi }{3};-\frac{\pi }{6};0;\frac{\pi }{6};\frac{\pi }{3}\right\}.$
It is sufficient to include only one of the endpoints of the 180 degree angle as the other extreme is on the same axis apart from the sign that precisely indicates the direction. The estimation of the direction of motion as performed with n=6, provides for selecting the estimate of the angle associated to said axis among those included in said discrete set of angles apart from the ambiguity that is resolved at a later step of the method according to the invention. This estimate is therefore affected by an uncertainty in the determination of the said direction equal to

(33) $\frac{\pi }{n}=\frac{\pi }{6}=30^{\circ },$
i.e. it is not allowed to estimate different directions if these directions are smaller changes

(34) 0$\frac{\pi }{6}$
radians.
Basic Steps of the Estimation Algorithm

(35) With particular reference to the above figures, the numeral 400 globally indicates a method for estimating the direction axis and direction of the walk of an individual who carries a device adapted to measure of inertia and that, according to the present invention, has a particular peculiarity in that it comprises: A first step 410 of initialization, during which the representative interval relative to one of the straight angles, for example

(36) $\left[-\frac{\pi }{2};\frac{\pi }{2}\right),$
as extractable by a round angle, is divided into n discrete angles, with n positive integer, which lay down the level of accuracy

(37) $\Delta \theta =\frac{\pi }{n}$
by which the final estimate of the angle of the direction of motion, according to this method, is determined. A second step 420 of measurement of the acceleration and the other inertial quantities useful to the attitude determination of said carried device adapted to measure the inertia. A third step 430 of calculation of the linear acceleration to which said transported device is subjected during the motion of said individual. A fourth step 440 of updating a floating data window containing said linear acceleration measurements rotated through said discrete angles, therefore calculating n rotations. A fifth step 450 of estimation of the discrete angle of the individual's axis of motion θ.sub.dir. A sixth step 460 of ambiguity resolution of the direction of said given axis as determined in said fifth estimation step so as to obtain the sought value θ starting from θ.sub.dir.
Measurements Made by the Measuring Device

(38) The measuring device, or any device equipped with inertial sensors, such as accelerometer, magnetometer (sensor adapted to measuring the attitude of said device with respect to the North) and, possibly, gyroscope, for example a smartphone, is used in the present invention for accumulating inertial data in a chosen time window.

(39) Preferably, said inertial sensors are triaxial sensors.

(40) The method according to the invention, in addition to the discretization of a straight angle extracted from a round angle, above, also makes a time discretization in time windows of predetermined length, managed by a floating window mechanism and in which the values acquired from the above sensors are stored and one calculates the angle θ representing the direction axis and direction of motion. In this way, one has a real-time estimate of that angle while the individual moves.

(41) Preferably, in said step of initialization, the number of measurements is fixed or the number of seconds k relative to the measurements from said inertia measuring device that are accumulated within a time window AQ and used to provide an estimate of said angle θ representing the direction of motion of said individual.

(42) In each time window AQ, one acquires 420 the vector acceleration acc.sub.d acting on the device and expressed in the reference system fixed with the same device, and other measurable quantities from said inertial sensors of 1s the device adapted to measure the inertia.

(43) Preferably, said second step 420 may comprise the application of digital filtering techniques (low pass, band pass, etc.) to the measured acceleration signals so as to reduce the effect of any noise on the measurement of the process of θ angle estimation.

(44) Computation of the Device Attitude

(45) To calculate the attitude, any of equivalents representations based on rotation matrices, Euler angles or unitary quaternions may be used, arriving at the same end result. For example, in the case in which the device is equipped only with accelerometer and magnetometer, using the acceleration vector acquired as above and a conventional method such as that described by Sergiusz Łuczak et al. in “Sensing Tilt With MEMS Accelerometers” IEEE Sensors Journal, Volume 6, Issue 6, Pages 1669-1675, December 2006, ISSN 1530-437X, one can derive the two angles (r, p) of roll and pitch shown in FIG. 6.

(46) Using the measurement of the geomagnetic vector m=(m.sub.x, m.sub.y, m.sub.z) with respect to the North (NED absolute frame of reference) and the roll and pitch angles, it is possible to calculate the yaw angle ψ shown in FIG. 2, by the following relation given in the just mentioned publication of Sergiusz uczak et al.:

(47) $\begin{array}{cc}\tan \left(\psi \right)=\frac{m_x\cos \left(r\right)+m_y\sin \left(r\right)\sin \left(p\right)+m_z\sin \left(r\right)\cos \left(p\right)}{m_z\sin \left(p\right)-m_y\cos \left(p\right)}& \mathrm{EQ}.3\end{array}$

(48) The roll and pitch angles can be calculated through traditional methods that rely only on the accelerometric measurements or also on the measurements obtained by a triaxial gyroscope, which may be optionally provided in the system of the invention. In the latter case, all available information from the accelerometer, magnetometric (for yaw) and gyroscopic data, as well as other traditional methods, such as the method described by said SOH Madgwick et al., are used.

(49) Thanks to these traditional methods, according to the invention, the attitude A.sub.d of said device (representable as shown in F. L. Markley, J. L. Crassidis (2014). Fundamentals of Spacecraft Attitude Determination and Control. Springer) with respect to said absolute reference system equivalently expressed by the Euler angles (for instance roll, pitch, yaw) or by a rotation matrix R∈.sup.3×3 or by a unitary quaternion q.sub.d={q.sub.0, q.sub.1, q.sub.2, q.sub.3} is obtained.

(50) Linear Acceleration Calculation

(51) In said preferred embodiment, said third step 430 of linear acceleration calculation includes: A first step, which uses the attitude information A.sub.d to rotate (see for example F. L. Markley, J. L. Crassidis (2014). Fundamentals of Spacecraft Attitude Determination and Control. Springer) the three-dimensional vector acc.sub.d acceleration measured in the horizontal plane X.sub.aY.sub.a of said absolute reference system, obtaining a new vector of acceleration acc.sub.h. A second step of elimination of the gravitational component from the vector acc.sub.h, to obtain the linear acceleration acc.sub.l:

(52) $\begin{array}{cc}ac{c}_l=\left[\begin{array}{c}ac{c}_{l}x\\ ac{c}_{l}y\\ ac{c}_{l^Z}\end{array}\right]=a{\mathrm{cc}}_h-gg=\left[0,0,g_{\mathrm{Za}}\right]& \mathrm{EQ}.4\end{array}$

(53) The value acc.sub.l represents the linear acceleration measurement obtained from the original acceleration as measured by said acceleration sensor of the device, rotated onto the horizontal plane of said absolute reference system to which the gravitational component has been subtracted. The gravitational component measured in the absolute reference system g=[0,0,g.sub.Za], if expressed in SI units, is equal to about g.sub.Za=9.81 m/s.sub.2.

(54) Preferably, in said third step 430, said first step provides for the use of roll and pitch angles for rotating the acceleration vector acc.sub.d onto the horizontal plane of said absolute reference system:

(55) $\begin{array}{cc}ac{c}_h=\left[\begin{array}{c}ac{c}_{h}x\\ ac{c}_{h}y\\ ac{c}_{h}z\end{array}\right]=\hspace{1em}=\left[\begin{array}{ccc}\cos \left(p\right)& 0& \sin \left(p\right)\\ 0& 1& 0\\ -\sin \left(p\right)& 0& \cos \left(p\right)\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& \cos \left(r\right)& -s\mathrm{in}\left(r\right)\\ 0& \sin \left(r\right)& \cos \left(r\right)\end{array}\right]\mathrm{ac}{c}_d& \mathrm{EQ}5\end{array}$
in which acc.sub.hx, acc.sub.hx, acc.sub.hx are the components of acc.sub.h in the reference system X.sub.aY.sub.aZ.sub.a.

(56) In the following, to avoid ambiguity problems related to the use of the Euler angles, (ambiguity known as “gimbal-lock”), in said third step 430, said first step provides for the use the said formalism of quaternions, thus defining quaternion q.sub.d that defines the attitude of said device.

(57) In this formalism, one proceeds as follows: Calculation of the axis of rotation vector axis obtained from the cross product of the gravitational vector g=[0,0,g.sub.Za] expressed in said absolute reference system and a new gravitational vector g.sub.s obtained by transforming, by the application of a quaternion q.sub.d that defines the structure A.sub.d, the gravitational vector g from the device reference system to the absolute reference system:
g.sub.s=q.sub.d.Math.g.Math.q.sub.d.sup.−1
axis=g×g.sub.s The term q.sub.d.sup.−1 is the conjugate of the quaternion q.sub.d, the operator represents the product between quaternions wherein the three-dimensional vector to be rotated has been increased by a dimension by inserting a 0 component. calculation of angle ϕ between the two vectors g and g.sub.s (FIG. 6), which integrates the information of roll and pitch, calculable for example, without limitation, as follows:
ϕ=−a tan 2(∥g×g.sub.s∥),<g,g.sub.s>) wherein the operator ∥.Math.∥ is the norm of a vector, the operator <.Math.,.Math.> is the scalar product operation between two vectors and the function a tan 2(y,x) represents the function

(58) $\mathrm{atan}\left(\frac{y}{x}\right)$
purged of the ambiguity of 180 degrees thanks to the analysis of the signs of the numerator y and denominator x, i.e.:

(59) $\mathrm{atan}2\left(y,x\right)=\{\begin{array}{cc}\mathrm{atan}\left(\frac{y}{x}\right)& \mathrm{if}x>0\\ \mathrm{atan}\left(\frac{y}{x}\right)+\pi & \mathrm{if}x<0\mathrm{and}\gamma \ge 0\\ \mathrm{atan}\left(\frac{y}{x}\right)-\pi & \mathrm{if}x<0\mathrm{and}y<0\\ \frac{\pi }{2}& \mathrm{if}x=0\mathrm{and}y>0\\ -\frac{\pi }{2}& \mathrm{if}x=0\mathrm{and}y<0\\ \mathrm{undefined}& \mathrm{if}x=0\mathrm{and}y=0\end{array}.$ Construction of a new quaternion q.sub.n from the values of axis and ϕ calculated in the previous steps, according to known mathematical methods (J. B. Kuipers, Quaternions & Rotation Sequences: A Primer with Applications to Orbits, Aerospace and Virtual Reality, Princeton University Press, 1999). Rotation of the original vector acc.sub.d onto the horizontal plane X.sub.aY.sub.a of said absolute reference system by said new quaternion q.sub.n:
acc.sub.h=q.sub.n.sup.−1.Math.acc.sub.d.Math.q.sub.n  EQ. 6
Update of the Floating Window

(60) In said preferred embodiment, said fourth step 440 to update the floating window includes a matrix representation of the floating window, in 1s which each column of the matrix represents a set of linear acceleration values rotated according to the set angles AS, To a time instant set, to be stored within the floating window itself. In particular, for the floating window management simplicity, the data is stored for time instants in the flow increasing from left to right the columns of the matrix AQ adopted to store the data of the floating window. It provides for the construction and update of floating window: A first calculation operation, for each angle θ.sub.j∈AS, with j=1, . . . , n, therefore belonging to the n angles of said set AS of discrete angles as determined in said first initialization step, of the following function of the vector acc.sub.l:
acc.sub.θ.sub.j=acc.sub.lx*cos(θ.sub.j)+acc.sub.ly*sin(θ.sub.j)  EQ. 7 A second appending operation of a (new) column c containing at least said values acc.sub.θ.sub.j, j=1, . . . , n, and the component acc.sub.lz of said vector acc.sub.l within the time window AQ, initialized as AQ=[ ](empty array).
c=[acc.sub.θ.sub.1, . . . ,acc.sub.θ.sub.n,acc.sub.lz].sup.T
AQ=[AQ,c]  EQ. 8 A third calculation operation, in which the variances {σ.sup.2(acc.sub.θ.sub.1), . . . , σ.sup.2(acc.sub.θ.sub.n), σ.sup.2(acc.sub.lz)} of the quantities stored in the rows of said time window AQ are calculated.

(61) It is here specified that the appending of a column takes place regardless of whether the relative measured accelerations correspond to a specific pattern such as the recognition of a walking step.

(62) Preferably, in said second appending operation, there is provided a control operation of the length of the time window AQ that includes eliminating the old column from the time window if the length (number of column) of said window exceeds a predetermined value, understood as the number of samples or as a number of seconds (then conveniently as a function of the columns number, because the number of columns is related to the sampling frequency of the accelerometer).

(63) Calculation of Angle θ.sub.dir Relative to the Axis of the Direction of Motion

(64) Advantageously, to optimize computational resources, said third operation of calculation of the variances of said fourth step 440 may include the use of traditional methods for on-the-fly and incremental calculation of the variance relevant to each angle θ.sub.j and to the signal relating to component acc.sub.lz. Some of these traditional methods are described by Chan, Tony F.; Golub, Gene H.; LeVeque, Randall J. (1983). “Algorithms for Computing the Sample Variance: Analysis and Recommendations”, The American Statistician 37, pages 242-247 and Ling, Robert F. (1974). “Comparison of Several Algorithms for Computing Sample Means and Variances”. Journal of the American Statistical Association, Vol. 69, No. 348, pages 859-866.

(65) In said preferred embodiment, one applies specifically to the present invention the physical principle of maximum variability of the linear measurement of the acceleration along the various possible directions. To do this, in said fifth step 450 of estimate of the angle θ.sub.dir, it is selected as an angle θ.sub.dir associated with the direction of the individual's motion the angle value θ.sub.j∈AS which has the maximum variance between said selected variances calculated in said third calculation operation of said fourth step of updating 440:
θ.sub.dir=arg max.sub.θ.sub.j.sub.∈ASσ.sup.2(acc.sub.θ.sub.j).

(66) Said fifth step of estimating the angle associated with the axis of the direction of motion 450 allows to obtain a direction θ.sub.dir selected from said set AS represented in FIG. 5; however, it does not yet allow to establish the travel direction along said direction axis. For example, if the selected direction axis is the one corresponding to the angle null, then it is not yet possible to determine which is the travel direction, which, with reference to FIG. 7, could be from the circumference center toward the indication of the arrow corresponding to the angle zero or from center of the circumference towards the indication of the corresponding arrow on the angle π. In a similar way, this holds for the other angles belonging to said set AS.

(67) Calculation of Angle θ (without Ambiguity on the Direction) Relevant to the Direction of Motion

(68) In said preferred embodiment, said sixth step (460) of ambiguity resolution of the direction is characterized by comprising: A first calculation operation of an index of “cross-correlation” (or “correlation with lag” that provides an correlation index as a function of the lag, which is defined as the deviation of a signal with respect to another) between two sequences extracted from the said time window AQ and respectively composed by the linear acceleration z component acc.sub.lz taken as the reference signal and the sequence of values acc.sub.θ.sub.dir, as obtained by extracting the row associated with θ.sub.dir from the reference time window AQ identified in the previous step; with respect to this row, lags are calculated; A second search operation that includes performing two searches starting from zero lag and respectively moving to the left by identifying the maximum correlation value of the left and moving to the right by identifying the maximum correlation value of the right. A third verification operation in which it is controlled in what position relative to zero lag the maximum of the left and right of the maximum are placed as found in said second search operation, and behaving as follows: If the maximum correlation value closest to zero lag is found on the left, then the direction θ.sub.dir chosen in the previous step correctly identifies also the direction, which corresponds to the value θ.sub.dir and proceeding to confirm θ=θ.sub.dir as estimate of the angle representing the axis and direction of the individual's motion. If the maximum correlation value closest to zero lag is found on the right, then the direction, related to θ.sub.dir, chosen in the previous step corresponds in the direction opposite to the direction of user's actual motion; in this case, one proceeds in estimating said direction by selecting as result the angle opposite to the angle θ.sub.dir estimated in said fourth step of calculation of the discrete angles; it follows that θ=θ.sub.dir+π.

(69) Preferably, the cross-correlation index is calculated using the Pearson correlation coefficient.

(70) Example of Determination of the Direction of Motion

(71) By way of non-exhaustive and non-limiting example, during the execution of an embodiment in accordance with the present method, the behavior of said fifth step of determining the direction is described when applied to a walk of an individual, who carries said device so that the angle θ to be estimated, represented in FIG. 4, is equal to +85° degrees or equal to +17π/36 radians. It is assumed that said initialization method of said embodiment example uses a value of n=18 by dividing accordingly the straight angle into 18 sectors of amplitude π/18 radians. It is also assumed that the size of said AQ time window is fixed to 2 seconds.

(72) FIG. 9a shows a graph of the acceleration detected on the X.sub.a, Y.sub.a, Z.sub.a axes as measured by said inertia measuring device during said exemplary walk. FIG. 9b shows said linear accelerations after execution of said third step 430 of linear acceleration calculation, wherein in said previous steps a digital low pass filter has been applied to the original signal measured by said device. As also described in the traditional method of WO2014/025429, it is evident that, with reference to the Z.sub.a axis of said figures, it is immediate to associate with each minimum of acceleration signal along Z.sub.a the occurrence of an event related to the beginning or the end walking step of said individual.

(73) FIG. 10a shows such acceleration values acc.sub.θ.sub.dir, i.e. the linear acceleration along the direction defined by the angle θ.sub.dir, and acc.sub.lz contained in a time window AQ.sub.0, during the execution of the method according to the invention, as observed immediately after the execution of said third operation of determining the direction of motion. Said figure compares the acc.sub.θ.sub.dir trend represented by the dotted line, which corresponds, in accordance with the present method, to the axis of the direction of motion defined by the angle

(74) ${\theta }_{dir}=\frac{4\pi }{9},$
with the linear acceleration measured along the Z.sub.a axis acc.sub.lz. In accordance with the execution of the method according to the invention, FIG. 10b shows the cross-correlation index calculated on two acceleration signals, respectively along the direction of motion as given by angle θ.sub.dir and along the Z axis of the reference system taken as absolute, highlighting with an arrow the point at which the maximum correlation is detected and with a vertical hatch the correspondence with the zero lag. In this case, said second verification operation of said fifth step of direction calculation confirms the estimate equal to

(75) $\frac{4\pi }{9}$
also for the direction of the axis of motion of said individual.

(76) Similarly, FIG. 11a shows said linear acceleration values acc.sub.θ.sub.dir, and acc.sub.lz contained in another specific exemplary time window AQ.sub.1 also observed immediately after the execution of said third operation of determining the direction of motion. Even in this case, said FIG. compares the trend represented with dotted line of the linear acceleration relative to said discrete angle that has a maximum variance and which corresponds in this example, in accordance with the present method, to the angle associated with the direction axis of motion

(77) 0${\theta }_{\mathrm{dir}}=-\frac{\pi }{2},$
with the linear acceleration acc.sub.lz measured along the Z.sub.a axis. In accordance with the execution of the method according to the invention, FIG. 11b shows the cross-correlation index calculated on the two acceleration signals respectively along the axis defined by the angle

(78) ${\theta }_{\mathrm{dir}}=-\frac{\pi }{2}$
and along the gravitational axis, highlighting with an arrow the point at which the maximum correlation is detected and with a vertical hatch the correspondence with the zero lag. In this case, said second verification operation of said fifth step of calculation of the direction detects a reversal of direction and produces as a result no longer the value of

(79) $-\frac{\pi }{2},$
but the opposite value equal to

(80) $\frac{\pi }{2}.$

(81) Other embodiments of the method 400 may provide for direction verification steps 460 different but referable to the same concept, i.e. to controlling how much phase difference is found between said two signals acc.sub.θ.sub.dir and acc.sub.lz. In this case, if said phase shift indicates a delay between said first sequence and said second sequence in the order, then θ=θ.sub.dir, otherwise θ=θ.sub.dir+π.

(82) In a possible invention embodiment, this has been advantageously verified by the use of the calculation of the cross correlation. The test of possible phase difference between said two signals corresponds substantially to the confirmation (if the direction is confirmed) or invalidation (if instead it is necessary to choose the opposite direction) of the hypothesis concerning the fact of using a right-handed reference frames for all of said absolute reference systems, individual and device identified in the present method.

(83) Also other embodiments of the method 400 may provide linear acceleration steps of calculation 410, of the window 440 update, for example by rotation of the raw acceleration value by a rotation matrix or by the use of a verification criterion of the direction, different from the use of cross-correlation, or by the use of other equivalent techniques preserving the general technical concept of the invention.

(84) Specific Technical Advantages in Using the Formalism of Quaternions

(85) Experimentally, the Inventors have verified that in some specific configurations of roll and pitch, the levelling by the defined quaternion q.sub.n further maximizes the variance along the direction of motion (decreasing the noise in the process of projection of the measured acceleration along the axis of motion), thereby improving the estimation of the angle θ.sub.dir according to the method proposed by the present invention.

(86) In particular, let us consider for example the numeric data, generated by a simulation using the simulation environment Matlab/Simulink and the 6DOF block (Quaternion) in aerospace Blockset™, relevant to a walk with the phone constantly kept in the attitude

(87) $\mathrm{conditions}\mathrm{roll}=\frac{\pi }{6}\mathrm{rad},\mathrm{pitch}=\frac{\pi }{4}\mathrm{rad}\mathrm{and}$$\mathrm{yaw}=\frac{\pi }{2}\mathrm{rad}$
with motion axis angle equal to

(88) ${\theta }_{dir}=\frac{\pi }{3}\mathrm{rad}.$
Using a discretization of straight angle with resolution

(89) $\Delta \theta =\frac{\pi }{10}\mathrm{rad},$
the estimation error on the angle θ.sub.dir, by using, respectively, in the step of acc.sub.h acceleration determination, a standard rotation according to attitude A.sub.d and a rotation linked to the proposed quaternion q.sub.n, comes out to be smaller in the second case (see FIG. 12). This experimental result shows the maximization of the variance along the direction of motion using the quaternion q.sub.n to obtain the acceleration acc.sub.h.
Novelty and Advantages

(90) In extreme synthesis, compared to the described conventional methods, the method subject-matter of the present invention presents the following differences: It allows to estimate the direction of motion of an individual or of a subject without tying said estimate to particular phenomena such as the recognition of a walking step or other type of pattern; It allows to set the desired degree of accuracy by defining the number of sectors in which the straight angle is split; It does not use heuristics for determining the direction of the actual direction of motion as other traditional methods do; The possibility of setting a priori a level of accuracy in the estimation allows also to adjust the computing power necessary for the execution of the method.

(91) From these differences the following advantages stem: Independence of the correctness of the estimate with respect to the proper functioning of other methods such as those relating to the recognition of patterns; Possibility to adjust the desired precision for the estimation of the direction and consequently the computational resources to be used; Estimate is not depending on heuristics that in certain situations could lead to incorrect results.

(92) The method of the above-mentioned patent document US2012136573A1 could provide a not exact projection onto the horizontal plane (because it is based on the hypothesis, not verified on-the-fly by the method, of coincidence between the gravitational axis and the axis of maximum variance in the accelerometer data). Such a problem does not occur in the method proposed in the present invention since the horizontal plane over which the individual moves is determined in an analytical way starting from the inertial data. In addition, the operation is computationally less expensive than the double exhaustive research proposed in the patent reported as prior art.

(93) As regards the disambiguation of the direction, the proposed method is able to provide an estimate for each time window analyzed unlike the prior art that in many cases fails to update the estimate.

(94) The thus conceived invention is susceptible to numerous modifications and variants, all falling within the scope of the appended claims.

(95) Moreover, all the details may be replaced with other technically equivalent elements.

BIBLIOGRAPHY

(96) SOH Madgwick et al. in “Estimation of IMU and MARG orientation using a gradient descent algorithm” published in proceeding of the IEEE International Conference on Rehabilitation Robotics (ICORR), 2011, ISBN 978-1-4244-9863-5. F. L. Markley, J. L. Crassidis, Fundamentals of Spacecraft Attitude Determination and Control, Springer, 2014. J. B. Kuipers, Quaternions & Rotation Sequences: A Primer with Applications to Orbits, Aerospace and Virtual Reality, Princeton University Press, 1999.

(97) 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 a limitation to the interpretation of each element identified, purely by way of example, by such signs and reference numbers.