WO2021010832A1

WO2021010832A1 - A method, a system and a computer program product for estimating positions of a subject's feet and centre of mass relative to each other

Info

Publication number: WO2021010832A1
Application number: PCT/NL2020/050466
Authority: WO
Inventors: Mohamed Irfan MOHAMED REFAI; Petrus Hermanus Veltink; Bernhard Johannes Frederik VAN BEIJNUM; Jacob Hilbert BUURKE
Original assignee: Universiteit Twente
Priority date: 2019-07-15
Filing date: 2020-07-15
Publication date: 2021-01-21

Abstract

The invention relates to a method for estimating positions of a subject's feet and centre of mass relative to each other during gait. The method comprises a step of collecting measurement data from a first inertial measurement unit located at a first foot or first shank of the subject, from a second inertial measurement unit located at a second foot or second shank of the subject, and from a third inertial measurement unit located at a pelvis of the subject. Further, the method includes a step of evaluating relative positions of the first and second foot and the center of mass of the subject over time, using the measurement data of the first, second and third inertial measurement unit. The method may include a step of applying the assumption that a moment around a center of mass of the subject vanishes for determining an estimation of relative foot positions.

Description

Title: A method, a system and a computer program product for estimating positions of a subject’s feet and centre of mass relative to each other

The invention relates to a method for estimating positions of a subject’s feet and centre of mass relative to each other during gait.

Ambulatory estimation of gait measures is useful in understanding gait patterns in healthy subjects, and also recovery in people with gait impairment. One possible ambulatory method is to use Inertial

Measurement Units, IMUs. IMUs include 3D accelerometers and 3D gyroscopes, and are small and wearable. They can be used to estimate full body kinematics, and also kinetics, if a full body suit of IMUs is used. A full body suit of IMUs provides information about the individual body segments, and applying biomechanical constraints as described e.g. in the article entitled“Xsens MVN : Full 6DOF Human Motion Tracking Xsens MVN ; Full 6DOF Human Motion Tracking Using Miniature Inertial Sensors”, Tech, Rep. 3, 2009 by Roetenberg et al. improves estimation of kinematics of individual segments. Drift caused by strap down inertial navigation is prevented by using biomechanical constraints regarding the measured segments. However, when reducing the number of sensors used, the number of possible biomechanical constraints is reduced and the influence of drift increases. Another approach to ambulatory estimation of gait is the use of ForceShoes™ (reference: Weenk et al. 2015 and van Meulen et al. 2016).

This specially designed shoe consists of an IMU, and 3D force and moment sensor for each foot. An ultrasound sensor system can be added to the setup.

It is noted that patent publication US 2011/0152727 Al discloses a system for fall prevention for a user, comprising one or more sensors for attachment to respective portions of the body, each sensor being adapted to measure movement of the respective portion of the body and to translate the movement into a signal; and a controller adapted to receive the signal or signals and to determine a risk of falling by estimating a trajectory of the centre of mass of the user relative to a trajectory of a lower portion of the body from the signal or signals.

It is further noted that patent publication US 2008/0285805 Al discloses a system for capturing motion of a moving object via a plurality of motion sensor modules placed on various body segments. The sensor modules capture both 3D position and 3D orientation data relating to their respective body segments, thereby gathering motion data having six degrees of freedom with respect to a coordinate system not fixed to the body.

Traditionally, reliable drift reduced estimations of relative distances between each foot and the pelvis, where the center of mass is assumed to lie, are estimated by either using a full body suit such as depicted in US 2008/0285805 Al, which has knowledge of all individual body segments, or using sensors that measure relative foot distance, such as infrared or ultrasound sensors as in the ForceShoe™ (Weenk et al. 2015). These estimations are used for measuring several spatio-temporal gait parameters, including the centre of mass, CoM, base of support, BoS, and Margin of Stability, MoS, while walking. These parameters are required for estimating gait and balance quality in people.

It is an object of the present invention to provide an improved method for estimating positions of a subject’s feet and center of mass relative to each other during gait. Thereto, according to the invention, a method is provided for estimating positions of a subject’s feet and center of mass relative to each other during gait, comprising the steps of:

- collecting first measurement data from a first inertial measurement unit located at a first foot or first shank of the subject;

- collecting second measurement data from a second inertial measurement unit located at a second foot or second shank of the subject;

- collecting third measurement data of a center of mass from a third inertial measurement unit located at a pelvis of the subject, and - evaluating relative positions of the first and second foot and the center of mass of the subject over time, using the first, second and third measurement data.

By evaluating measurement data from IMUs at the feet and the pelvis of the subject, most preferably using only three IMUs, relative positions of the feet and the center of mass of the subject, as well as velocities thereof, can be estimated over time. This is usually performed using zero velocity update, ZVU, (Skog et al. 2010) for improving the estimations from a strapdown inertial navigation method. This provides change of positions and velocities of the feet over time during a step. There is eventual issue of drift between the feet and CoM, which can be reduced or eliminated using special biomechanical constraints. This is e.g. feasible by applying the theory of zero moment point, ZMP, which provides distance relations between said center of mass and the feet without the need for additional sensing. Generally, position and/or velocity estimations are in 3D. Further, based on the estimations, spatio-temporal gait parameters, and also balance metrics such as the Base of Support, BoS, and/or the Margin of Stability, MoS, can be measured. The relative positions may be expressed in a global reference frame that changes with the direction of gait.

Further, based on only the third measurement data, i.e. data from the inertial measurement unit at the pelvis of the subject, force vector components exerted at the center of mass of the subject can be evaluated.

By using a special cahbration method an initial estimate of orientation of the feet and pelvis can be improved. Optionally, an error extended Kalman filter, EEKF, estimates the orientation of the body during gait, which can be used for estimating the Ground Reaction Forces, GRF using just the third measurement data.

The positions of the feet, CoM, and Ground Reaction Forces, GRF, can be estimated in a so-called current global frame, i.e. a fixed frame of reference for the body. Different methods of estimating current global frame are feasible including a first case of defining the frame based on the pelvis orientation and the average change in direction over a cycle. A second case could be based on defining the current global frame based on the change in translation of the foot over the ground before and after a step. In both cases, the vertical axis of the current global frame is determined using inclination of gravity. The information for this current global frame can either be obtained from only the third measurement for the first case, or the first and second measurement data for the second case. The current global frame can also be referred to as the current step frame, and the terms are used interchangeably .

Based on the information described above, foot kinematics and kinetics can be estimated, e.g. using filtering techniques such as extended Kalman filters, for estimating foot and center of mass positions over time. The estimates can be improved based on inter segmental distances measured with application of the theory of zero moment point above other commonly used biomechanical constraints, including zero velocity update. Then, human spatial and temporal gait parameters, specifically, position and velocity of a subject’s feet and position and velocity of center of mass during gait, and balance metrics can be estimated, especially using only three inertial measurement units. Balance metrics stated here may include Base of Support, BoS, and Margin of Stability, MoS.

The invention also relates to a system, comprising a first, second and third inertial measurement unit for location at the feet and pelvis respectively, as well as a processing unit for receiving respective

measurement data from the inertial measurement units and evaluating relative positions of the first and second foot and the center of mass of the subject over time.

Further, the invention relates to a computer program product. A computer program product may comprise a set of computer executable instructions stored on a data carrier, such as a CD or a DVD. The set of computer executable instructions, which allow a programmable computer to carry out the method as defined above, may also be available for

downloading from a remote server, for example via the Internet, e.g. as an app.

Further advantageous embodiments according to the invention are described in the following claims.

It should be noted that the technical features described above or below may each on its own be embodied in a system or method, i.e. isolated from the context in which it is described, separate from other features, or in combination with only a number of the other features described in the context in which it is disclosed. Each of these features may further be combined with any other feature disclosed, in any combination.

The invention will now be further elucidated on the basis of a number of exemplary embodiments and an accompanying drawing. In the drawing:

Fig. 1 shows a flow chart of a method according to the invention;

Fig. 2 shows a system according to the invention carried by person P to be monitored;

Fig. 3 shows the person P shown in Fig. 2, and two right foot positions before and after a step, and the different sources for position estimates;

Fig. 4 shows trajectories of zero moment point and centre of pressure during a gait cycle;

Fig. 5 shows trajectories of the feet and CoM in a top-down view;

Fig. 6A shows a pelvis based current global frame;

Fig. 6B shows a foot based current global frame;

Fig. 7 shows the transformation processes to transform data from sensor frame to current global frames ;

Fig. 8 shows an EEKF filter used to track pelvis orientation, and measure the GRF ; Fig. 9 shows the GRF estimated using only three IMUs,;

Fig. 10 shows a brief overview of estimating the velocity of center of mass using the three IMUs;

Fig. 11 shows an example of the two sources for center of mass velocity: VCOM_hf and VCOMi_l, that are fused;

Fig. 12 shows an example of an estimated instantaneous velocity of center of mass (VCOM_est) compared with a reference;

Fig. 13 shows a complete overview of estimating feet and center of mass positions from the three measurement data, and applying the mentioned constraints and updates;

Fig. 14 shows a trajectory of CoM Height estimated from the pelvis IMU for a subject performing a walk and turn task, and compared with reference optical data;

Fig. 15 shows a step-wise comparison of foot and CoM positions tracked by the system described in this patent, with reference from optical data;

Fig. 16 shows comparison of step lengths measured by the system described in this patent, with reference from optical data;

Fig. 17 shows comparison of step widths measured by the system described in this patent, with reference from optical data;

Fig. 18 shows comparison of CoM widths measured by the system described in this patent, with reference from optical data;

Fig. 19 shows comparison of distributions of right and left step lengths for a normal walk task, measured by the system described in this patent, with reference from optical data, and

