STEP-BASED POSITIONING

Abstract

In a positioning system (100) for estimating the position of a mobile device (7), a processing system (9) receives external-range data, representative of a range between the mobile device and an external unit (2, 3, 4, 5), and acceleration data representative of acceleration of the mobile device due to its movement as it is carried by a person (6). The acceleration data is processed in a step-detection algorithm to determine step-distance data representative of a time series of step-data-based distances travelled by the mobile device, and step-distance data is processed to determine a step-data-based position estimate for the mobile device. A position estimate for the mobile device is determined by solving an optimisation problem comprising a first cost term based on distance to positions located at said range from the external unit, and a second cost term based on distance to the step-data-based position estimate or to positions located at a step-data-based distance from said step-data-based position estimate.

Claims

1. A method of estimating the position of a mobile device, the method comprising: receiving external-range data representative of a range between the mobile device and an external unit, wherein the external-range data is determined using a signal travelling between the external unit and the mobile device; receiving acceleration data, determined from an accelerometer of the mobile device, wherein the acceleration data is representative of acceleration of the mobile device due to a movement of the mobile device as it is carried by a person; processing the acceleration data in a step-detection algorithm to determine step-distance data representative of a time series of step-data-based distances travelled by the mobile device; processing the step-distance data to determine a step-data-based position estimate for the mobile device; and determining a position estimate for the mobile device by solving an optimisation problem for an objective function comprising a first cost term that depends on distance to positions located at said range from the external unit, and a second cost term that depends on distance to said step-data-based position estimate or to positions located at a step-data-based distance from said step-data-based position estimate.

2. (canceled)

3. (canceled)

4. (canceled)

5. (canceled)

6. (canceled)

7. (canceled)

8. (canceled)

9. The method of claim 1, wherein the second cost term depends on distance to positions located at a step-data-based distance from said step-data-based position estimate, and wherein the step-data-based distance is determined from one or more preceding step-data-based position estimates.

10. The method of claim 9, wherein the second cost term is weighted in the objective function at least in part in dependence on a variance of the step-data-based distance over time.

11. The method of claim 1, wherein the second cost term depends on distance to said step-data-based position estimate, and wherein the second cost term is equal or proportional to a Mahalanobis distance to the step-data-based position estimate.

12. The method of claim 1, wherein determining the step-data-based position estimate comprises processing an initial position estimate, determined from an earlier solving of an optimisation problem, to update the initial position estimate based on the time series of step-data-based distances.

13. The method of claim 1, comprising combining the heading data for the mobile device with the step-distance data to determine a motion vector or series of motion vectors, and using the motion vector or series of motion vectors to determine the step-data-based position estimate.

14. The method of claim 1, wherein the acceleration data comprises a time series of acceleration values, and wherein processing the acceleration data in the step-detection algorithm to determine the step-distance data comprises fitting a periodic function to the time series of acceleration values to determine a characteristic frequency representative of a step frequency of a person carrying the mobile device.

15. (canceled)

16. (canceled)

17. A method of determining a distance travelled by a mobile device, the method comprising: receiving first acceleration data, determined using an accelerometer of the mobile device, representative of a movement of the mobile device as it is carried by a person, the first acceleration data comprising a time series of first acceleration values; fitting a periodic function to the time series of first acceleration values to determine a characteristic frequency; receiving second acceleration data, determined using the accelerometer of the mobile device, representative of a further movement of the mobile device, the second acceleration data comprising one or more second acceleration values; determining, from the characteristic frequency and the second acceleration data, a step-completeness value representative of a non-integer fraction of a step taken by the person carrying the mobile device; and processing the step-completeness value and a step length parameter for determining a distance travelled by the mobile device due to the further movement.

18. The method of claim 17, wherein the periodic function is a sine function.

19. The method of claim 17, wherein the step length parameter represents a typical or recent step length of the person carrying the mobile device.

20. The method of claim 17, comprising filtering the first acceleration data with a low-pass frequency filter before fitting the periodic function, wherein the filter has a threshold frequency that corresponds to an expected minimum frequency of steps taken by the person carrying the mobile device.

21. The method of claim 17, comprising processing the first acceleration values to remove a contribution due to an effect of gravity on the mobile device before the periodic function is fitted to the acceleration data.

22. The method of claim 17, comprising calculating an acceleration parameter from the periodic function and determining the step-completeness value from the acceleration parameter, wherein calculating the acceleration parameter comprises evaluating the periodic function for a time instant associated with the second acceleration data.

23. The method of claim 17, wherein determining the distance travelled by the mobile device due to the further movement comprises calculating the product of the step completeness value and the step length parameter.

24. The method of claim 17, comprising calculating further values of distance travelled by the mobile device in response to receiving further acceleration data from the accelerometer due to additional movement of the mobile device.

25. The method of claim 17, comprising calculating a sum of a succession of step-completeness values to determine a total fractional step count, representative of an accumulated number of steps taken by the person carrying the mobile device.

26. The method of claim 25, comprising: using step-completeness values calculated at successive times to determine a time series of total fractional step count values; filtering and/or smoothing the time series of total fractional step count values; and using the time series of total fractional step count values to determine a time series of accumulated distances travelled by the mobile device.

27. The method of claim 26, comprising applying a linear regression operation, over time, to successive sliding time windows of the time series of total fractional step count values to generate a series of linear-regression values, and determining a maximum value of the characteristic frequency, for use when fitting the periodic function, from the slope of the series of linear-regression values.

28. The method of claim 17, comprising determining an accumulated distance travelled by the mobile device by calculating a sum of a succession of values of distances travelled by the mobile device, or as the product of a total fractional step count, representative of an accumulated number of steps taken by the person carrying the mobile device, and the step length parameter.

29. A system for determining a distance travelled by a mobile device, the system comprising a processing system configured to determining the distance travelled by the mobile device by: receiving first acceleration data, determined using an accelerometer of the mobile device, representative of a movement of the mobile device as it is carried by a person, the first acceleration data comprising a time series of first acceleration values; fitting a periodic function to the time series of first acceleration values to determine a characteristic frequency; receiving second acceleration data, determined using the accelerometer of the mobile device, representative of a further movement of the mobile device, the second acceleration data comprising one or more second acceleration values; determining, from the characteristic frequency and the second acceleration data, a step-completeness value representative of a non-integer fraction of a step taken by the person carrying the mobile device; and processing the step-completeness value and a step length parameter for determining a distance travelled by the mobile device due to the further movement.

30. (canceled)

Description

BRIEF DESCRIPTION OF THE DRAWINGS

[0103] Certain preferred embodiments of the invention will now be described, by way of example only, with reference to the accompanying drawings, in which:

[0104] FIG. 1 is a perspective diagram of a positioning system embodying the invention;