Fig. 20 shows comparison of distributions of right and left step lengths for the asymmetrical walk task measured by the system described in this patent, with reference from optical data. It is noted that the figures show merely preferred embodiments according to the invention. In the figures, the same reference numbers refer to equal or corresponding parts.

In evaluating the theory of Zero moment point, ZMP, it is assumed that for a stable gait pattern, the moments around the center of mass, CoM, are zero. The ZMP follows the center of pressure, CoP. Solving the ZMP assumption, relations between ZMP, CoM, and distances between the feet can be derived. This could be translated to relative distances between each foot and the CoM. This could be potential information that can reduce drift in position estimations. Therefore, the assumptions of ZMP are possible as biomechanical constraints that could be useful as a measurement update for a proposed sensor fusion filter that tracks the velocities and positions of the feet and center of mass. The assumptions of ZMP may also be used to estimate relative segmental distances and reduce drift using methods other than sensor fusion. Note that the ZMP is considered similar to a Centroidal Moment Pivot, CMP point for this apphcation, and the terms will be used interchangeably (Popovic et al. 2005)

Figure 1 shows a flow chart of a method according to the invention. The method is used for estimating positions of a subject’s feet and CoM relative to each other during gait. The method 100 comprises a step of collecting 110 first measurement data from a first inertial measurement unit located at a first foot of the subject, a step of collecting 120 second measurement data from a second inertial measurement unit located at a second foot of the subject, a step of collecting 130 third measurement data from a third inertial measurement unit located at a pelvis of the subject, and a step of evaluating 140 relative positions of the first and second foot and the center of mass of the subject over time, using the first, second and third measurement data.

The step of evaluating relative positions can be performed using dedicated hardware structures, such as FPGA and/or ASIC components. Otherwise, the method can at least partially be performed using a computer program product comprising instructions for causing a processor of a computer system to perform the above described steps. The step can in principle be performed on a single processor. However it is noted that at least a sub-step can be performed on a separate processor, e.g. a sub-step of collecting measurement data.

Further, a step of building a sensor fusion filter as seen in Fig. 13, for fusing information from strapdown inertial navigation including zero velocity update, ZVU, and relative foot positions, based on zero moment point, ZMP, to improve foot tracking as a function of time can be

implemented in the computer program product. Also, the computer readable code may include code causing a processor to perform the further step of displaying the information in a frame of reference defined with respect to the pelvis or the direction of gait as described above as a current global frame.

Generally, in strapdown inertial navigation, the position, including velocity, of a body of interest is estimated from integrating the accelerations measured on the body, after orienting it appropriately to an outside reference frame using information from the body’s angular velocity. The accelerations and angular velocities are measured using an accelerometer and (rate) gyroscope or angular velocity sensor, respectively, placed on the body. Optionally, the method may also use magnetometer information. The sensor that incorporates the accelerometer(s), and gyroscope(s) or angular velocity sensor(s), is called an inertial measurement unit or IMU. A combined sensor unit including an IMU and a magnetometer is called an inertial and magnetic measurement unit. Generally, when referring to an IMU, the IMU might be extended to incorporate information from a magnetometer that is separate or included in a combined sensor unit.

Figure 2 shows a system 1 according to the invention carried by person P to be monitored. The system includes three IMU’s, e.g. a first IMU 2a mounted at a first, right foot, a second IMU 2b mounted at a second, left foot, and a third IMU 2c mounted at the pelvis of the person P. CoM is denoted as a blue circle and its projection on the ground is GP. XCoM is a virtual point on the ground which encodes the position of CoM including its velocity. The centre of pressure, CoP, a yellow circle, is the point of contact of the body on the ground. The ground reaction force, GRF, acts from CoP and is directed towards the CoM as is shown as F. Tp is the trajectory of the CoM over time. The system 1 further includes an on-body processing unit 3 that is arranged for receiving measurement data from the respective IMUs 2, and for evaluating relative positions of the first and second foot and the center of mass CoM of the subject over time, based on said data. Optionally, the processing unit 3 can be arranged for generating on-body feedback signals for guiding / coaching the person P based on the measurements performed. In the shown embodiment, the processing system 3 is carried by the person P. Alternatively, the processing unit 3 is at least partly implemented as a separate stand-alone system, such as a stationary server, for performing at least a number of the above-mentioned steps, i.e. the step of evaluating the relative positions. Here, the system 1 includes merely three IMUs 2a-c.

In this respect it is noted that the first IMU 2a may be mounted to the first shank of the person P instead of the first foot of said person P. Similarly, the second IMU 2b may be mounted to the second shank of the person P instead of the second foot of said person P. Also, the third IMU 2c may be mounted anywhere at the pelvis of the person P or anywhere at a circumferential position around the hip.

When applying the above method, first, second and third measurement data are collected from the respective inertial measurement units IMUs, located at the feet and the pelvis of the subject to be evaluated. In an embodiment, the data is received wirelessly to the processing unit 3 located on or outside the body. All the three measurement data should be synchronized in time and stored on a storage system that may be

implemented in the processing unit 3 for further processing. This processing may include relative positions of the first and second foot and the center of mass of the subject evaluated over time, using the first, second and third measurement data.

In principle, the method may further include a step of evaluating the velocity of the center of mass of the subject over time, preferably at any moment of time, from the third measurement data from the third inertial measurement unit and may also be based on a position estimation of the center of mass as a function of time. Alternatively or additionally, the velocity of the center of mass can be estimated at any moment in time from the velocity estimates of the feet which have been corrected with

biomechanical constraints such as zero velocity update. The estimate of velocity of center of mass may be improved by fusing different sources of information, which may include the average velocity from the velocity estimates of the feet, and instantaneous velocity from strapdown integration of accelerations at center of mass, as shown in Fig. 10. Such improvement of the estimate of velocity of center of mass is e.g. illustrated in the third exemplary implementation.

Also, force vector components exerted at the center of mass of the subject can be evaluated, using the first, second and third measurement data. Similarly, force vector components exerted at the center of mass of the subject can be evaluated, using only e.g. the third measurement data. Then, referring to Fig. 2, a single IMU setup could be applied, the person P then carrying a single IMU 2c at the pelvis, without IMUs at the feet. An orientation of the inertial measurement unit can be estimated, preferably, using an error extended Kalman filter, as shown in Fig. 8.

Typically, position, velocity and/or force are evaluated in three- dimensional space, 3D. Figure 3 shows the center of mass CoM of the person P shown in Fig. 2, and two right foot positions. Specifically, a first position P1 of the right foot, before making a step, as well as a second position P2 of the right foot, after making the step, is shown. A trajectory T that is made during the step is also depicted, connecting the first and second positions P1, P2. Tp is the trajectory of the CoM during this step. The figure shows a vertical line L extending from the CoM vertically downwardly to a point GP at the ground. GP is the projection of the CoM on the ground. The second right foot position P2 is offset from the ground point GP with a distance OFF. The distance OFF is also defined as the distance between the contact foot and CoM in the horizontal or ground plane.

Generally, a walking pattern repetitively includes a first phase wherein a first foot is resting on a support structure supporting the first foot, and a second foot is moving, followed by a second phase wherein the second foot is resting on a support structure supporting the second foot, and a first foot is moving.

When assessing the first phase, the method may include the step of determining a first estimation of the second foot position relative to the first foot, after a step of the second foot, applying the ZMP assumption that a moment around a center of mass of the subject vanishes. Further, the method may then include a step of determining a second estimation of the second foot position relative to the first foot, after the step of the second foot, by applying strap-down navigation estimation of displacement of the second foot during the step of the second foot, combined with initial and final conditions that follow from zero velocity update, and an estimated position of the second foot relative to the first foot, before the step of the second foot. This results in a trajectory similar to T as shown in Fig. 3 for the second foot.

In this context it is noted that the system 1 may include a sternum sensor such as an accelerometer, gyroscope and/or angular velocity sensor for estimating a movement of the upper body of the person P, thereby further improving estimates made using the ZMP assumption.

Also, the method may include a step of improving an estimation of the second foot position relative to the first foot, after the step of the second foot, by combining the first and the second estimation of the second foot relative to the first foot, preferably taking into account uncertainties of the first and second estimation, preferably using a specifically designed sensor fusion filter.

Similarly, when assessing the second phase, the method may include the steps of performing a first estimation, a second estimation and an improved estimation of the first foot position relative to the second foot, after the step of the first foot, by performing the steps defined above wherein the first foot and the second foot have interchanged.

As indicated above, the method may also comprise of first estimate of determining the velocity of the center of mass of the subject, e.g. based on the third measurement data from the third inertial measurement unit located at a pelvis of the subject, preferably including a strap-down navigation estimation. A second estimate of the velocity of the center of mass may be derived as an average velocities of the feet while moving. A fusion of first and second estimate may advantageously be used to get an improved estimate of velocity of center of mass. Further, by fusing the velocities for obtaining an improved estimate may highly improve the estimates for positions later on.

Also, the position of the center of mass of the subject as a function of time relative to the first foot when the second foot is moving and relative to the second foot when the first foot is moving may be estimated by applying the ZMP assumption that a moment around a center of mass of the subject vanishes. This leads to the estimate OFF in Fig. 3. The estimation of the position of the center of mass may be improved by including the third measurement data from the third inertial measurement unit located at a pelvis of the subject, preferably including a strap-down navigation

estimation.

The method may also comprise a step of estimating the position of the center of mass of the subject during double stance based on the third measurement data from the third inertial measurement unit located at a pelvis of the subject, preferably including a strap-down navigation

estimation. This is shown as Tp in Fig. 3.

Advantageously, the position estimation of the center of mass can be improved by fusing all information about center of mass movements available over time, preferably using a specifically designed sensor fusion filter

Preferably, a location, in particular the height of the center of mass of the subject, can be estimated based on measurement data from an inertial measurement unit located at the pelvis of the subject, preferably by applying strap-down inertial estimation and using biomechanical

constraints regarding the height of the center of mass.

Foot kinematics and kinetics can be determined relative to a frame of reference that is stationary relative to the pelvis or stationary relative a direction of gait of the subject. The first case of defining a reference frame or current global frame relative to the pelvis is performed by utilizing an information about the average heading of the pelvis as the forward direction. The second case of defining a reference frame relative to direction of gait may include a step of estimating foot kinematics and kinetics over time, based on measurement data from a first inertial measurement unit located at a first foot, and measurement data from a second inertial measurement unit located at a second foot. The change in foot positions before and after a step can be used as the forward direction to estimate a new reference frame, or current global frame or current step frame. All foot and center of mass kinematics and kinetics can be expressed in the current global frame that is continuously changing. In both first and second cases of defining this current global frame, the vertical provides one of the other axes of the current global frame.

In principle, balance metrics can then be estimated, based on the kinematics of the foot and the center of mass of the subject, e.g. by

performing the steps of:

- estimating a base of support defined by boundaries of the feet while in contact with the ground during walking;

- estimating the base of support defined by projected boundaries of the feet if either one is not in contact with the ground during walking;

- estimating an extrapolated center of mass, which includes the information of velocity and direction of walking, and

- estimating a margin of stability, which is the distance between the extrapolated center of mass and the base of support.

The step of estimating foot kinematics and kinetics over time may include applying a filter method including an extended Kalman filter to measurement data.

Foot kinematics can be updated over time using the first distance between the position of the first foot and the center of mass of the subject, and/or the second distance between the position of the second foot and the center of mass of the subject, which are derived from the ZMP. The foot kinematics can also be updated over time using other biomechanical constraints.

It is noted that the step of evaluating the velocity of the center of mass of the subject over time can be performed in combination with the method as defined in claim 1, but also more generally, with a method for estimating positions of a subject’s feet and centre of mass relative to each other during gait, comprising the steps of collecting first measurement data from a first inertial measurement unit located at a first foot or first shank of the subject, collecting second measurement data from a second inertial measurement unit located at a second foot or second shank of the subject, collecting third measurement data from a third inertial measurement unit located at a pelvis of the subject, and evaluating relative positions of the first and second foot and the center of mass of the subject over time, using the first, second and third measurement data.

It is also noted that the step of evaluating force vector components exerted at the center of mass of the subject, using the first, second and third measurement data, or using only the third measurement data can be performed in combination with the method as defined in claim 1, but also more generally, with a method for estimating positions of a subject’s feet and centre of mass relative to each other during gait, comprising the steps of collecting first measurement data from a first inertial measurement unit located at a first foot or first shank of the subject, collecting second measurement data from a second inertial measurement unit located at a second foot or second shank of the subject, collecting third measurement data from a third inertial measurement unit located at a pelvis of the subject, and evaluating relative positions of the first and second foot and the center of mass of the subject over time, using the first, second and third measurement data.

Below, four exemplary implementations are described in more detail. The first exemplary implementation, as published at the 41^st Annual International Conference of the IEEE Engineering in Medicine and Biology Society (EMBC) as“Portable Gait Lab: Zero Moment Point for Minimal Sensing of Gait”, explains the theory of ZMP and validates it using high resolution optical tracking systems. The example also provides an idea of the sources of error regarding the assumptions of the theory of ZMP. The second exemplary implementation, as pubhshed in an article entitled “Portable Gait Lab: Estimating 3D GRF Using a Pelvis IMU in a Foot IMU Defined Frame” in IEEE Transactions on Neural Systems and

Rehabihtation Engineering, 28, pages 1308-1316, describes the feasibility of deriving GRF from the minimal sensor setup. The third exemplary implementation shows the estimate of CoM velocity by fusing information from foot velocity and CoM accelerations. The fourth exemplary

implementation shows the design of the the sensor fusion filter that would enable estimating the positions of a subject’s feet and centre of mass by fusing information from different sources, including ZMP.

First exemplary implementation

First, the assumptions of zero moment point are to be validated. As shown in Herr et al. 2008 (Angular moment in human walking), the ZMP (similar to the Centroidal Moment Pivot, CMP, in Herr et al. 2008), is a point where the ground reaction force would have to act to keep the horizontal component of the whole-body angular moment constant. When these moments around the CoM is zero, the ZMP coincides with the CoP. Therefore, first, the differences between ZMP and the CoP for walking in two conditions is studied: in a normal and a casted walking condition.

Following this, a relation derived from the ZMP assumption that could provide relative distances between each foot and CoM is explored. Further, an IMU example is illustrated to describe steps to implement the ZMP in IMU based sensor fusion approaches.

Seven healthy female subjects were asked to walk on the GRAIL (Motekforce Link, The Netherlands). The GRAIL consists of a split belt treadmill with force plates and ten VICON motion capture cameras to collect gait biomechanics. The setup measured the 3D ground reaction forces and also the 3D kinematics of the body positions. The subjects’ average age was 22.9±1.4 years, height was 1.78±0.06m, and weight was 73.4±5.4kg . The subjects were asked to walk for 5 minutes on the treadmill at 1.2 m/s. After this, a plaster technician casted the right foot of the subject and they were asked to walk again for 5 minutes at the same speed on the treadmill. Casting was done to induce asymmetry in gait.

It is known to use a complementary filter algorithm to estimate low and high frequency components of CoM from the 3D forces and moments measured from the ForceShoe™ . Here, we apply the algorithms to measurements from the GRAIL. The first stage estimates CoM from both foot kinetic and kinematic information by low pass filtering the CoP to estimate the position of CoM. The total body CoP is estimated as follows from the force measurements on the GRAIL (Eq. 1).

All variables in (Eq. 1) are a function of time. Equation (1) is depicted in the global frame with Z axis positive upwards along the vertical, and X positive along the walking direction. Here, subscript r and l stand for the right and left foot respectively, and subscript ax corresponds to either X or Y axes. F_z refers to the force in the Z axis. The CoP_ax is then low pass filtered at 0.4 Hz to obtain the CoM_ax;LF .

The second algorithm estimates CoM from kinetic information alone by double integration of the net forces based on Newton’s second law. The acceleration of the body mass m_body at the CoM is given as follows (Eq. 2).

Here, F is the net force vector acting on the body, and g is the gravitational acceleration. The change in CoM position over time was derived from integrating the a_CoM twice. This results in x_CoM;int which was high pass filtered with a cut off at 0.4 Hz to obtain CoM_ax;HF . This is the same cut off as that of CoM_ax;LF’s low pass filter. The CoM_ax;LF and CoM_ax;HF are fused using a complementary filter to obtain the trajectory of CoM.

The theory of ZMP assumes that the moment around the CoM is zero for a stable walking trajectory. Here, ZMP is a point on the ground, such that the cross product of a vector joining the CoM and ZMP and the ground reaction force vector (F) is zero. This gives us the following equations (Eq. 3, 4, 5):

r x F = 0

In (Eq. 4), ZMP_z is zero as it lies on the floor, and r_z can be estimated as height of pelvis from the GRAIL system.

Therefore, (Eq. 5) provides ZMP positions in X and Y axes. This is then compared with the CoP estimated from the treadmill force plates using (Eq. 1). The difference between ZMP and CoM gives the value of OFF.

The relation between ZMP and CoM as shown in (Eq. 5) can be utilized as additional information about relative distance between the feet and CoM. Therefore, they can be used as measurement updates for a sensor fusion filter, if the other variables are known. For example, during swing phase of the left foot, the CoP of the body will lie under the right foot.

Although, in an IMU based approach it is not straight forward to measure CoP of each foot, the foot positions can be estimated using an IMU on each foot, and CoM can be tracked using a pelvis IMU. We can provide an estimate for the right foot during left foot swing phase as (Eq. 6):

Here, pos_ax,r is position of the right foot, and subscript sl denotes instances of left foot swing phase. We have assumed that the differences in ZMP and foot position is trivial. Subsequently, we can derive an estimate for the left foot during the swing phase of the right foot (sr), Eq. 7.

We then compare pos_ax,r and pos_ax,l with the true foot positions at the respective instances (si and sr).

We compare the ZMP_ax measured in (Eqn. 5) with the CoP measured from the treadmill force plates in GRAIL. Then, we compare the error as a percentage of subject’s foot length ^{( b% )} with Herr et al. 2008. Following this, we compare pos_ax,r and pos_ax,l from (Eqn. 6) and (Eqn. 7) with foot positions measured by VICON in GRAIL. Following this, we plot the trajectories of the interesting parameters in an example with IMUs.

Fig. 4 shows the normalised trajectories of ZMP and CoP in the X and Y axes for both conditions: normal and casted. The graph shown is the normalised gait cycle averaged over all subjects. The shaded regions show the standard deviation of the trajectories. The first column denotes the normal condition and the second column denotes the casted condition. Each row corresponds to one axis in the global frame. The mean absolute RMS of the differences between the ZMP and CoP over the complete cycle is shown in table I for both conditions. In table II, we compare the mean RMS of the distance between the ZMP and CoP across the gait cycle normahsed by foot length ^{( b% )} with that of Herr et al. 2008. Further, in table III, the mean RMS of the differences between the pos_ax,r and pos_ax,l and respective foot positions from GRAIL is shown for normal and casted walking conditions.

Figure 5 shows different trajectories in a top-down view. The dark green dotted line and blue dotted lines are the left and right foot trajectories respectively, estimated from the algorithm of Weenk et al. 2015 from the foot IMUs. The red dotted line is the reference centre of mass (CoM) estimated by the ForceShoe^TM The solid yellow lines are possible estimates for left foot during right swing phase, and solid green lines are possible estimates for the right foot during left swing phase. These are estimated using the CoM and zero moment point assumptions. The triangles denote the starting position of the two feet.