[0105] FIG. 2 is a schematic drawing of a static transmitter unit and a mobile device of the positioning system;

[0106] FIG. 3 is a flow chart showing operations carried out by a step detection algorithm embodying the invention;

[0107] FIG. 4 is a plot of acceleration measurements determined by the mobile device as part of the step detection algorithm;

[0108] FIGS. 5a and 5b illustrate an accumulation of steps using the step detection algorithm, and the associated acceleration measurements determined by the mobile device;

[0109] FIGS. 6a and 6b illustrate an accumulation of steps using the step detection algorithm;

[0110] FIG. 7 is a flow chart showing operations carried out by the positioning during a heading estimation process;

[0111] FIG. 8 is a flow chart showing operations carried out by the positioning according to a pedestrian dead reckoning process;

[0112] FIG. 9 is a flow chart showing operations carried out by the positioning as part of a synthetic-range based position estimation process; and

[0113] FIG. 10 is a flow chart showing operations carried out by the positioning as part of a Mahalanobis distance based position estimation process.

DETAILED DESCRIPTION

[0114] FIG. 1 shows part of a positioning system 100 that may be used in, for example, a shopping mall in order to determine the locations of shoppers within the shopping mall. Of course, this is just one example environment, and the positioning system could also be used in warehouses, hospitals, domestic homes, vehicles, etc.

[0115] FIG. 1 shows a room 1, to the walls of which are affixed four static transmitter units 2, 3, 4, 5. A person 6 in the room is carrying a mobile device 7. A network cable 8 connects each transmitter unit 2, 3, 4, 5 to a server 9, which is typically located in another room or in another building. These components cooperate to provide a positioning system 100, capable of estimating the three-dimensional location of the mobile device 7 within the room 1. In practice, the system 100 may have further similar transmitter units, installed throughout a building or series of rooms, and a plurality of similar mobile receiver units attached to, or incorporated into, people, animals, vehicles, robots, stock, equipment, etc.

[0116] FIG. 2 shows more detail of the mobile device 7 and a representative one 2 of the transmitter units 2, 3, 4, 5. The transmitter unit 2 has an ultrasonic sounder 213, a controller 215 for causing the transmitter device 2 to transmit ultrasonic signals, and a battery 217 for supplying power to the transmitter unit 2. The mobile device 7 has a microphone 203 capable of receiving ultrasonic signals from the transmitter unit 2, a microcontroller unit (MCU) 205, a battery 207 for powering the mobile device 7, a separate memory 208, and an inertial measurement unit (IMU) 209.

[0117] The inertial measurement unit 209 of the mobile device 7 comprises an accelerometer 210 and a gyroscope 211, with associated control circuitry. It is configured to sense acceleration along three axes, and rotation around three axes.

[0118] The microcontroller 205 executes software instructions stored in a memory of the microcontroller 205 or on the off-chip memory 208. It is configured to receive analogue or digital electronic signals representing acoustic signals received at the microphone 203, as well as rotation and acceleration signals from the IMU 209, and to process these signals.

[0119] The transmitter units 2, 3, 4, 5, and the mobile device 7 may have further standard electronic components such as radio transceivers, wired-network interfaces, display screens, buttons, etc. In some embodiments, the mobile device 7 is a tablet or mobile telephone (cellphone) such as an Apple or Android smart phone.

[0120] The microcontroller units 205, 215 can include one or more processors, DSPs, ASICs and/or FPGAs. They can include memory for storing data and/or for storing software instructions to be executed by a processor or DSP. The microcontroller 205, 215 or units 2, 7 more generally can include any other appropriate analogue or digital components, including oscillators, ADCs, DACs, RAM, flash memory, etc.

[0121] The transmitter units 2, 3, 4, 5 are configured to transmit, at regular intervals (e.g. every one second), a locating signal comprising a signature unique to that transmitter unit. The transmissions may be coordinated, e.g. by the server 9, to minimise negative interference between nearby transmitter units. The transmitted signals may be received at the microphone 203 of the mobile device 7. The received signals may then be demodulated and processed by the microcontroller unit 205 of the mobile device 7 to identify signatures that have been transmitted by transmitter units 2-5 within audible range of the mobile device 7 and to determine one or more times of arrival for each signature.

[0122] The receiver unit 7 uses an internal clock to timestamp the received signatures. The receiver unit 7 may be synchronised with each of the transmitter units 2-5 (and optionally with the server 9) and can therefore also calculate a time of flight (TOF) for each signature it receives. Based on the TOF measurement, and the known location of the transmitter units 2-5, a range between the mobile device 7 and the identified transmitter unit may can be determined. By combining three or more TOF measurements from known transmitter locations, the position of the mobile receiver unit 7 can be determined using geometric principles of trilateration (also referred to as multilateration). If the receiver unit 7 is not synchronised with the transmitters 2-5, a time difference of flight (TDOF) method can still be used to determine the position of the mobile device 7; however, in this case, more transmitter units 2-5 may have to be in range of the mobile device 7 to get an accurate position estimate.

[0123] New TOF distances may be obtained repeatedly over time, which can result in an over-specified problem. Individual measurements may be inaccurate, e.g. due to noise, and can also become out-of-date when the mobile device 7 is in motion. One approach that the microcontroller 205 uses to estimate the position of the mobile device 7 is to use a minimization algorithm to generate a position estimate, r, of the mobile device 7 using a cost function, J.sub.1, with the following form:

[00001]$J_1\left(r;r_i,R_i\right)=\frac{.Math.w_i\left({\left(r-r_i\right)}^2-R_i^2\right)}{.Math.{\mathrm{rel\_w}}_i}$

where: [0124] the sums are over all the transmitter units, i, from which signals have been received, [0125] r is a vector representing the position estimate to be determined for the mobile device 7 (in three-dimensional space or two-dimensionally in a horizontal plane), [0126] r.sub.i is the known 2D or 3D position of the transmitter unit i, and [0127] R; is the latest measured range between the mobile device 7 and the transmitter unit i.

[0128] Each of the transmitter units 2, 3, 4, 5 i is given a weighting, w.sub.i. The weight, w.sub.i, may include an absolute component, e.g. which depends on the variance (or alternative error term) of the measurements for a single transmitter unit i. It may also optionally include a relative component, rel_w.sub.i, which may represent a confidence in the reliability of a measurement from a transmitter unit i relative to the other transmitter units; this may, for example, depend on the age of the measurement (i.e. with a lower relative confidence score being assigned to older measurements).

[0129] The cost function J.sub.1 is normalised for the relative weights, rel_w.sub.i, but not for the absolute weight components. This is because the absolute component (e.g. the variance in measurements) indicates how reliable each transmitter unit is, so should be taken into account globally in the cost function, whereas the relative weights only indicate which transmitter units are trusted more, relative to the other transmitter units, not the reliability of the measurements in an absolute sense.