In practice, (Eqn. 3) is not valid. The moments around the CoM oscillate around zero. Upper body angular rotations also cause moments around the CoM. This is a missing component in (Eqn. 3). However, here we look at how the ZMP and CoP agree, during straight walking, where the moments around the CoM may be really small.

Fig. 4 shows that there is close overlap between the trajectories of ZMP and CoP for the normalised gait cycle. The gait cycle begins with right heel strike and we can see the transition of the CoP from left to the right foot. The CoP falls completely under the right foot around 15% of the gait cycle. Following the left swing phase, we notice the left heel strike around 50% of the gait cycle, as the CoP starts to move towards the left foot. The trajectory continues to the next right heel strike which is the end of the gait cycle. Looking closer, we observe that the standard deviation of the trajectories (both ZMP and CoP) are smaller during the transition from one foot to the other. In both normal and casted conditions, the trajectory of ZMP is closer to the CoP in the Y axis during these transition (double stance) phases, when compared to the swing phases. This could suggest that the moments around the CoM are smaller during double stance phase, thereby showing lower differences in the ZMP and CoP trajectories during these instances. TABLE I

MEAN RMS OF THE DIFFERENCES B ETWEEN ZERO MOMENT POINT AND

CENTRE OF PRESSURE OVER A GAIT CYCLE

TABLE II

MEAN DISTA NCE BETWEEN ZERO MOMENT POI NT AND CENTRE OF

PRESSURE ACROSS GAIT CYCLE NORMALISED BY FOOT LENGTH (b%) |

j

TABLE III

MEAN RMS OF D I FFERENCE B ETWEEN pos_ax AND REFERENCE FOOT

POS ITIONS AT RESPECTI VE INSTANCES . {* : p < 0,05)

Table I shows these differences as mean RMS of the differences between ZMP and CoP over the whole gait cycle. It is seen that the casting increases the error margins of the differences. Casting could induce asymmetry, causing increased rotation of the upper body to compensate for the change in walking pattern, and therefore, we see the differences in Table I. Table II shows the mean distance between ZMP and CoP across the gait cycle, normalised by foot length of the subjects.

Table III shows the differences between foot positions estimated from CoM using ZMP assumptions and true foot positions from GRAIL system. The larger errors can be explained by the fact that we are actually comparing ZMP estimates from CoM with foot positions. We assume ZMP to lie close to CoP for each foot, and further, assume differences between CoP and foot positions to be trivial. Therefore, the table III shows the error margins associated with these assumptions. Looking closer, the table shows that the error margins are around 9.5 cm in X axis, and about 1.6 cm in the Y axis, for the normal walking condition. They show larger deviations in case of casting. These margins give us an idea of the feasibility of using ZMP based assumptions for estimating foot positions from CoM.

The minimal sensing system could consist of three IMUs; one on each foot, and one at the pelvis as seen in Fig. 2. The foot IMUs could track the movement of the feet in 3D. Measurement updates such as zero velocity update will minimize the drift in the X and Z directions. The CoM can be tracked using a pelvis IMU. Fig. 5 shows the movement of the feet and CoM for an exemplary walking measurement where the subject walks forward for about 10 metres. The ZMP assumptions could be used during left swing to estimate relative position of CoM relative to right foot. Additionally, during right swing phase, we can estimate CoM relative to left foot. If we fuse all this information, we can estimate the relative positions between the two feet during subsequent stance phases. This removes the need for full body sensing, or an inter-foot distance sensor.

Eqn. 5 also requires knowledge of the height of the CoM and forces in 3D. CoM height can be measured by the pelvis IMU, with appropriate measurement updates. Estimating forces in 3D could be solved by the descriptions in the second exemplary implementation. If we assume that the body is only in contact with the ground, then the accelerations of the pelvis could be similar to the accelerations at the CoM. This is simply the specific ground reaction forces in 3D.

The current method assumes that the CoM position is used as a reference, and the estimates of the two feet could be corrected based on (Eqn. 6) and (Eqn. 7). An alternative method is to assume the right foot to be a reference point and then estimate the CoM, and subsequently, left foot position.

The errors in table III are majorly present as we compare ZMP with reference foot positions. Therefore, a possible solution could be to measure CoP during walking, as they show lower errors with ZMP, as can be seen in table I. In an ambulatory sensing setup, these errors could be solved by using a pressure insole to measure CoP providing more accurate relative distances between CoM and either foot.

This study shows possible applications of using ZMP assumptions to reduce the lateral drift during minimal sensing of gait using IMUs. The next step is indeed to build a sensor fusion algorithm that can implement these updates iteratively. It is advised that the assumptions are studied for different walking patterns.

Second exemplary implementation

The ZMP assumption mentioned earlier requires the knowledge of whole body 3D GRF. The main goal of the second exemplary implementation is thus to evaluate the feasibility of estimating whole body 3D GRF (GRF) from a single pelvis IMU, for different walking patterns seen in daily life. The GRF can be estimated using Newton’s law, given the mass of the body and accelerations measured at the centre of mass (CoM), assuming the body is only in contact with the ground. A special calibration procedure, and an Error Extended Kalman Filter (EEKF) are used to estimate body orientation, in order to determine the 3D components of the GRF. The other IMUs of the three IMU setup are placed on each foot. Using these, a changing frame of reference, that depends on the direction of step being made, is employed to express the GRF. This is opposed to a commonly used fixed global frame of reference. A reference frame that is changing with the moving body is a better option than an arbitrary fixed global frame, to express kinetics of the moving body as it will always be relative to the direction of gait. The GRF is then compared with that measured by the reference system, which is the ForceShoe^TM and also compared with results from existing literature.

GRF are usually expressed in a fixed global frame of reference.

This could be attributed to the fixed nature of force plates. As shown in Fig.

6A and 6B, the global frame ^g has a predefined and fixed frame

throughout the measurement. However, ambulatory setups allow us to define frames that are associated with the moving or turning body. For instance, pressure profiles, and centre of pressure patterns are expressed in foot frames. Similarly, here we can define reference frames that are attached to the moving body. There are two possible options. One option is to have a frame of reference attached to the pelvis as shown in Fig 6A. This would result in a frame defined along the anterio-posterior, and vertical axes of the body. Alternatively, we could define reference frames based on the direction of steps being made. The second method is easier to define for both the foot and pelvis segments, and is therefore preferred for this study.

Figure 6A shows a pelvis based current global frame, while Figure 6B shows a step based current global frame. This current global frame is interchangeably referred to as the current step frame. Here, the first position P1 of the right foot before a step k is shown, as well as the second position P2 of the right foot, corresponding to the positions shown in Fig. 3. Similarly, the trajectory T that is made during the step k is depicted.

Further, a trajectory Tp of the pelvis or CoM is shown. The changing reference frame or current global frame

defined graphically in Fig. 6B depends on the direction of the step, and thus, changes for each step. Here, we use it interchangeably as current step frame . fine the X axis of this frame as positive in the forward direction, defined by the line between the beginning and end of a step. The Z axis is positive upwards along the vertical. This provides the

for the current step, and is redefined for each step

The sensors on the IMU, however, measure in their own frame of reference, ^s . This has to be transformed to the

per step. Figure 7 shows the transformation matrices required to transform from the sensor frame measured by the IMU to the special current step frame. Note that the frame ^{Y cg} is interchangeable with ^Ψcs First, each sensor was transformed to their respective segment ( seg ) frames using a simple calibration method.

The segments of interest in this study are the pelvis (p ), left foot (fl), and the right foot (fr ). Then, the orientation of the segments have to be expressed in

^{Yc (k )} of a given step k. In order to estimate this, the change in foot position over the step k has to be known. Therefore, first, the change in orientation of the segments are estimated in a reference frame

that was built in the previous step k-1. EEKFs are used to estimate this change in

orientation during this step. Subsequently, the change in position of the moving foot can be estimated at the end of the step k. Once the initial and final positions are known, the rotation between the reference frames of the previous k - 1 and current k steps can be estimated as

. This is used to transform to the ^Ycs(k) . This can be redefined for each step. In short, four frames of reference were used in this study: sensor frame ( ^Ys ) for each sensor, segment frames (pelvis, right foot and left foot), reference frame defined by the previous step (^{Ycs k - 1 )}), and a current global frame defined by each step

.The notations used in this study are tabulated in Table IV.

The 3D accelerometer in the IMU provides the following signal in the sensor frame:

where a is the hnear acceleration of the body, g is gravity, and eA is Gaussian white noise. The 3D rate gyroscope present in the IMU

measures angular velocity in the sensor frame, where _w is the angular velocity, b is a slowly varying offset, and

is the Gaussian noise.

In order to estimate the GRF, the accelerations at the CoM has to be known. Here, we assume that the pelvis moves with the CoM, and that all mass is concentrated at this point. Additionally, the feet are the only contact with the external world, and no additional load is carried by the body. Therefore, the pelvis IMU was used to estimate the accelerations at the CoM, and eventually the GRF. The foot IMUs were used to estimate the direction of the steps being made, to define the

.

Table IV: Notation used, shown for an arbitrary vector a

Foot Contact and Step Detection: The method of Skog et al. 2010 was used to estimate the foot contact instances for the two feet. Step detection is important in estimating the current global frame or current step frame. As the IMUs are synchronized in time, the double stance instances can be estimated. The time instance for a step was defined to start from the midpoint of a double stance event, and ends at the midpoint of the next double stance.

Orientation in the different reference frames: The transformation between different frames as seen in Fig. 7 are explained in the following sections. The static calibrations used to transform the sensor data to the respective segment frame are first described. Following this, the structure of the fusion filter, used to estimate changes in orientation of the segment in the previous global frame is described. This is finally followed by estimating and translating all information to the current step frame.

First,

has to be estimated using the mounting frame calibration techniques described by Bonnet et al. 2009. The inclination estimate was estimated from (Equation below) during an initial standing still phase, during which the 3D accelerometer is expected to measure only gravity. The Y axis of the pelvis can be estimated by asking the subject to bend forward. Principal component analysis was applied to the gyroscope output to find the axis measuring largest angular rotation. The X axis of the pelvis can be estimated using right hand thumb rule. The orientation was then estimated using the following equations.

On the other hand, for the feet, the X axis was taken to lie along the direction of velocity vector for the first step being made. The Y axis was estimated using right hand thumb rule, and finally the rotation matrix was estimated. Following this, the change in orientation during movement has to be estimated. For a given step k, the

has been defined. Note that the can only be defined at the end of the current step. Therefore, the change in orientation was first expressed in

using an EEKF. The

EEKF was used to track the

, i.e., the orientation of the segment with respect to

for a given instance i. Here, i denotes the samples including the start and end of the current step k. The EEKF filter used for the pelvis orientation is shown in Fig. 8 and was based on Weenk et al. 2015, and Kortier et al. 2014. Similar EEKFs were built for the other segments too. The states included in the state vector (x) of the EEKF were orientation error and gyroscope bias, and its covariance matrix was denoted as P. The advantage of using an EEKF for estimating orientation errors is that the inertial processes can be considered linear, if the errors are assumed to be small. The state vector can be written as:

The initial orientation error

was assumed to be zero, and the initial gyroscope bias

was estimated from gyroscope output during the initial standing still phase, as the angular velocity should be zero. An initial estimate of

has to be estimated. For this, the direction of the first step being made was measured, by estimating the direction of the velocity vector of the foot that moved first. This gave the direction of the X axis, and the Z axis was taken to be the vertical. After estimating the Y axis using cross product, the

was estimated using earlier equations. We started by making models for each of the states. First, the gyroscope bias error was modelled as a first-order Markov process :

Here,

is white Gaussian noise associated with the process.

Thus, the gyroscope bias was predicted as :

The gyroscope bias error can then be written as :

The three preceding equations gives us (Eqn. 8):

Orientation can be estimated from orientation error as (Eqn. 9):

Equation (9) is valid when orientation errors are assumed to be small.

For any

. Furthermore, we can derive orientation from angular velocity using (Eqn. 10):

The derivative of orientation error and its discretized form is (Eqn 11 and

12)

The Kalman filter prediction equation is given as (Eqn 13 and 14)

The covariance matrix is predicted using (Eqn 15)

where, Q is the process noise covariance matrix.

The measurement updates were apphed to the EEKF using the standard equations. The measurement used at any instance i to update the state vector is given in (Eqn. 16).

is the noise associated with this measurement. When H is known, the Kalman gain is estimated using (Eqn. 17). Then, the state matrix and the error covariance matrix were updated with (Eqn. 18) and (Eqn. 19) respectively, see Eqn 16-19 below.

It was assumed that on average, the accelerations at the pelvis measure inclination due to gravity. This can be true during forward walking or when making changes in walking direction. The orientation error was corrected based on the current prediction and expected inclination for the vertical axis. An estimate of the accelerometer output can be estimated using the estimate of the orientation matrix as Eq (20)

Then, (Eqn. 20) was applied to the equation to estimate inclination at instance i. The difference between measured and estimated inclination, also referred to as

was then used to update the orientation error as shown below:

The orientation estimate and gyroscope bias were updated using the resulting state vector.

Using the estimates of the pelvis orientation

in each iteration, the accelerations were estimated from (Eqn. 20). The GRF can be estimated from the product of acceleration at the pelvis and body mass.

The GRF estimated so far has been expressed in the previous global frame. This has to be transformed to

for the current step k. As the change in orientations have been estimated for the current step, the change in foot positions can be estimated. For this, a separate extended kalman filter (EKF) was developed. The EKF estimates the velocity and position of each foot. Zero velocity and zero height updates are used to improve these estimations. At the beginning of each step, the foot position and velocity are reinitialised as we are only interested in the change in position for a given step. The EKF then tracks the change in position and velocity for the current step. For example, in Fig. 6B, the right foot (shaded green) changes direction while making the step k. The change in position in the XY plane was measured between the start and end of this step. The forward direction, shown by the dotted hne is the X axis, and the vertical is the Z axis. This gives us the new

for the step k. This was redefined for every new step being made, and then we transform the

to the new .

The

was found to have noisy peaks around foot contact.

This was removed by identifying local maxima or minima around foot contact and using a polynomial function to smooth the peaks. Following this, the signal was low pass filtered using a zero-phase Butterworth filter of order 4. However, it was observed that the cut off frequency of the filter influenced the resulting errors, and therefore the cut off frequency was adjusted to evaluate the optimum value.

A three IMU setup was used as seen in Fig. 2. One Xsens^TM IMU was mounted at the pelvic region on the sacrum using a strap. The IMU is placed such that it was at the midway point between the line connecting the left and right posterior superior iliac spine. One IMU was placed on each foot at the midfoot region. The MT Manager software was used to read the data from the IMU wirelessly, which was sampled at 100 Hz.

The ForceShoe™ was used as the reference system. The data from ForceShoe™ were sent wirelessly to a PC, sampled at 100 Hz. It was then low pass filtered twice at 10Hz using second order Butterworth filter to ensure zero-phase lag. The GRF measured by ForceShoe™ on each foot is summed to obtain the GRF, which is then transformed to the

as defined above.

Seven healthy subjects were recruited for the study. The inclusion criteria included subjects with no history of impaired gait or leg injury. The subjects began by standing still for a few seconds, following which they were asked to bend the trunk forward thrice. This was used to calibrate the sensor to segment orientation for the pelvis sensor. The subjects were then asked to perform different walking tasks similar to daily life. They were instructed to begin with their feet placed parallel. Once the researcher gave the start sign, the subject walked along a given path. The time taken between start and stop of the walking was measured using a stopwatch.

This activity was repeated six times. The walking tasks are performed in the following different scenarios:

1) Normal Walk (NW): During this task, the subject was asked to walk at his preferred walking speed for 10 m.

2) L Walk (LW): During this task, the subject was asked to walk for 15 m and then turn right at a right angle and walk for another 10m.

3) Slow Walk (SW): During this task, the subject was asked to walk at a slower pace. They were guided by the use of a metronome beating at 50 beats per minute. Each beat corresponded to a heel strike. This frequency was used so that the subjects walked slower than 0.6 m/s.

4) Walk and Turn (WT): During this task, the subject was asked to walk for 10 m and then turn and walk back to start position.

5) Slalom Walk (SIW): During this task, the subject was asked to walk in a slalom pattern. Pylons were placed on the floor for this purpose. Asymmetric Walk (AW): During this task, the subject was asked to walk in an asymmetric manner. The instruction given was to induce a stiff left knee and abduct the hip as much as possible on the right side.

The GRF estimated will be used to estimate relative distances between either foot and CoM as given by (Eq. 21):

In the above equation, subscript ax stands for the axis of interest and is either X or Y. Note that (Eqn. 21) is true only when the foot of interest is the only point of contact with the ground. Here, we see that the ratio of the forces in either X or Y axis with the Z axis is required to estimate the position of the foot of interest. Thus, the root mean square of the differences ( rmsd ) between the ratios calculated from

and

was estimated. The rmsd was then expressed as percentage of error as rRatxz and rRaty z, by normalising it with the range of the ratio measured from the reference

_

Some trials had to be excluded from the analysis due to issues with the sensor systems. However, it was made sure that each subject had at least three walking trials per task. Fig 9., wherein an example of the estimated

in the changing reference frame is seen. is shown in blue, and the reference

is in red. Each row in the figure corresponds to an axis in the current step frame. In this, the subject performs the LW task, where a turn to the right side is made at 35 seconds. This moment is shaded light red.

In Table V we compare the estimated post-processed instantaneous

_and reference

, for all walking tasks. The table shows the root mean square of differences (RMS), correlations (CORR), and also difference in 2D GRF vector angle in the XY plane

between the and

. The percentage error of the rmse ( rRatxz , and rRaty z ) calculated from

and

is tabulated in Table VI.

The study aims to estimate GRF using a minimal IMU only setup. This has a lot of potential in ambulatory monitoring. Fig. 9 shows an example of the

and

for the LW task. Here, the subject was asked to make a right turn, and the moment this occurs is denoted by a shaded rectangle. We see that after the turn, the GRF look similar to their profiles in their respective axes just before the turn.

However, if a fixed global frame, denoted in Fig. 6B was used, then after the turn, the F_x would show a profile similar to that of F_Y before the turn, and vice-versa. Instead, as the frame used here is defined by the direction of steps being made, the profiles of F_x and F_Y remain unchanged. This representation of GRF represents the biomechanics of the body irrespective of the change in walking direction. However, minor deviations can be noted at the moment of the turn, particularly in F_Y , possibly due to a larger deviation in this axis when making the right turn.

In this study, an EEKF was tuned to be able to resolve the accelerations measured at the pelvis into 3D components of the GRF within a specific reference frame. This is a tricky task as the contribution due to gravity is quite large as compared to the shear GRF. This resulted in high frequency noise creeping in the X and Y axis.

TABLE V: Comparing the Estimated and Reference GRF values: Root Mean Square of the Differences (RMS), Correlations (CORR), and Difference in 2D GRF Vector Angle in the XY Plane ( qd ).

TABLE VI: Root Mean Square of the Differences Between the Estimated and Reference GRF of the Ratio of GRF in either X or Y with the Z Axis, expressed as Percent Error of the Range of the Reference GRF