[0130] In this way, the positioning system 100 shown in FIG. 1 may be used to determine the position of the mobile device 7 within rooms of a building in which static transmitter units 2, 3, 4, 5 are located.

[0131] While this approach may be used to estimate the position of the mobile device 7, the applicant has recognised that, if an initial position estimate of the mobile device 7 is known, such as from ultrasound-signal TOF data as described above, subsequent position estimates can be improved using additional data obtained from the inertial measurement unit (IMU) 209 of the mobile device 7. IMU data associated with the motion of the mobile device 7 may be processed using pedestrian dead reckoning techniques to determine additional range estimates, which may be combined with the position data determined based on ultrasound TOF measurements to improve the position estimate of the mobile device 7. As long as an initial positon estimate of the device is known, IMU data can also be used to determine an updated position of the mobile device 7 even when it is out of range of the transmitter units 2, 3, 4, 5.

[0132] Thus a first position estimate for the mobile device 7 may be generated when in range of the transmitter units 2, 3, 4, 5 using TOF measurements as described above, and updated position estimates of the mobile device 7 may be determined using one or more further range estimates determined from data output by the IMU 209.

[0133] In order to determine a range estimate based on measurements by the IMU 209, data from the accelerometer 210 and the gyroscope 211 of the IMU 209 may be processed using a step detection algorithm as will be described in the following.

[0134] Step Detection Algorithm Outline

[0135] When the mobile device 7 is being carried by a walking person 6, the IMU 209 is in a dynamic state as a result of the walking motion, and this motion can be detected by the accelerometer 210 of the IMU 209 through periodic fluctuations in the acceleration signal. The applicant has recognised that due to the periodic nature of the acceleration changes, if the contribution due to gravity is removed, the acceleration signal associated with a walking motion can be well approximated by a sine wave, characterised by a particular walking frequency, and that this can be tracked using a step detection algorithm, executed by the microcontroller 205 of the mobile device 7.

[0136] The step detection algorithm disclosed herein can be used to determine the number of steps taken by a user 6 by approximating the norm of the gravity-free acceleration vector as a sine wave characterised by a walking frequency .sub.. A left-leg step and a right-leg step are here considered two individual steps, and the walking frequency is the frequency of these individual steps.

[0137] FIG. 3 shows a process flow diagram depicting the steps taken by the step detection algorithm.

[0138] In step 300, the microcontroller 205 samples acceleration data (in the form of an acceleration vector) measured by the accelerometer 210 of the IMU 209 to be processed using the step detection algorithm. The acceleration vector may be sampled at regular time intervals (e.g. every 30 ms), or sampling may be carried out whenever new data is received at the microcontroller 205 from by the IMU 209 (i.e. based on the update frequency of the IMU 209).

[0139] In step 302, the step detection algorithm applies a low pass filter to the acceleration vector measured by the accelerometer 210 of the mobile device 7.

[0140] The power spectral density of the acceleration vector is approximately zero at lower frequencies, otherwise the walking user 6 would be in constant acceleration. Except for during transition periods (i.e., at the start of a walking motion), when true acceleration or deceleration occurs, a lower bound on the walking frequency can therefore be applied in the frequency domain, for example using a low pass filter. The lower bound on frequency can be determined by considering the average duration of each step taken by a person 6 carrying the mobile device 7. Starting from an initial position, the walking motion of the person 6 carrying the mobile device 7 can be described by a step cycle consisting of two steps (i.e. left foot then right foot). At the end of each cycle, the person's posture returns to a state approximately equivalent to that of the person's initial state. As a result, the acceleration profile measured by the IMU 209 approximately repeats every step cycle, i.e. every two steps. A lower bound value for the walking frequency .sub. can therefore be set as the reciprocal of half the duration of the step cycle. However, lower frequencies might still occur in practice due to slight changes in walking velocity and sensor movements unrelated to the walking motion.

[0141] A low pass filter is applied along each of the three accelerometer axes (referred to in the following as x, y and z axes) to determine acceleration components a.sub.x, a.sub.y, a.sub.z as a function of time t.

[0142] In step 304, by calculating the norm of the acceleration value along each axis, and removing the local gravity contribution, g, from the measured acceleration vector, the gravity-free accelerator norm, a.sub.n(t) is determined.

[0143] The local gravity contribution can be separated into a magnitude scalar g.sub.m and a unit vector (referred to in the following as the gravity vector), represented in the local frame of the device 7, that relates the device's frame of reference to the horizontal plane of the North-East-Down (NED)-frame in local tangent plane coordinates. The gravity-free accelerator norm, a.sub.n(t) can therefore be given by

a.sub.n(t)={square root over (a.sub.x(t).sup.2+a.sub.y(t).sup.2+a.sub.z(t).sup.2)}g.sub.m.Math.

[0144] FIG. 4 shows an example of the gravity-free accelerometer norm 401 (and its derivative with respect to time 403), over time, determined using acceleration measurements from the accelerometer 210 of a mobile device 7 being carried by a walking person. It can be seen in FIG. 4 that the gravity-free accelerometer norm is periodic in nature, and can hence be approximated using a periodic function, e.g. a sine function.

[0145] This approximation is performed in step 306, using the function

a.sub.n(t)A sin(2.sub.t)

where .sub. is an estimated walking frequency and A is an unknown amplitude. The estimated walking frequency may initially be set at a value close to a typical average walking frequency of the person 6 carrying the mobile device 7 (typically between 0.5 and 2 Hz), but can be updated over time as more acceleration data is measured by the accelerometer 210 as will be explained in the following.

[0146] The derivative of the gravity-free accelerator norm, {dot over (a)}.sub.n, is also periodic in nature, and can be approximated by the cosine function

{dot over (a)}.sub.n(t)2.sub.A cos(2.sub.t)

which can be scaled towards the approximated gravity-free accelerator norm to ensure a smooth gradient in the following calculations as

[00002]${\stackrel{.}{a}}^{*}\left(t\right)=\frac{1}{2f_{}}\stackrel{.}{a}\left(t\right)$

[0147] The magnitude of the gravity-free accelerator norm, {dot over (a)}.sub.n is compared to a threshold acceleration value (e.g. 0.98 m/s.sup.2) in step 307. If the magnitude of the gravity-free accelerator norm, {dot over (a)}.sub.n is found to be below the threshold value, the acceleration measurement is not processed using the step detection algorithm.