Table V shows that the RMS between the and

is low across the different axes. Additionally, CORR shows good agreement, especially in the X and Z axis. The Y axis shows lower correlations, which seems to worsen for the AW task. The AW task was meant to simulate gait impaired walking and is not a standardized test. The subjects were given instructions about how to walk asymmetrically, but each of them chose a unique pattern. Additionally,

shows quite some differences in the XY plane, which could be attributed to the low correlations in the Y axis. A possible solution to improve the estimations in the Y axis, is to further analyse the

biomechanics of the CoM for additional constraints or updates.

The RMS error margins could also be attributed to disagreements in transforming the measurement and reference system to the changing step frame.

In this example, we show that GRF can be estimated with good confidence. The ZMP equation further requires a ratio between two axes, and Table VI shows that on average there is a 13% error while measuring this ratio across walking tasks. This has to be acknowledged while solving (Eqn. 21). Looking closely, the AW task reflects the larger errors as seen earlier. Further, if we are able to measure the CoM kinematics, relative foot distance can be estimated.

Finally, it has to be noted that although foot IMUs were used to estimate the current step frames per step, they are not necessary for estimating the GRF during gait. It is possible to resolve the 3D components of GRF using only the pelvis IMU. As a measurement update, if an estimate of the forward CoM movement is known, it could be used to express the GRF in the pelvis based reference system, as also described above in Fig. 6A.

If this setup is to be used in people with impaired gait, a separate validation study has to be performed on the population of interest. The AW task performed in this study may not hold true for gait patterns exhibited by gait impaired populations. The same argument holds for rapid walking or running scenarios.

These results show that it is feasible to monitor GRF in an ambulatory manner using simply an IMU on the pelvis. Foot IMUs are used to provide a reference frame that allows us to measure the GRF with respect to the moving and turning body. This study is the first step in developing a minimised and portable gait lab.

Third exemplary implementation The third exemplary implementation shows the estimation of instantaneous of CoM velocity by fusing information from foot velocity and CoM accelerations.

Estimation of Centre of Mass (CoM) velocity (VCOM) has several practical applications, including measuring gait parameters, and balance measures such as Extrapolated Centre of Mass (XCOM). Ambulatory estimations are commonly done using Inertial Measurement units (IMUs). However, studies have either measured the magnitude of VCOM, or apply machine learning techniques which requires an additional training step.

Here, a setup of three IMUs is described for estimating the VCOM; one IMU at the pelvis, and one on each feet. High pass filtered strapdown integration of pelvic accelerations can provide an estimate of cychcal changes in VCOM. Additionally, on average foot velocities encode

information about the VCOM. These two sources of VCOM can be fused using a complementary filter method, resulting in drift free instantaneous estimates of VCOM. Here we assume that the pelvis IMU accelerations are similar to the CoM accelerations (ACOM).

The following sections describe the methods used to obtain the instantaneous 3D VCOM in a special current step frame (

) based on the principle described, and explains the performance of the method for variable gait.

First, the method for strapdown integration of ACOM. Following this, the estimations of an average VCOM from foot velocities is explained, and then the fusion of the two information sources using a complementary filter.

In this implementation, a changing reference frame was defined to enable expression of kinematics in a body centric frame, irrespective of changes in walking direction. Fig. 10 gives a brief overview of the processes described in this exemplary implementation, where we fuse information from the pelvis IMU and the feet IMUs to obtain a drift free VCOM. The ACOM were first transformed to the pelvis frame by asking the subjects to bend forward, and using this axis of rotation as the Y axis of the pelvis frame. Then a dedicated fusion filter was used to transform to the intermediary

integrating the angular velocities at the pelvis IMU. To reduce drift in this orientation, it is assumed that on average the ACOM measures gravity.

Finally, the ACOM were transformed to

using the pre

determined

Hence, for

, the X axis is positive in the forward direction, and Z is positive upwards along the vertical. These steps have been described in detail in the second exemplary implementation.

As seen in Fig. 10, gravity was removed from ACOM in

(acc^cs(k)), and then strapdown integrated using the direct and reverse integration method (DRI) to obtain velocity (VCOM_sdi). The velocities at the beginning and end of trial were set to 0, as required by the DRI method. VCOM_sdi is a time varying velocity estimate with drift that accumulates over time. A high pass 2^nd order zero lag Butterworth filter was applied to obtain the VCOM_hf., thus providing CoM Velocity using strapdown integration.

On average, the movement of the feet encode information about the movement of CoM. Hence, a low frequency estimate of the VCOM can be approximated from foot velocities. Note that foot velocity estimates are drift free due to zero velocity constraints. As seen in Fig. 10, the velocities of left and right feet (vel_fl and vel_fr) were summed, and averaged. A low pass 2^nd order zero lag Butterworth filter was applied to obtain the VCOM_lf.

The cut off frequencies used for the VCOM_hf and VCOM_lf were the same and after a preliminary analysis, the optimal values were found to be 0.5, 0.2, and 1.4 Hz for X, Y, and Z axes respectively. The VCOM_hf and VCOM_lf were then fused to obtain the instantaneous VCOM_est., thereby obtaining CoM Velocity from Foot Velocities.

The subjects were asked to perform each of the following walking tasks four times in their preferred walking speed:

1. Normal Walk (NW): The subject walked at for 5 m.

2. L Walk (LW): The subject walked 3m and then turned 90° to the right and walked another 2m.

3. Walk and Turn (WT): The subject walked 5m and then turned and walked back to start position.

4. Walk and Turn Twice (WT2): The subject performed WT and then turned again and walked another 5m.

5. Slalom Walk (SW): The subject walked in a slalom. Two pylons were placed to guide them.

6. Asymmetric Walk (AW): The subject walked in an asymmetric manner. The instruction given was to induce a stiff left knee and make shorter steps using the right side.

A three IMU setup was used as seen in Fig. 2: One Xsens™ IMU was mounted on the sacrum using a strap, and one was placed on each foot on the midfoot region. A MT Manager was used to read the data from the IMU wirelessly, which was sampled at 100 Hz. A VICON^© motion capture system was used as the reference. Markers were placed on the following locations on both the left and right hmbs: anterior superior iliac spine, posterior iliac spine, the second and fifth metatarsal, and heel. One marker was also placed on each IMU. The VICON^© was sampled at 100 Hz.

The CoM position obtained from VICON^© (PCOM_ref) was assumed to lie at the centroid of the pelvis markers. The PCOM_ref was differentiated and low pass filtered with a 2^nd order zero phase Butterworth filter of cut off

10 Hz to obtain the 3D VCOM_ref. These were transformed to the

, determined independently using the foot position data of the VICON^®. The data includes walking trials by three healthy males. The mean height, weight, and age was 1.74 ± 3 m, 79.3 ± 9 kg, and 25 ± 3.5 years respectively. Leg length was 94 ± 3 cm. All participants signed an informed consent before the experiment. The study was conducted in accordance with the Declaration of Helsinki, and the protocol was approved by the Ethical Committee of the faculty.

Fig. 11 shows an example of VCOM_hf 21 (blue line) and VCOM_lf 22 (red line) is shown for the WT2 task. The plot includes gait initiation, turning and termination. Following this, Fig. 12 shows an example of the estimated instantaneous VCOM_est 23 (blue line) compared with the reference VCOM_ref 24 (red line). Here, we also depict the strapdown integrated VCOM_sdi (thin black line), which is seen to clearly drift. In both Fig. 11 and 12, the subject makes two 180_l turns which are denoted by shaded light red rectangles. Also, in both figures, each subplot corresponds to an axis of the

Table VII compares the average root mean square of the error between the VCOM_est and VCOM_ref across all subjects for different walking tasks as both absolute, and percentage error normalized to the range of VCOM_ref.

Table VII. Absolute and percentage root mean square error between VCOM_est and VCOM_ref

The current method shows the feasibility of estimating 3D VCOM during variable gait. Although, the complementary method is similar to Sabatini et al. 2016, they used machine learning to estimate the average VCOM. Our approach discards the need for a training step. Fig. 11 shows the complementary information present in the VCOM_hf and VCOM_lf. VCOMi_f derived from the foot velocities encodes the trend and VCOM_hf has information regarding a drift free change in this trend. In Fig. 11 and 12 the kinematics are expressed in

, and hence, the velocities remain positive in the X axis even as the subject makes two 180⁰ turns, during the shaded regions. Note that in Fig. 12, the drift in the vertical VCOM_sdi is quite limited compared to the other axes, but exists, as evident during the last few steps. Further, the VCOM_ref shows more drastic jumps during the turns, as compared to VCOM_est, more clearly seen in the Y axis due to transformation to the . As we account for the

per step k, we can represent the kinematics in a fixed global frame, or the frame of any other required step.

Table VII shows that the errors are on average 13.1 ± 2.2 % of the range of VCOM across all axes and walking tasks. The errors seem to be largest for the SW task. The error margins are quite low overall, about less than 19% of the range in the worst case. The algorithm has lower errors for variable walking when compared to the results of Sabatini et al. 2016. The apphcability of the method however, would dependent on the application, and user’s acceptable error margins. Note that the cut offs used in the complementary filter was optimized across all subjects. The errors found could be further lowered if this was optimized per subject. A drawback of this method is that it employs a DRI method for integration, which requires knowledge of the final state of the velocities.

Note that the use of foot IMUs is two-fold: for defining the

as well as obtaining a low frequency VCOM information. The fused VCOM_est can further improve estimates of the position of COM, and XCOM. A three IMU setup could therefore be sufficient to obtain foot and CoM kinematics.

Fourth exemplary implementation In this exemplary implementation we describe how the data estimated in the previous implementations can be used to track the feet and CoM positions. A brief overview of the steps are shown in Fig. 13.

An extended Kalman filter EKF is built to track the velocity and position of the three segments in the current global or alternatively current step frame

. For the pelvis IMU, we assume that the pelvis moves with the CoM, and that all mass is concentrated at this point. The EKF can therefore track the position and velocity of the CoM. Therefore, the pelvis segment will be hereto referred as the CoM. The state vector of the EKF is denoted as x and its covariance matrix as P. The states required for each segment are its 3D position p, and 3D velocity v. The state vector can be thus written as:

The following text expands on process overview shown in Fig. 13. It describes the a-priori estimate determined using strap down integration of the accelerations measured at each segment. Then, the special

biomechanical constraints applied to each segment are shown. The implementation of the CMP assumptions to reduce the drift between the three segments is also described.

1) Initialization: Before applying the EKF, the states for each segment and their covariance noises have to be initialized. The right foot is assumed to be the origin. The location of the CoM, and the left foot is input using the values from the reference system. They could also alternatively be measured using tape measure. All initial velocities were set to zero.

The initial noise values are set to be an arbitrary value. 2) Prediction: We start by making models for each of the elements in the state vector. The acceleration for each segment is given as

The position and velocity can be estimated from acceleration using integration as follows:

The Kalman filter prediction equation is given as

and the covariance matrix is predicted using

where, Q is the process noise covariance matrix.

3) Measurement Update: The measurement updates used to reduced drift in the estimation of position and velocity of the feet and CoM are as follows. A summary of when these updates are applied is shown in Fig. 13. The measurement used at any instance k to update the state vector is given as

The measurement updates used are as follows.

• Zero Velocity Update : During foot contact instances, as estimated using Skog et al. 2010, the velocity of the feet are assumed to be zero. This can be implemented as a measurement update.

In the above equations, H transforms the state vector to a measurement prediction

and Z denotes the actual measurement

(Welch et al. 1995).

denotes the error associated with this measurement. The same notations are used in the following equations.

• Zero Height Update : During the foot contact instances, we also have information regarding the height of the foot from the floor. Here, we assume that gait occurs over a flat ground.

The measurement matrix H has only one row as it transforms only the Z axis of each foot. Here,

is the initial height and the superscript

(f) denotes either the left (fl) or right foot (fr). This update was only applied to the feet.

· Minimum Height Update : A minimum height check is also performed during the iteration of the EKF. The height of the feet are constantly assessed to see if it is above ground level. In case it drops below the ground, the height is reset to the initial height.

• CoM Velocity Update : The ZV update ensures that the velocity of the feet do not drift due to integration errors. However, the CoM is constantly moving. Therefore, an estimate of the CoM velocity was derived by fusing two complementary sources of information (third exemplary implementation). A high frequency information was derived from an optimally filtered direct and reverse strapdown integration (Zok et al. 2004) of the CoM accelerations using a cut off of 0.6 Hz. Then, low frequency information of the CoM velocity was derived from low pass filtering the average of the foot velocities using the same cut off used for the earlier high pass filter. The two sources were fused to get estimates of the instantaneous 3D CoM velocity:

• CoM Height Update : The CoM position may also drift as there are not many constraints available. However, we can assume that the height of the pelvis will not go beyond certain ranges given an upright gait. An‘average’ height of the pelvis is estimated from the first two steps made. Then, a boundary of ±0.20% is assumed. If the estimate of the pelvis exceeds this boundary, the pelvis height is reset to be the average height.

Alternatively, the height of the CoM can be estimated using a complementary filter method. An optimally filtered direct and reverse strapdown integration (Zok et al. 2004) of vertical CoM velocity was used to obtain the changes in CoM height during gait using a cut off of 0.3 Hz to obtain the instantaneous change in height. Then, as the subject does not crouch or jump while walking, the height of the CoM should oscillate around an offset. Assuming an average walking CoM height as 98% of the height during quiet standing showed least errors when validating this method. The average height and the instantaneous change in height were fused to get an estimate of CoM height during walking.

Zero Moment Positions Update or Centroidal Moment

Pivot Update : The CMP update was used to restrict the feet and CoM from drifting apart or towards each other. Fig. 13 summarizes the steps involved; the CMP update was used to first estimate the horizontal position of stance foot with respect to the movement of the CoM. After correcting discrete changes at the start of each swing phase, the CoM trajectory was updated for the corrected foot positions.

The above equation is a mathematical expression of the CMP theory (Popovic et al. 2005). When the 3D components of GRF are known, an estimate of the distance between the CoM and a virtual CMP point under the stance foot can be estimated (first exemplary implementation, Mohamed Refai et al. 2020). Here, F is the 3D GRF in a certain axis, and ax denotes either X or Y axes. The 3D components of GRF estimated from the CoM accelerations (Mohamed Refai et al. 2020) were used to estimate the ratio F_(ax)/F_Z. Further, the height of the CoM was already estimated earlier, and

is the estimate of CoM in the X or Y axis. Note that there are a few assumptions regarding the above equation. For instance, we assume that the virtual CMP position coincides with the stance foot positions tracked by the IMU. Secondly, we assume that the moment of inertia around the trunk is negligible while walking. Thus, during single stance phase, the equation provides the relation between CoM and cmp, the stance foot. This can be used as measurement updates for either foot as follows:

During left swing:

with,

During right swing:

^' with.

The measurement matrices transform only the X and Y axes of the foot positions in the state vector. This update corrects the drift in relative positions between the CoM and stance foot during swing phase. However, this may cause a discrete jump in relative foot distances at the start of the swing phase. To have a smooth change in relative foot distances between subsequent steps, knowledge of the relative foot distances at the end of the preceding step was used to update the relative foot distances at the beginning of the subsequent swing phase.

Start of left swing:

with,

Start of right swing:

with,

In above equa

t

ions,

are the respective foot positions at the end of the preceding step in the axis ax. As in earlier equations, the measurement matrices in the above equation transform only the X and Y axes of the foot positions. This correction of relative foot distances requires a final update of the CoM position following the CMP theory. The CMP theory equation was adapted to obtain the CoM from foot estimates as follows:

During left swing,

represents the right foot, and vice- versa for the right swing. During these instances, the CoM position was improved using the following measurement update:

4) Reinitialising for step k+1: After the state vector has been updated using the prediction and measurement updates, the trajectory of the segments are known for the current step k. The

was adjusted using the improved estimates of the foot positions. For the next step k+1, accelerations are transformed to the current step frame

, and the sensor fusion filter is re-iterated.

Numerical results are shown in the following Figures 14-20. Figure 14 shows a trajectory of CoM Height for a subject performing a WT task estimated from the pelvis IMU 25 seen as solid blue line. The dotted red line is the VICON reference 26. Figure 15 shows a step-wise comparison of foot and CoM positions by the current system and VICON reference in the frame. The subject performs a WT task, with the turning step highhghted with a light red shaded background 27. The darker red 28, blue 29, and green 30 circles points denote the PGL values, and the lighter red 28’, blue 29’, and green 30’ circles points denote the reference values. In either case, positions of the right foot, left foot, and CoM are shown using red 28, 28’, blue 29, 29’, and green 30, 30’ circles. The step numbers are mentioned within each plot. Figure 16 shows comparison of step lengths measured using the current system with that of the optical reference using correlation and Bland- Altman plots. Red circles 31 are turning steps. *p < 0:05. Figure 17 shows comparison of step widths measured using the current system with that of the optical reference using correlation and Bland- Altman plots. Red circles 32 are turning steps. *p < 0:05. Figure 18 shows comparison of CoM widths measured using the current system with that of the optical reference using correlation and Bland-Altman plots. Red circles 33 are turning steps. *p < 0:05. Figure 19 shows comparing distributions of right and left step lengths for the NW task, and Figure 20 shows comparing distributions of right and left step lengths for the AW task. *p < 0:05.

First, an example of the CoM height estimated using the complementary filter approach in is shown in Fig. 14 as the solid blue line 25. The dotted red line 26 shows the measurement by the VICON. Here, the subject is performing a WT task and makes the 180 degree turn around 25 seconds. Table VIII shows the errors (rRatX , and rRatY ) in estimating the ratio of forces and that of CoM height (RCoMZ ) for the different walking tasks. The CoM height estimations were significantly (p < 0.05) correlated with the reference values for all tasks.

Fig. 15 shows a graphical step-wise comparison of the feet and CoM positions estimated by the PGL compared with the reference VICON for the same trial as shown in Fig. 14. In Fig. 15, the darker circles 28, 29, 30 denote estimations by the PGL, and the lighter circles 28’, 29’, 30’ denote the measurements by VICON . For both systems, the trajectory of the right foot is shown in red 28, 28’, the CoM in green 30, 30’, and the left foot in blue 29, 29’. Each subplot is a top-down view of a step expressed in the current step frame _¾rcs. The X axis of each plot corresponds to the X axis of the _¾rcs, and similarly for the Y axis. As each step is represented in its own _¾rcs frame, they all progress to the right, even during turns. In the first subplot, the left foot moves first. Then, the subject can be seen to make consecutive steps, until step 6, where they prepare for the 180 degree turn. The turning step after step 6 has a shaded light red background 27. Although the PGL shows deviations during the turn when compared to the reference, it converges to the reference values two steps after.

Table VIII displays the errors (RightX , RightY , LeftX , LeftY , CoMX , and CoMY ) in estimating the horizontal feet and CoM positions for each step. These are an average across all steps in a walking task excluding the turning steps. Turning steps were those that made a 60 degree or larger change in direction when compared to the preceding step. Table IX also summarizes the difference in relative distance between the feet (ED) at the end of each step. Across all walking tasks, this was found to be 8.6 ± 1.5 cm on average.