[0148] In step 308, for each acceleration measurement, a.sub.k=a(t.sub.k), recorded by the accelerometer 210, a step completeness estimate, .sub.k.sup.c, is calculated using a complex number transformation, where the real part of the complex number is taken to be the magnitude of the acceleration, a.sub.n, and the imaginary part of the complex number is given by the scaled derivative of the acceleration, {dot over (a)}*. A phase value, corresponding to the step completeness estimate, .sub.k.sup.c, is then calculated using the two-argument arctangent function:

[00003]${}_k^c=\frac{1}{2}\arctan 2\left(\frac{a_k}{{\stackrel{.}{a}}_k^{*}}\right)$

which provides a step completeness estimate constrained to the interval [0.5, 0.5]. The step completeness estimate is a value which represents the fraction of the step taken by the person 6 carrying the mobile device 7 at the time at which the acceleration measurement is recorded by the accelerometer 210.

[0149] Having calculated the step completeness estimate associated with a given acceleration measurement, the distance travelled by the person 6 carrying the mobile device 7 over multiple accelerometer measurements can be estimated by calculating a total fractional step count, .sub.k, representing an accumulation of steps having both an integer component and a fractional component.

[0150] The total fractional step count, calculated in step 310, is an accumulated value updated using the step completeness estimate, .sub.k.sup.c calculated in response to each accelerometer measurement, i.e. the total fractional step count can be calculated according to:

.sub.k=.sub.k.sup.c[.sub.k.sup.c.sub.k1](1)

where .sub.k1 is the value of the total fractional step count calculated in response to the previous accelerometer measurement, and the square parenthesis indicate rounding to the nearest integer. In this way, the total fractional step count, .sub.k, is continually updated as more accelerometer measurements are recorded as the person 6 carrying the mobile device 7 continues to walk or run. In general, the value of the total fractional step count increases as the person 6 continues to take steps, however it is not prohibited from decreasing a small amount in response to certain accelerometer measurements, as shown in FIGS. 5a and 5b.

[0151] FIG. 5a shows a second example of the gravity-free accelerometer norm 401 over time, determined using acceleration measurements from the accelerometer 210 of a mobile device 7 being carried by a walking person 6. Its derivative with respect to time is shown as a dashed line 403. The gravity-free accelerometer norm 401 can be seen to be periodic in nature, and, in the example shown, comprises two complete step cycles 405a, 405b, each corresponding to two steps being taken by the walking person 6.

[0152] FIG. 5b shows the estimated total fractional step count (.sub.k) 409 over the same time frame as FIG. 5a, calculated based on the measured acceleration data using the step detection algorithm, calculated as discrete values at successive time instants (the larger diamonds marks in FIG. 5a). It can be seen in FIG. 5b that, as successive acceleration measurements are recorded, the total fractional step count is generally increasing. It can also be seen that for each step cycle, i.e. for every two steps taken by the user, the total fractional step count increases by approximately two. The total fractional step count shown in FIG. 5b can also be seen to be discontinuous in nature. Such discontinuities in the total fractional step count occur in the event that the gravity-free accelerator norm, {dot over (a)}.sub.n associated with an accelerometer measurement is found to be below a threshold value, and hence is not processed, as described above in relation to FIG. 3.

[0153] Also shown in FIG. 5b, using smaller diamond marks, is a filtered total fractional step count 407. This is calculated from the unfiltered count 409 by applying linear regression over the preceding N samples (for some appropriate integer N, such as N=sample rate4 seconds), such that each data point of the line 407 is calculated by applying a respective linear fit to the last N data points of the line 409. This process is illustrated further in FIG. 6a.

[0154] FIG. 6a shows a later portion of the total fractional step count, calculated using acceleration vectors a.sub.k at sequential times t, represented by the line 501. This line 501 is a continuation of the line 409 from FIG. 5b.

[0155] Also shown in FIG. 6a is a single linear fit line 503 for the total unfiltered fractional step count 501, determined for the latest time point, t=30. This linear fit is calculated by storing the last N values of the total fractional step count [.sub.i, .sub.i+1, . . . .sub.k] (where i=kN1) in a buffer in the memory 208 of the mobile device 7, and applying a linear fit (shown as the line 503 in FIG. 6a) to (t) over time using the fitting parameters {tilde over ()}.sub.0 and {tilde over ()}.sub., i.e.

[.sub.i,.sub.i+1, . . . .sub.k].fwdarw.(t)={tilde over ()}.sub.0+{tilde over ()}.sub.(t)

[0156] In this particular example, N is chosen to span four secondsi.e. N=sample rate4 secondsbut other values may be used instead. The value of this single linear-fit line 503, at time t=30, is then used as the value for the filtered total fractional step count 407 for time t=30. Successive points of the filtered count 407 over time are determined similarly, i.e. by performing a corresponding succession of linear-fit operations.

[0157] In this way, a series of total fractional step count values may be used to determine an updated walking frequency {tilde over ()}.sub. (in optional step 312 of FIG. 3) as the slope of a linear fit, which can be used in the next iteration of the calculation of the total fractional step count. An improved calculation of the step completeness estimate can then be given by .sub.k={tilde over ()}.sub.0+{tilde over ()}.sub.t.sub.k (where the bar symbol denotes an estimated or predicted value), which may be used as the value of the total fractional step count in the next iteration of calculation of the total fractional step count.

[0158] FIG. 6b shows a comparison of a calculated value of the filtered total fractional step count 505 (calculated as described for the filtered count 407 above), determined for a person 6 carrying a mobile device 7 configured to execute the step detection algorithm described above, and a plot 507 of true steps taken by the person 6, as recorded directly using a motion-capture system. The true step is based on movement of the foot of the person 6 immediately before breaking contact with the floor, at which point the step is recorded. Such a motion-capture system would not normally be used in practice; this comparison is provided here, rather, to demonstrate the efficacy of the disclosed step counting methodology. It can be seen in FIG. 6b that, by approximating the acceleration of the mobile device 7 as a periodic function, and applying a continuously-updated linear-fit process to the total fractional step count over time, near real-time or even predictive step estimation can be achieved. This is in contrast to prior art approaches, in which steps are not counted until after the motion has finished.

[0159] In some embodiments the step detection algorithm may be further improved by setting appropriate boundary conditions such as limiting the minimum or maximum amplitude based on the detected acceleration (e.g. by setting A.sub.max={square root over (a.sup.2+{dot over (a)}.sup.2)}), or by setting a maximum value for the estimated walking frequency .sub. determined using the linear fit.

[0160] In order to estimate the distance travelled by the person 6, the total fractional step count, .sub.k, is combined with step length datae.g. multiplied by a standard step length. A default step length may be used, at least initially (e.g. based on the height of the user), or an average step length may be calculated by the device 7 based on external measurementse.g., by counting the number of steps detected when the mobile device 7 is in range of the ultrasound transmitters 2, 3, 4, 5, and dividing this by the distance travelled as determined by the ultrasound TOF measurements. After an extended walking period, small changes in step length may be estimated by the device 7 from the changes in the acceleration amplitude over time.

[0161] By fusing distance estimates determined using the step detection algorithm with heading data determined using gyroscope 211 of the IMU 209, it is possible for the controller 205 to determine an updated position estimate of the person 6 carrying the mobile device 7 after each update from the IMU 209 using a pedestrian dead reckoning approach as will be described further below.

[0162] Heading Estimation Outline

[0163] In order to estimate the heading of a person 6 carrying the mobile device 7 using gyroscope measurements, it is assumed that the change in the horizontal rotation of the mobile device 7 is equivalent to the change in heading of the person 6 carrying the mobile device 7. This assumption is valid so long as the relative position of the mobile device 7 and the person 6 carrying the mobile device 7 is constant, such as when the mobile device 7 is carried within a pocket, or on the wrist of the person 6. Although the absolute value of the horizontal orientation of the mobile device 7 (i.e. its yaw) may be different to that of the person 6 carrying the mobile device 7, the change in the yaw of the mobile device 7 and in the heading of the person 6 carrying the mobile device 7 coincide. Depending on how it is carried, the mobile device 7 may be subject to certain additional rotations (e.g. due to the motion of the person's leg in the case of the mobile device 7 being held in a pocket), however the overall change in direction corresponds to the change in heading of the person 6 carrying the mobile device 7. In this way, if an initial absolute orientation of the mobile device 7 is known, the heading of the person 6 carrying the mobile device 7 can be determined.

[0164] FIG. 7 illustrates how the heading of the person 6 carrying the mobile device 7 is measured over time using a heading estimation process. The user heading estimation process is based on determining, for each new measurement, k, received from the IMU 209, the horizontal component of the rotation rate .sub.k as measured by the gyroscope 211, where horizontal indicates the plane orthogonal to the Down vector in the North-East-Down (NED)-frame.

[0165] The horizontal rotation rate of the mobile device 7, i.e. the rate at which the yaw of the mobile device 7 changes, can be extracted from the raw gyroscope data by applying the gravity vector at the time of measurement, .sub.k, to each incoming measurement of the rotation rate .sub.k. The rate of change of the heading, .sub.k, of the person 6 carrying the mobile device 7 (i.e. the user heading) can therefore be determined by calculating the dot product of the rotation rate, .sub.k and the gravity vector .sub.k for each measurement carried out by the IMU 209.

[0166] The heading estimation process consists of an initialization phase, a propagation phase and a correction phase, as will be described in the following in relation to FIG. 7.

[0167] In the initialisation phase (step 600) an initial user heading .sub.0 (with variance P.sub..sub.0), an initial gravity vector .sub.0 (with 33 covariance matrix P.sub.0), and a sensor noise term Q=I.sub.3.sub.t are set (where I.sub.3 is the 33 identity matrix, .sub. is a vector containing the standard deviation over the three axes of the gyroscope 211, and t is the interval between successive gyro measurements or IMU updates).

[0168] The initial user heading .sub.0 may be estimated from on a pair of position estimates, taken at different times, which may be determined from an external system, such as by processing ultrasound signals received from the transmitters 2-5, or from a compass, or GPS data, or in any other way. An initial variance P.sub..sub.0 (an error term) may be determined similarly, or may be set to infinity if unavailable. In some embodiments, the initial user heading may instead be set to zero if it is sufficient to track only relative heading over time.

[0169] In some preferred embodiments, the initial gravity vector .sub.0 is determined by averaging data from the accelerometer 210 over a prolonged static periodi.e. a time when the device 7 is stationary. More advanced methods to estimate the gravity vector under dynamic conditions may be used in some embodiments.

[0170] The standard deviation of the gyroscope 211 over three axes may be determined empirically, e.g. by the controller 205 analysing the gyroscope output over time.

[0171] In the propagation phase, starting in step 602, the microcontroller 205 samples the rotation rate .sub.k measured using the gyroscope 211 of the IMU 209. The heading vector and the gravity vector are then estimated for each incoming rotation rate measurement. This is achieved by applying the gravity vector .sub.k to the measured rotation rate .sub.k and carrying out a forward integration over the time period t between incoming measurements to determine the change in user heading .sub.k, in step 604, i.e.

.sub.k=.sub.k.sup.T.sub.kt

[0172] The variance of the change in user heading .sub.k is based on the sensor noise Q and the covariance of the gravity vector, and can be given by

P.sub..sub.k=.sub.k.sup.TQ.sub.k+.sub.k.sup.TP.sub..sub.k.sub.k

[0173] Having calculated the change in user heading and the variance of this change, an updated heading estimate .sub.k+1 of the person 6 carrying the mobile device 7 and an updated variance in the user heading P.sub..sub.k+1 are calculated, in step 606, by

.sub.k+1=.sub.k+.sub.k

P.sub..sub.k+1P.sub..sub.k+P.sub..sub.k(2)

[0174] In order to update the gravity vector, .sub.k, a rotation matrix, R.sub.k is generated based on the rotation rate .sub.k, and is calculated according to:

[00004]$R_k=I+\sin \left(.Math.{}_k.Math.t\right)\left[\frac{{}_k}{.Math.{}_k.Math.}\right]+1-{\cos \left(.Math.{}_k.Math.t\right)\left[\frac{{}_k}{.Math.{}_k.Math.}\right]}^2$

[0175] By applying the rotation matrix R.sub.k to the gravity vector, the updated value of the gravity vector, .sub.k+1, can be determined (in step 608) from the previous value .sub.k as

.sub.k+1=R.sub.k.sub.k(3)

and the covariance can be updated as

P.sub..sub.k+1=R.sub.kP.sub.kR.sub.k.sup.T+[.sub.k].sub.xQ[.sub.k].sub.x.sup.T

where [.sub.k].sub.x is the skew-symmetric matrix of [.sub.k], i.e. for [.sub.k]=(.sub.k.sub.1,.sub.k.sub.1,.sub.k.sub.1).sub.T, [.sub.k].sub.x can be given by:

[00005]${\left[{\stackrel{}{}}_k\right]}_{}=\left[\begin{array}{ccc}0& -{\stackrel{}{}}_{k_3}& {\stackrel{}{}}_{k_2}\\ {\stackrel{}{}}_{k_3}& 0& -{\stackrel{}{}}_{k_1}\\ -{\stackrel{}{}}_2& {\stackrel{}{}}_{k_1}& 0\end{array}\right]$

[0176] In this way, for each new measurement of rotation rate .sub.k output by the gyroscope 211, the user heading estimate .sub.k and the gravity vector .sub.k are updated, allowing the heading of the person 6 carrying the mobile device 7 to be determined and updated over time. Updated error terms, in the form of a covariance matrix, are also calculated. Although in the embodiment described in relation to FIG. 7, steps 606 and 608 are shown as being sequential, in some embodiments steps 606 and 608 may instead be run in parallel.

[0177] In some embodiments, the estimated user heading and gravity vector may be further adjusted or corrected (in step 610) by combining the updated values with additional data determined by the mobile device 7, using a Kalman filter. For example, heading data may be corrected using magnetometer data (e.g. from a magnetometer in the IMU 209) or from Doppler-shift data obtained from ultrasound signals received by the microphone 203 from one or more transmitter units 2-5.

[0178] Similarly, gravity vector corrections may be determined, based on additional data from the accelerometer 210, where low-variance measurements (e.g. obtained during static conditions) can give more accurate results.

[0179] In the event that such a correction is used, the corrected user heading and gravity vector values are used in the next iteration of the heading estimation process in place of the original updated values.

[0180] Pedestrian Dead Reckoning Model

[0181] Having outlined how step data is determined using IMU measurements and the step detection algorithm, and that a user heading estimate .sub.k may also be calculated using IMU measurements, the application of this data to a pedestrian dead reckoning (PDR) framework based on an extended Kalman filter (EKF) will be explained in the following.

[0182] In the PDR framework, a state vector, x, is defined, which contains 2D position p=(x,y), heading and step parameters. The state vector may be described by

[00006]$x=\left[\begin{array}{c}x\\ y\\ \\ \\ L\end{array}\right]$

where x, y represent the latest estimated 2D position of the mobile device 7, is the latest user heading, is the gravity vector and L is the step length of the person 6 carrying the mobile device 7.

[0183] The controller 205 may store data representing this state vector in the memory 208, and calculate an update to the vector at regular intervals (e.g. every time a new measurement is available from the IMU 208).

[0184] An initial estimate of the 2D position may be determined using ultrasound signals received from the transmitter units 2-5, with subsequent values updated using IMU data, as well as further ultrasound signals, as will be described in the following.

[0185] Since measurements from the IMU 209 come in at discrete time intervals t, the propagation of the state vector constitutes a discrete, non-linear estimation problem of the form:

x.sub.k+1=(x.sub.k,.sub.k,.sub.k,v.sub.k)

where the dynamic model is a function of the current state x.sub.k, a rotation rate measurement .sub.k, an acceleration measurement .sub.k, and sensor noise v.sub.k.

[0186] The dynamic model implemented by the controller 205 updates the step-cycle and heading estimates using the step-detection algorithm and heading-estimation process described above. It also evaluates a position-update equation incorporating distance values determined using the step detection algorithm.

[0187] The change in 2D position, as a function of total fractional step count and user heading, is given by

x.sub.k=.sub.kL.sub.k cos(.sub.k)

y.sub.k=.sub.kL.sub.k sin(.sub.k)(4)

where .sub.k=.sub.k+1.sub.k is the increment in the total fractional step count.

[0188] In this way, the dynamics of the state vector can be described by combining the outputs of the step-detection algorithm (implementing Equation 1), user heading process (implementing Equations 2, 3), and position-update algorithm (implementing Equation 4), resulting in the following dynamic model setup:

[00007]$f\left(x_k,{}_k,{}_k,v_k\right)=\left[\begin{array}{c}{x}_k\\ y_k\\ {}_k\\ 0\\ L_k\end{array}\right]+\left[\begin{array}{c}{}_k{L}_k\cos \left({}_k\right)\\ {}_k{L}_k\sin \left({}_k\right)\\ {}_k^T{}_{k}t\\ R\left({}_{k}t\right){}_k\\ 0\end{array}\right]$

[0189] In the dynamic model setup, the position and user heading equations are implemented directly, and .sub.k is the result of the step detection algorithm which determines the total fractional step count, based on the last N measurements provided by the accelerometer 210 as described above.

[0190] Having established the state vector, x, and the dynamic model, , the extended Kalman filter (EKF) setup of the pedestrian dead reckoning algorithm, implemented by the controller 205, will now be described.

[0191] FIG. 8 provides an overview of the steps performed.

[0192] The EKF setup of the PDR algorithm consists of an initialization phase, a propagation phase and a correction phase.

[0193] In the initialisation phase (step 700) initial values of the state vector, x, and a covariance matrix P.sub.0, based on uncorrelated noises of the initial state parameters (except for the 2D position), are set, based on measurements from the IMU 209, as:

[00008]${\stackrel{}{x}}_0=\left[\begin{array}{c}{\stackrel{}{x}}_0\\ {\stackrel{}{y}}_0\\ {\stackrel{}{}}_0\\ {\stackrel{}{}}_0\\ {\stackrel{}{L}}_0\end{array}\right],P_0=\left[\begin{array}{cccc}{P}_{p_0}& & .Math.& 0\\ & {}_{{}_0}^2& & .Math.\\ .Math.& & P_{{}_0}& \\ 0& .Math.& & {}_{L_0}^2\end{array}\right]$

where P.sub.p.sub.0 is a 2-by-2 covariance matrix of 2D position, .sub..sub.0.sup.2 is the user heading variance, P.sub..sub.0 is a 3-by-3 covariance matrix of the gravity vector, and .sub.L.sub.0.sup.2 is the step length variance. These initial values may be set in a similar way to the initialisation processes described above. An alternative way to initialize the state is to set the initial position {circumflex over (p)}.sub.0 and user heading {circumflex over ()}.sub.0 to zero, in which case the filter will estimate the relative position and user heading over time.

[0194] In the propagation phase (step 702), the state vector is updated as discrete data is received from the IMU 209 at regular time intervals t. The IMU data comprises a measured rotation rate .sub.k from the gyroscope 211 and an acceleration measurement .sub.k from the accelerometer 210.

[0195] The expected value of the propagated state, x.sub.k+1 as well as its covariance, P.sub.k+1, (i.e. an updated error term) can then be determined by applying the dynamic model (in step 702) as

x.sub.k+1=(x.sub.k,.sub.k,.sub.k)

P.sub.k+1=F.sub.kP.sub.kF.sub.k.sup.T+G.sub.kQG.sub.k.sup.T

where: [0196] i. F.sub.k is the Jacobian F of evaluated at x.sub.k. This Jacobian is composed of its individual components described in relation to the step detection (1), user heading (2, 3) and position estimation (4) equations, giving F.sub.k as:

[00009]$F_k=\left[\begin{array}{ccccc}1& 0& -{\stackrel{}{y}}_k& 0^T& {\stackrel{}{x}}_k/L_k\\ 0& 1& {\stackrel{}{x}}_k& 0^T& {\stackrel{}{y}}_k/L_k\\ 0& 0& 1& {}^{T}t& 0\\ 0& 0& 0& R\left({\stackrel{}{}}_{k}t\right)& 0\\ 0& 0& 0& 0^T& 1\end{array}\right]$ where 0 is a 31 vector of zeros. [0197] ii. G.sub.k is a matrix describing the relationship between the process noise and the state vector given by:

[00010]$G_k=\left[\begin{array}{c}{0}^T\\ 0^T\\ {\stackrel{\_}{}}_k^T\\ {\left[{}_k\right]}_{}\\ 0^T\end{array}\right]$ [0198] iii. Q is the 3-by-3 covariance matrix of the gyroscope 211 where [0199] Q.sub.11=(.sub..sub.xt).sup.2, Q.sub.22=(.sub..sub.yt).sup.2, Q.sub.33=(.sub..sub.zt).sup.2

[0200] In some embodiments in which the state vector is measured directly, the standard Kalman filter equations can optionally be applied for state corrections (in optional step 704). Examples of direct measurements include position measurements based on trilateration of ultrasound ranges (giving an xy-position {circumflex over (p)} with covariance P.sub.{circumflex over (p)}), average step length measurements based on travelled distance

[00011]$\left(e.g.\stackrel{}{L}=\frac{.Math.p_2-p_1.Math.}{{}_2-{}_1}\right),$

or user heading measurements ({circumflex over ()}) based on magnetic measurements and/or ultrasound Doppler-shift velocity vectors.

[0201] Thus, for each measurement output by the IMU 209, the state vector may be updated, providing an updated position estimate and heading of the person 6 carrying the mobile device 7, along with associated error terms, by implementing the step detection algorithm described above.

[0202] Combining Range Estimates

[0203] Having described how the pedestrian dead-reckoning state vector is used to provide a position estimate of the mobile device 7, the combining of this position estimate with TOF-based position measurements to more accurately determine the position of the mobile device 7 will now be described.

[0204] As described above, a position of the mobile device 7 may be determined, without using the IMU 209, by solving an optimization problem for range measurements between the mobile device 7 and the individual ultrasound transmitter units 2, 3, 4, 5 based on TOF measurements of ultrasound signals using a cost function, J.sub.1:

[00012]$J_1\left(r;r_i,R_i\right)=\frac{.Math.w_i\left({\left(r-r_i\right)}^2-R_i^2\right)}{.Math.{\mathrm{rel\_w}}_i}$

[0205] Combining such TOF position data with position data estimated using the step detection algorithm and pedestrian dead reckoning state vector described above can provide improved positioning accuracy. Moreover, a time series of position estimates may potentially then be updated as frequently as the update frequency of the IMU 209, which may be substantially higher than the rate of ultrasound signal transmissionse.g. more often than once per second (however, such faster updating is not essential, and the optimization calculation may be performed less often than on every IMU update in some embodiments). Incorporating dead-reckoning can also be useful for continuing to track the mobile device 7 accurately when it moves out of range of some or all of the ultrasound transmitters 2, 3, 4, 5 of the system 100.

[0206] The ultrasound positioning and the PDR process can be regarded as two separate processes, each based on its own set of sensor data and producing independent results. Therefore, finding the optimum position requires some form of data fusion between the two processes. This could be achieved, using a nave approach, by calculating a first independent position estimate from the ultrasound positioning system, and a second independent position estimate using the PDR process, and fusing the two position estimates together using the standard Kalman filter equations.

[0207] However, in embodiments of the invention, an improved position estimate may instead be determined by combining external-range data obtained using the external transmitters 2-5 with dead-reckoning data obtained from the internal IMU measurements in a single optimization process, e.g. by incorporating data from the PDR state vector described above in the ultrasound positioning cost function J.sub.1.

[0208] Two different exemplary approaches for doing this are described below.

[0209] The first approach is based on converting PDR data into a synthetic range measurement which can be directly incorporated into the ultrasound positioning cost function, J.sub.1, as defined above, resulting in a single augmented cost function J.sub.1* containing an additional range term, associated with a virtual transmitter, with properties that depend on values determined using the PDR process.

[0210] The second approach is to incorporate the latest position estimate determined using the PDR algorithm into an optimisation routine involving the ultrasound positioning cost term, J.sub.1, in the form of an additional cost term, J.sub.2, that encourages the overall position estimate to be close to the latest dead-reckoning position estimate. This term may be weighted by a scaling factor, . An improved position estimate for the mobile device 7 may then be determined, e.g. following each IMU update, by minimising the cost function, J.sub.3=J.sub.1+J.sub.2.

[0211] Approach IConverting to a Synthetic Range Measurement

[0212] Some embodiments implement a first method to combine position estimates determined using ultrasound measurements with those determined using the PDR approach described above, which involves converting the travelled distance, as estimated by the PDR algorithm, into a range measurement centred on the initial position of the current state vector. In effect, this creates a synthetic transmitter located at an origin corresponding to the initial position of the state vector. The step detection dead-reckoning may be used to determine a synthetic range, which may then be treated identically to the range estimates for the ultrasound transmitter units 2, 3, 4, 5 determined using TOF measurements.

[0213] FIG. 9 shows the steps carried out as part of such a synthetic range based position estimation process.

[0214] To determine the synthetic range estimate, an initial position estimate, p.sub.1, derived from the corrected state of the PDR algorithm, recorded at a time t.sub.1, is stored in the memory 208 of the mobile device 7 in step 800. At a second later time, t.sub.2 (e.g. the next time at which the IMU 209 updates), an updated position estimate p.sub.2 and corresponding covariance P.sub.p.sub.2 (corresponding to the most recent corrected state of the PDR algorithm) are determined, in step 802, by the controller 205 of the mobile device 7 using the PDR algorithm described above. An estimated range, r, between the initial position and the updated position, and an approximate range variance .sub.R.sup.2 are then determined, as will be described in the following, in step 804.

[0215] The initial position, p.sub.1, of the mobile device 7 is incorporated into the cost function as if it were the position of a synthetic transmitter (i.e. equivalent to another transmitter position r.sub.i in the ultrasound positioning cost function J.sub.1 above). A range measurement, r is also incorporated, as if it the range, R.sub.i, between the synthetic transmitter and the mobile device 7 in the positioning cost function J.sub.1. This value is calculated by the controller 205 as the distance between the synthetic transmitter position at time t.sub.1 and the estimated position at time t.sub.2, i.e. as r=p.sub.2p.sub.1.

[0216] The range variance, .sub.R.sup.2, can be used as an error term for weighting the contribution from the synthetic dead-reckoning transmitter in the augmented cost function. It is computed from the position covariance matrix at time t.sub.2, which is obtained by a first order Taylor expansion of the range calculation:

.sub.R.sup.2=HP.sub.p.sub.2H.sup.T

where H is the Jacobian defined by:

[00013]$H=\frac{1}{\stackrel{}{r}}\left({\stackrel{}{p}}_2-{\stackrel{}{p}}_1\right)$

[0217] The synthetic transmitter position, p.sub.1 and the synthetic range, r, are then directly incorporated into the cost function J.sub.1 given above (in step 806), as an additional transmitter position and an additional range between the transmitter device and the mobile device 7, with a corresponding weight, w.sub.p given by the range variance .sub.R.sup.2 and a scaling factor .

[0218] The resulting cost function, J.sub.1*, incorporating the synthetic range r and the weighted range variance .sub.R.sup.2 is therefore:

[00014]$J_1^{*}\left(r;r_i,R_i\right)=\frac{.Math.w_i\left({\left(r-r_i\right)}^2-R_i^2\right)+w_p\left({\left(r-{\stackrel{}{p}}_1\right)}^2-{\stackrel{}{r}}^2\right)}{.Math.{\mathrm{rel\_w}}_{i,p}}$

where rel_w.sub.i,p is an optional relative component representing the relative confidence in the reliability of the measurement from each transmitter unit 2-5, indexed by i, and supplemented by the synthetic transmitter, p.

[0219] This augmented cost function J.sub.1* can then be minimized over r, to determine an updated position estimate, r, taking into account the additional synthetic transmitter.

[0220] This approach allows a synthetic range term to be included in the cost function J.sub.1* directly, without introducing more complexity to the cost function.

[0221] Thus in some embodiments the PDR algorithm used by the controller 205 to calculate a synthetic transmitter position and a synthetic range, which are weighted by the range variance and added to an ultrasound positioning cost term in a cost function, J.sub.1*. The cost function J.sub.1* may then be minimized to determine an improved position estimate for the mobile device 7.

[0222] By solving a minimization problem involving this cost function, in step 808, using any suitable technique such as the Newton-Raphson method or a quasi-Newton method, the controller 205 can calculate an updated position estimate.

[0223] Including a synthetic range term in the positioning cost function J.sub.1* as described here can be done without relying on absolute heading. The PDR process may therefore, in some embodiments use only relative headings (i.e., determined relative to a default initial heading that may be set to zero).

[0224] However, although no heading term is included in the cost function J.sub.1*itself, changes in heading are nevertheless taken into account when calculating the positions estimates p.sub.1, p.sub.2 using the PDR algorithm. The augmented cost function J.sub.1* includes a term centred around the earlier dead-reckoning position estimate, p.sub.1, and hence the updated overall position estimate, r, is dependent on changes in heading as the user walks with the mobile device 7, determined from IMU measurements. By updating the augmented cost function J.sub.1* at, or closer to, the IMU update frequency, or at least more often than the rate of ultrasound transmissions, an updated position can be determined at a higher frequency, and with greater accuracy, than is possible using range estimates determined from TOF measurements alone.

[0225] In some variant embodiments, if the time between updates is small enough (e.g. less than the duration of one user step), then a reasonable assumption is that the user heading remains approximately constant between updates. In this case, the heading estimation may optionally be removed from the PDR dynamic model, and the synthetic range and the range variance to be added to the cost function may be reduced to terms involving only step length and total fractional step count, as well as their associated covariances:

r=L

.sub.R.sup.2=.sup.2P.sub.L+L.sup.2P.sub..sub.

[0226] While an initial position, p.sub.1, is still necessary as a simulated beacon position, the synthetic range, r, need no longer depend on changes in heading, but only on the magnitude of distance travelled by the mobile device 7 (i.e. step length multiplied by the value of the total accumulated steps taken). This can reduce the computational load on the controller 205.

[0227] Approach IIUsing Mahalanobis Distance

[0228] An alternative method, implemented in some embodiments, is to incorporate the PDR-derived position into the optimization routine by adding a second cost term, J.sub.2, that increases with distance from the latest position estimate determined by the PDR algorithm.

[0229] FIG. 10 outlines the steps of an exemplary implementation of this method.

[0230] In contrast to the first approach, this second approach does not rely on a synthetic transmitter position and a range, determined from two different position estimates determined from the PDR state vector, but instead relies on a single dead-reckoning position estimate and the dead-reckoning covariance matrix P.sub.p.sub.k.

[0231] To achieve this, the most recent position estimate p.sub.k=(x.sub.k, y.sub.k) and estimated 2-by-2 covariance matrix P.sub.p.sub.k are determined, in step 900, using the PDR algorithm as described above. The distance between the most recent position, determined using the PDR algorithm, p.sub.k, and the estimated position, r, that is to be determined through the cost optimization process (i.e. the same r used in the ultrasound positioning cost term J.sub.1), is then given by a second cost term J.sub.2 (established in step 902). This distance may, in some embodiments, be given by the Mahalanobis distance, i.e.:

J.sub.2(r;p.sub.k)={square root over ((rp.sub.k).sup.TP.sub.p.sub.k.sup.1(rp.sub.k))}

[0232] The cost term J.sub.2(r;p.sub.k) is then weighted, using a scaling factor , and added to the cost term J.sub.1 defined above, to establish, in step 904, a cost function J.sub.3 for determining the position of the mobile device 7, as follows:

[00015]$J_3\left(r;r_i,R_i,p_k\right)=\frac{.Math.w_i\left({\left(r-r_i\right)}^2-R_i^2\right)}{.Math.{\mathrm{rel\_w}}_i}+\sqrt{{\left(r-p_k\right)}^T{P}_{p_k}^{-1}\left(r-p_k\right)}$

[0233] Solving a minimization problem based on this cost function, in step 906, using any suitable technique such as the Newton-Raphson method or a quasi-Newton method, then enables the controller 205 to calculate an updated position estimate.

[0234] This method produces an improved heading estimate due to the added constraints on the optimized position. The heading correction is an efficient process based on the weighted means estimated by the ultrasound algorithm and PDR process, where the weights are determined by their respective variances.

[0235] However, in other embodiments of this second approach, alternative methods may be used to extend the ultrasound cost term. For instance, in some embodiments, the position estimate determined using the PDR algorithm may be converted into a matrix of distinct sigma points (as used in an unscented Kalman filter) each weighted by a respective position covariance, and a second cost term, J.sub.3, may be defined based on the weighted distance to each of the sigma points.

[0236] More generally, it will be appreciated by those skilled in the art that the invention has been illustrated by describing one or more specific embodiments thereof, but is not limited to these embodiments; many variations and modifications are possible, within the scope of the accompanying claims.