The estimates of SL, SW, and CW for all tasks except AW are compared against the reference using Fig. 16, 17, and 18 respectively. Some steps had particularly large SL or CW values measured by the VICON than the average, thereby skewing the distribution as outliers. They were removed based on the interquartile range of the distribution of the differences between PGL and VICON estimates for each parameter (SL, SW, and CW). Fig. 16, 17, and 18 distinguish between a step with relatively straight heading, denoted as a blue filled square 31’, 32’, 33’ and turning steps denoted by a red filled circle 31, 32, 33. In each of the three Figures 16, 17, and 18, the left subplot shows the correlation between the PGL estimates and reference. The dotted gray hne 34 is the identity line, the sohd blue hne 35 is the linear fit for the straight steps, and sohd red line 36 is the hnear fit for the turning steps. The legends in each figure indicate the correlation between the PGL and reference and its significance denoted by a star. The average root mean square of the errors is also shown. Further, in Fig. 16, 17, and 18, the right subplot shows the Bland-Altman plot. The sohd hnes 37, 37’ denote the median difference between the two systems, and the dotted lines 38, 38’, 39” denote the 95% limits of agreement (LoA). The LoA for SL were found to be [-20.1 14.9] cm and [-19.9 21.9] cm for the straight and turning steps respectively. Similarly, for SW it was [-15.3 12.3] cm and [-36.6 17.6], and finally for CW it was [-9.4 2.7] and [-17.7 18.1] cm.

Finally, the feasibility of PGL in differentiating symmetric from asymmetric gait is shown in Fig. 19 and 20. These figures compare the boxplot distribution of SL on the left and right side. Fig. 19 plots this comparison for the NW task, where the subjects walked symmetrically, and Fig. 20 plots this for the asymmetric AW task. In both figures, the left subplot shows the distributions measured by the reference system, and the right subplot shows that of the PGL. Fig. 20 shows that the both, the reference and PGL, find significant differences in SL between the right and left side for the AW task.

Table VIII: Root Mean Square of the differences in estimating ratio of forces rRatx, and rRaty, and CoM height RCoMz

Table IX: Average Root Mean Square of the errors in horizontal positions of the feet and CoM, and the differences (ED) in relative foot distances at the end of each step.

The EKF mentioned above is to be applied in an iterative manner per step. Updates have been formulated on biomechanical constraints for each segment. This will allow the estimation of drift free positions of both feet and CoM over time using only three IMUs.

The invention is not restricted to the embodiments described above. It will be understood that many variants are possible.

These and other embodiments will be apparent for the person skilled in the art and are considered to fall within the scope of the invention as defined in the following claims. For the purpose of clarity and a concise description features are described herein as part of the same or separate embodiments. However, it will be appreciated that the scope of the invention may include embodiments having combinations of all or some of the features described.

Claims

1. A method for estimating positions of a subject’s feet and centre of mass relative to each other during gait, comprising the steps of:

- collecting (110) first measurement data from a first inertial measurement unit (2a) located at a first foot or first shank of the subject (P);

- collecting (120) second measurement data from a second inertial

measurement unit (2b) located at a second foot or second shank of the subject;

- collecting (130) third measurement data from a third inertial measurement unit (2c) located at a pelvis of the subject, and

- evaluating (140) relative positions of the first and second foot and the center of mass of the subject over time, using the first, second and third measurement data characterized in that the method further comprising the steps of:

- assessing a first phase during gait wherein a first foot is resting on a support structure supporting the first foot, and a second foot is moving;

- determining a first estimation of the second foot position relative to the first foot, after a step of the second foot, applying the assumption that a moment around a center of mass of the subject vanishes.

2. A method according to claim 1, further comprising a step of evaluating the velocity of the center of mass of the subject over time.

3. A method according to claim 1 or 2, further comprising a step of evaluating force vector components exerted at the center of mass of the subject, using the first, second and third measurement data.

4. A method according to claim 1 or 2, further comprising a step of evaluating force vector components exerted at the center of mass of the subject, using only the third measurement data.

5. A method according to claim 1, further comprising a step of determining a second estimation of the second foot position relative to the first foot, after the step of the second foot, by applying strap -down

navigation estimation of displacement of the second foot during the step of the second foot, combined with initial and final conditions that follow from zero velocity update, and an estimated position of the second foot relative to the first foot, before the step of the second foot,.

6. A method according to claim 5, further comprising a step of improving an estimation of the second foot position relative to the first foot, after the step of the second foot, by combining the first and the second estimation of the second foot relative to the first foot, preferably taking into account uncertainties of the first and second estimation, preferably using a specifically designed sensor fusion filter.

7. A method according to any of the preceding claims, further comprising the steps of:

- assessing a second phase during gait wherein the first foot is moving and the second foot is resting on the support structure supporting the second foot;

- performing a first estimation, a second estimation and an improved estimation of the first foot position relative to the second foot, after the step of the first foot, by performing the steps defined in preceding claims 6-8 wherein the first foot and the second foot have interchanged.

8. A method according to any of the preceding claims, further comprising a step of estimating the position of the center of mass of the subject as a function of time relative to the first foot when the second foot is moving and/or relative to the second foot when the first foot is moving applying the assumption that a moment around a center of mass of the subject vanishes.

9. A method according to claim 8, further comprising a step of improving the estimation of the position of the center of mass by including the third measurement data from the third inertial measurement unit located at a pelvis of the subject, preferably including a strap-down navigation estimation.

10. A method according to any of the preceding claims, further comprising a step of estimating the position of the center of mass of the subject during double stance based on the third measurement data from the third inertial measurement unit located at a pelvis of the subject, preferably including a strap -down navigation estimation.

11. A method according to any of the preceding claims 8-10, wherein the position estimation of the center of mass is improved by fusing all information about center of mass movements available over time, preferably using a specifically designed sensor fusion filter.

12. A method according to any of the preceding claims 8-11, further comprising a step of estimating a velocity of the center of mass at any moment time based on a position estimation of the center of mass as a function of time and optionally the third measurement data from the third inertial measurement unit.

13. A method according to any of the preceding claims, further comprising a step of estimating a velocity of the center of mass from the velocity estimates of the feet.

14. A method according to the preceding claims 12 or 13, wherein the step of estimating a velocity of the center of mass is improved by fusing information from the third measurement data and the velocity estimates of the feet.

15. A method according to any of the preceding claims, wherein a location, in particular the height of the center of mass of the subject is estimated based on measurement data from an inertial measurement unit located at the pelvis of the subject, preferably by applying strap-down inertial estimation and using biomechanical constraints regarding the height of the center of mass.

16. A method according to claim 15, wherein a ground reaction force exerted on the center of mass of the subject is estimated from acceleration data measured by the inertial measurement unit located at the pelvis of the subject.

17. A method according to claim 16, further comprising a step of estimating an orientation of the inertial measurement unit, preferably using an error extended Kalman filter.

18. A method according to any of the preceding claims, wherein foot kinematics and kinetics are determined relative to a frame of reference that is stationary relative to the pelvis or a direction of gait of the subject.

19. A method according to any of the preceding claims, further comprising a step of estimating foot kinematics and kinetics over time, based on measurement data from a first inertial measurement unit located at a first foot, measurement data from a second inertial measurement unit located at a second foot and measurement data from a third inertial measurement unit located at the pelvis of the subject.

20. A method according to claim 19, further including a step of estimating balance metrics, based on the kinematics of the foot and the center of mass of the subject, by

- estimating an extrapolated center of mass, which includes the information of speed and direction of walking, and

21. A method according to claim 19, wherein the step of estimating foot kinematics and kinetics over time includes applying sensor fusion filters, such as an extended Kalman filter to measurement data.

22. A method according to claim 19 or 20, wherein the foot kinematics over time are updated using the first distance between the position of the first foot and the center of mass of the subject, and/or the second distance between the position of the second foot and the center of mass of the subject.

23. A method according to any of the preceding claims, wherein the foot kinematics over time are updated using biomechanical constraints.

24. A system (1) for estimating a subject’s feet relative to each other during gait, comprising a first inertial measurement unit (2a) for location at a first foot or first shank of a subject (P), a second inertial measurement unit (2b) for location at a second foot or second shank of a subject (P), and a third inertial measurement unit (2c) for location at a pelvis of the subject, the system further comprising a processing unit (3) arranged for:

- collecting (110) first measurement data from the first inertial

measurement unit (2a);

- collecting (120) second measurement data from the second inertial measurement unit (2b);

- collecting (130) third measurement data form the third inertial

measurement unit (2c), and

- evaluating (140) relative positions of the first and second foot and the center of mass of the subject over time, using the first, second and third measurement data, characterized in that the processing unit is further arranged for:

25. A computer program product for estimating a subject’s feet relative to each other during gait, the computer program product comprising computer readable code for causing a processor (3) to perform the steps of: - collecting (110) first measurement data from a first inertial measurement unit (2a) located at a first foot or first shank of the subject (P);

- collecting (120) second measurement data from a second inertial

measurement unit (2b) located at a second foot or second shank of the subject (P);

- collecting (130) third measurement data from a third inertial measurement unit (2c) located at a pelvis of the subject (P), and

- evaluating (140) relative positions of the first and second foot and the center of mass of the subject over time, using the first, second and third measurement data, characterized in that the computer program product also comprises computer readable code for causing the processor to perform the steps of:

26. A computer program product according to claim 25, comprising computer readable code for causing a processor to perform the further step of:

- building a sensor fusion filter for fusing information from strapdown inertial navigation and relative foot positions to improve foot tracking as a function of time.

27. A computer program product according to claim 26, comprising computer readable code for causing a processor to perform the further step of:

- displaying the information in a frame of reference defined with respect to the pelvis or the direction of gait.