Orientation Estimation Through Magneto-Inertial Sensor Fusion: A Heuristic Approach for Suboptimal Parameters Tuning

Magneto-Inertial Measurement Units (MIMUs) are a valid alternative tool to optical stereophotogrammetry in human motion analysis. The orientation of a MIMU may be estimated by using sensor fusion algorithms. Such algorithms require input parameters that are usually set using a trial-and-error (or grid-search) approach to find the optimal values. However, using trial-and-error requires a known reference orientation, a circumstance rarely occurring in real-life applications. In this article, we present a way to suboptimally set input parameters, by exploiting the assumption that two MIMUs rigidly connected are expected to show no orientation difference during motion. This approach was validated by applying it to the popular complementary filter by Madgwick et al. and tested on 18 experimental conditions including three commercial products, three angular rates, and two dimensionality motion conditions. Two main findings were observed: i) the selection of the optimal parameter value strongly depends on the specific experimental conditions considered, ii) in 15 out of 18 conditions the errors obtained using the proposed approach and the trial-and-error were coincident, while in the other cases the maximum discrepancy amounted to 2.5 deg and less than 1.5 deg on average.


I. INTRODUCTION
Instrumented movement analysis is used to provide a quantitative kinematic description of the locomotor apparatus in subjects with and without motor impairments [1]. Optical stereo-photogrammetry (SP) is widely accepted as the gold standard for human movement analysis applications; however, SP systems are expensive, and their use generally limited to paradigmatic motor tasks performed in motion analysis laboratories.
In recent years, the use of low-cost Magneto and Inertial Measurement Units (MIMUs) in movement analysis has gained traction thanks to some of their remarkable features. Miniaturized MIMUs integrate a triaxial accelerometer, a triaxial gyroscope and a triaxial magnetometer in a small device, they can be used to record movement data both indoors and outdoors and are inexpensive. Nonetheless, some limitations must be considered when MIMUs are used for human movement applications [2]. Whereas, specific acceleration and angular velocities can be directly measured, the absolute orientation of a MIMU can be only estimated by combining the accelerometer, gyroscope and magnetometer readings [3]. This sensor fusion process is typically performed in two steps. The first step involves the time integration of the differential kinematic equation, which links the orientation update with the angular velocity, to obtain a first approximation estimate of the orientation change. The resulting drift, caused by the integration of the slow-varying bias affecting the gyroscope readings, is then corrected by using the accelerometer and magnetometer data which provide a gyro-free estimate of the absolute orientation. The accelerometer is used for estimating the inclination by sensing the Earth's gravity in its local coordinate frame, whereas the magnetometer can be employed to estimate the orientation of the MIMU with respect the Earth's magnetic North, acting as a compass. However, the inclination estimate is accurate only when the MIMU is stationary, since when the MIMU is in motion, the accelerometer signals are the result of the combination of gravity and MIMU's acceleration and the effect cannot be separated unless additional information are used. Moreover, the heading information estimated from the magnetometer output may be corrupted by the magnetic field distortions, which may arise from the surrounding metallic objects and electric appliances, therefore limiting its indoor use. Finally, electronics noise of the MIMU sensors, misalignment, nonorthogonality of the sensor axes, sensitivity to changes in temperature, further affect the accuracy of the estimates [4].
To compensate and minimize the effects of the various sources of errors, several sensor fusion algorithms for orientation estimate, from both Kalman and complementary filter families, have been designed in the last two decades [5]- [19]. Recently, also machine and deep learning approaches have been proposed to estimate the orientation [20], [21]. In order to work properly, each sensor fusion algorithm requires the tuning of the values of a given number of parameters. It has been demonstrated that the specific choice of the parameters values could greatly affect the orientation estimate accuracy [21], [22]. Ludwig and Burnham in [23] stated that the algorithm from Mahony et al. [13] exhibited good results (i.e. 1 deg for inclination and 4.2 deg for heading) as long as the optimum parameter set (i.e. the values which provide the lowest absolute errors) was used, and this set appeared to be unique to each data set or type of motion. It is widely recognized that a major problem when implementing any sensor fusion algorithms is to define the most suitable values for the parameters required [2], [21]- [26]. In fact, effective parameter values selection depends on different intrinsic and extrinsic factors; among which the most important are the type and amplitude of motion (which are reflected in the magnitude of the body accelerations) [12], [22], [27], [28], sensors noise character-istics [15], sensor fusion algorithm convergence speed [29], and ferromagnetic disturbances [12]. Recent literature suggests that no well-established solutions for the definition of the parameter values have been found yet [2], [23], [26], [30]. From a practical point of view, "appropriate" set of values are commonly chosen either following the suggestions provided by the proposers of a sensor fusion algorithm in their original implementations based on their dataset or by minimizing the overall errors between the estimated and errorfree orientation provided by a gold standard technique [18], [30]. The latter approach, also known as trial-and-error (or grid search), was also adopted by Madgwick and colleagues in [15]. Unfortunately, the trial-and-error approaches are timedemanding, require a good level of expertise, and do not guarantee for generalization. In a recent paper, Ludwig and Jiménez proposed a genetic algorithm to find the parameter value to optimally weight the orientation computed from the gyroscope and from accelerometer-magnetometer [25]. Even in that case, the optimization was performed by exploiting the orientation reference acquired from the proprietary Xsens filter which should rather be considered a silver standard. In addition, also the neural network proposed by Seel et al. was trained using a marker-based optical ground-truth [21].
In this study, we proposed an original approach for suboptimal parameter values identification designed to work without relying on any reference data. The approach assumes that two MIMUs rigidly connected are expected to show no orientation difference during motion. This feature may be convenient when developing MIMU-based methods for monitoring human movement outside the laboratory. The validity of the proposed "rigid-constraint" approach was tested and assessed on the popular open source complementary filter (MAD) by Madgwick and colleagues [15], which requires to set only one parameter and using the orientation provided by a SP system as gold standard. The proposed approach was tested under different experimental scenarios by varying the commercial product employed (Xsens, APDM, Shimmer), the angular rate (slow, medium, fast), and the motion dimensionality (2D vs 3D). We hypothesized that different values of the filter parameter values would be most appropriate for different experimental scenarios, as already anticipated by Mahony and Ricci in [13] and [2], respectively.

II. MATERIALS AND METHODS A. Experimental Set-Up
To test the robustness of the proposed approach to different sensor noise characteristics, we selected three different commercial MIMUs: Xsens -MTx, APDM -OPAL, and Shimmer Sensing -Shimmer3. Two MIMUs for each product were aligned on a wooden board (Fig.1). All lines were drawn using a T-square to ensure a high accuracy in the MIMU and marker positioning. The alignment error due to orthogonal tolerance of the instrument was estimated to be lower than 0.2 deg. Relative distance between MIMUs was 50 mm. The board was also equipped with eight reflective spherical markers (dia = 14 mm, minimum inter-distance of 85 mm). Three markers positioned centrally on the board were used to define the SP Local Coordinate System (LCS) aligned with the MIMUs LCS. Marker redundancy was exploited to minimize the effects of the stereo-photogrammetric errors on the orientation estimation by means of the Singular Value Decomposition technique (SVD) [31]. A Vicon T20 (software Nexus 2.7) SP with 12 cameras was used to obtain the gold standard orientation. One force platform (AMTI, sampling frequency = 1000 Hz) integrated with the Vicon system was used to synchronize MIMUs and SP system by generating a series of mechanical shocks.

B. Experimental Protocol
To limit the temperature effects on the sensor readings, a ten-minute instrument warm-up was carried out before measurement. Measures were performed in static and dynamic conditions. The static measure consisted in a 60-second acquisition with the board fixed horizontally on an aluminum tripod positioned over the force platform. The dynamic recordings consisted in a sequence of rotations of the board. An operator moved the board to manually change its orientation, first by rotating it about the three axes (x, y, x) one at a time from 0 to 180 degrees and back to 0 for five times (2D motion) and then by performing a free rotation to span the three rotational degrees of freedom simultaneously. The duration of the 3D movements was approximatively the same as that of 2D. The board was then repositioned on the tripod for 60 s. The board was hit at the beginning and at the end of the recording for synchronization purposes. The protocol was executed at three different angular rate conditions (whose rms values were assessed in post processing): slow (rms(ω) = 120 deg/s for a total of 70 s), medium (rms(ω) = 260 deg/s for a total of 45 s), and fast (rms (ω) = 380 deg/s for a total of 30 s). The experiments were conducted at about 20 • C in an approximately 1 m 3 volume. Data from Xsens (MT Manager Version 1.7) and Shimmer (Consensys v.1.5.0) systems were recorded at 100 Hz, whereas data from OPAL (Motion Studio Version 1.0.0.201712300) were recorded at 128 Hz. Marker trajectories were recorded at 100 Hz. The full processed MIMU dataset together with the SP orientation (computed as described in section II.E) and the videos have been uploaded on IEEE DataPort at http://dx.doi.org/10.21227/b23b-rx94. This protocol was designed to explore both 2D and 3D motions. The assessment of planar movements is of particular interest in gait analysis as leg motion mainly occurs in the sagittal plane. Multi-axial motions are usually encountered when assessing upper limb movements during functional activities [32].

C. MIMU-Based Orientation Algorithm
The "rigid-constraint" approach for determining the parameter values was applied and tested on the MAD sensor-fusion filter proposed by Madgwick and colleagues [15]. The MAD algorithm belongs to the complementary filter (CF) family: the orientation obtained from the gyroscope readings is fused in the frequency domain with that computed from accelerometer/magnetometer readings according to their complementary spectral characteristics. A high-pass filter is adopted to reduce the effect of the slow-varying bias affecting the gyroscope, whereas the high-frequency noise which corrupts the accelerometer and magnetometer is low-pass filtered. The cut-off frequency value must be the same for both filters and represents a trade-off between the two preserved bandwidths [18]. The final orientation is the weighted sum of the two filtered signals. The CF, in general, does not include any statistical description of the noise to be considered into the algorithm model.
MAD addresses the problem of the orientation estimate in the quaternion form instead of orientation matrix or Euler angles, because it has a lower computational demand and it allows to avoid singularities such as the Gimbal Lock [30]. The quaternion is a four-terms parametrization, given an angle of rotation θ and the rotation axis n = [n x , n y , n z ] T around which the rotation occurs; accordingly, MAD expressed the orientation of the Global Coordinate System (GCS), defined to have one axis aligned with the gravity and one to the magnetic North LCS with respect to the LCS, as in (1): The first step of any sensor fusion algorithm involves the numerical integration of the angular velocity over time to obtain the change of orientation. The kinematic equation relating the angular velocity expressed in the LCS of a moving body, and the temporal derivative of the orientation GC S LC Sqω,t is presented in (2): where GC S LC S q t −1 is the orientation estimated at the previous time-step, and ⊗ represents the mathematical operator for quaternion multiplication. The quaternion describing the change of orientation GC S LC S q ω,t can be computed according to II-D under the two hypotheses of a constant angular velocity within sampling period t, and of a sufficiently small t (first order approximation): Time-integration errors affecting GC S LC S q ω,t due to the gyroscope bias can be corrected using a quaternion GC S LC S q ∇,t computed with gradient descent algorithm based on the accelerometer and magnetometer readings. In particular, the magnetometer readings are projected on the horizontal plane and then used to estimate the heading component of the orientation (the rotation around the MIMU's LCS z-axis, also known as yaw angle). In this way, the ferromagnetic disturbances could affect only the heading without influencing the estimate of the attitude (the inclination of the MIMU, also jointly known as the combination of roll and pitch angles).
The fusion process between GC S LC S q ω,t and its correction quaternion GC S LC S q ∇,t gives GC S LC S q t (the orientation at the current time-step) and it is governed by the weighting factor β as follows: A larger value of β gives more weight to the orientation computed from the accelerometer and magnetometer, so as to limit the orientation drift; however, this makes the resulting orientation more sensitive to body acceleration and ferromagnetic disturbances. It has to be said that in MAD it is not possible to decouple the weight given to the contribution from the accelerometer and magnetometer, as opposed to other works (e.g. [14] and [18]

D. Rigid-Constraint Approach for the Identification of the Suboptimal Parameter Value
The purpose of this section is to describe a method, and the underlying assumptions, for finding a suitable value of β starting from the orientation of two aligned MIMUs without using the ground-truth data.
In absence of errors, when motion data is recorded from two aligned MIMUs attached to the same rigid body at a relative distance r (the distance between the relevant LCS origins), no orientation difference would be observed between the two MIMUs throughout the motion. If the value of β that minimizes the relative orientation difference between two MIMUs during the recorded motion, also guarantees small absolute orientation errors then the same value for β could be used without requiring any orientation reference to be available.
The hypothesis relies on the assumption that the sources of orientation errors which affect the sensors embedded in two different MIMUs are different. In the following paragraphs this hypothesis will be discussed separately for two gyroscopes, two accelerometers, and two magnetometers.
In this regard, it is well known that low-cost MEMS gyroscopes are affected by non-stationary bias and can exhibit non-negligible run-to-run changes [33], [34] due to the fact that variations due to temperature changes are not completely compensated for. The bias corrupting the gyroscope readings may be defined as the gyroscope output in the absence of motion [35] and its integration leads to an orientation drift which grows unbounded. To limit the variation of the gyroscope bias due to non-linear temperature influence, it is recommended to perform a sensor warm-up prior to experiments [22], [33], [36]- [38]. Another effective solution is to estimate the bias during a static trial just before the experiments (especially if these experiments need to be repeated at intervals of few hours), and subtract it from the gyroscope readings during motion [18], [30]. Unfortunately, gyroscope bias is not stationary, and some residual bias cannot be eliminated this way. It is however reasonable to assume the residuals of two similar gyroscopes to be independent and uncorrelated sources of error, being mounted on different chips, with different sensing elements, read-out circuits and signal conditioning circuits [39].
As previously mentioned, the complementary information provided by the accelerometer and the magnetometer are used to compensate for the drift generated during the angular velocity integration. However, within the CF framework, the accelerometer provides a reliable estimate of the MIMU inclination (pitch and roll) only in static conditions. Specifically, the accelerometer measures the specific force (a s ) which is the difference between the coordinate acceleration (a) and gravity vectors (g), as reported in (5).
Hereinafter, all quantities are expressed in the sensor LCS (superscripts are omitted for the sake of brevity), except for g which is expressed in the GCS, as stated by the superscript in (5). It is worth noting that during movement a and the projection of g cannot be separated using accelerometer measurements only.
Considering two accelerometers attached on the same rigid body at a distance r (Fig.2) the linear acceleration of the LCS origins of the two sensors are related by: where ω andω are the body angular velocity and acceleration. It is worth underlining that (6) assumes that the LCSs of the two accelerometers are perfectly aligned. This may not be entirely true in reality, due to a residual misalignment between the sensor axes and the case axes and the uncertainty of the manual alignment of the two cases. However, these two contributions may be neglected due to their limited magnitude, as reported in section II.A and in Table. IV. By combining (5) and (6), in case of null relative orientation difference, it yields: The difference between the two accelerometer readings increases with the distance between the sensor's origins and with the angular velocity and acceleration of the rigid body. According to (7), the two accelerometers would sense the same specific force only if the LCSs share the same origin.
The magnetometer is used to measure the direction of the local magnetic field to correct the MIMU heading. Without ferromagnetic disturbances, two magnetometers aligned on the same rigid body would sense the same Earth magnetic field (Fig. 3a). When magnetic disturbance (electrical motors and metal objects) is present, the magnetic field sensed by the two magnetometers may be different depending on their relative position with respect to the source of magnetic disturbance (Fig. 3b). In addition, during the rigid body movement, a stationary source of magnetic disturbance is sensed by the two magnetometers as a time-variant and unknown source of uncertainty [14]. It is reasonable to expect that, in the presence of a distorted magnetic field, the readings from the two magnetometers will differ by a quantity which also depends on their relative distance r [40].
The magnetic field measured by each magnetometer could be expressed as follows: Given h 1 and h 2 , the magnetic fields measured by the two MIMUs, and provided that h Ear t h has the same components (in the same geographic area) when the MIMUs have the same orientation, their difference is expressed as: It is possible to assess that h 1 and h 2 are coincident in two cases: 1) when the terms h 1ex t and h 2ex t are negligible due to the absence of external source of magnetic disturbances and 2) when r approaches to zero so that h 1ex t is close to h 2ex t even when magnetic disturbances are present.
Finally, the electronic noise of the accelerometer, gyroscope, and magnetometer can be also considered an uncorrelated source of uncertainty for the two MIMUs [39].
In summary, the proposed approach for the parameter tuning relies on the assumption that the errors and disturbances affecting the two MIMUs attached on the rigid body can be different during the recorded motion. Whereas for two different gyroscopes the errors can be considered uncorrelated, when the rotation component could not be neglected and ferromagnetic disturbances are present, for accelerometers and magnetometers differences between relevant signals are expected to grow as the relative distance vector r increases.

E. Data Processing and Error Computation
Data processing was performed in MATLAB, software release 2019a (The MathWorks Inc., Natick, MA, USA). Marker trajectories, after being labelled and gap-filled in Nexus, were low-pass filtered using a zero-phase Butterworth filter of the 6 th order (cut-off frequency set to 6 Hz as suggested in [30]) to remove high frequency noise. The MIMU signals and the marker trajectories were first delimited by finding the two acceleration and force peaks recorded by the vertical axes of the accelerometer and force plate, respectively. Then, the signals were resampled at 100 Hz using a linear interpolation technique. The synchronization was refined by aligning the signals of each MIMU with the SP data after determining the delay (by means of the cross-correlation technique) between each MIMU angular velocity and that obtained from the marker trajectories following the approach presented in [41].
The gold standard orientation (q S P_G ) was obtained by computing the LCS orientation with respect to the SP GCS with the SVD technique [31]. The errors in SP orientation can be assumed to be equal to 0.5 deg and this information can be derived from the trigonometry by considering the actual size of the central marker cluster used to define the SP LCS, and that the expected errors on marker position are in the order of 0.1 mm [42]. Data were analyzed only between the synchronization events.
The accuracy of the accelerometer calibration was assessed by aligning each axis of the case along the vertical direction (the alignment was verified with the spirit level embedded in the tripod). Accelerometer measurements were averaged and compared to g = 9.81 m/s 2 . Since the maximum difference never exceeded 0.02 m/s 2 , the accelerometers were considered as properly calibrated [19].
The 60 seconds static acquisition performed at the beginning of each trial was used to characterize the noise of all sensors and the bias of all gyroscopes. The latter was computed as the mean value of the gyroscope readings of each axis [30].
The accuracy of the magnetometer calibration was assessed by performing a movement to cover as much as possible a sphere [43]. Since the maximum difference never exceeded Fig. 3. Two different magnetometers sense the same magnetic field in absence of ferromagnetic disturbances (a). When an external magnetic field is superimposed to the Earth's magnetic field (b) the two magnetometers sense two different magnetic fields due to their relative distance r. 0.1 μT, the magnetometers were considered as properly calibrated.
The procedure used for the data processing is detailed below and reported in the pseudocode (Fig. 4).
Initially, the bias of each gyroscope was removed from the dynamic readings. To minimize the convergence time, the orientation of each MIMU was initialized by means of an algebraic quaternion obtained with the algorithm proposed in [18]. The absolute orientations q 1_G and q 2_G of the two MIMUs were computed separately for every β value (76 values from 0 rad/s to 1.5 rad/s); β = 0 rad/s, the orientation is obtained from the gyroscope only, β >0 the orientation determined from the accelerometer and magnetometer comes into play more as β increases. To compare the orientation of the MIMUs with that obtained with the SP, it was necessary to refer these orientations to a common GCS (the GCS of the SP and MIMU were not coincident on the horizontal plane). To this purpose the accurate alignment of the LCS of each system was exploited: q 1_G , q 2_G , and q S P_G were separately referred to their initial frame to obtain q 1 , q 2 , and q S P , respectively.
The absolute orientation errors q abs1 and q abs2 were computed in the quaternion form as follows: The relative orientation difference between MIMUs pair (q r el ) was computed in the quaternion form as follows: To obtain a compact representation of the error, q abs1 , q abs2 and q r el were converted into their corresponding rotations (θ abs1 , θ abs2 and θ r el , respectively) by inverting (1) to compute the scalar part of their quaternion. Then, the two absolute error values θ abs1 and θ abs2 were averaged to obtain θ abs . Lastly, the rms value of θ abs and θ r el were computed only during the dynamic portions of the recording to obtain e β and δ β . This procedure was repeated for each value of β to obtain the two vectors e and δ which contain the absolute error and relative difference for each of the 18 experimental condition (3 commercial products x 3 angular rates x 2 motion directions).

F. Metrics Used for Testing the Validity of the "Rigid-Constraint" Approach
The approach validity was tested for each of the 18 experimental conditions. The accuracy of the MIMU orientation obtained using the estimated β value (β * ) was assessed against the best available MIMU orientation (β opt ) and the reference orientation computed using the SP system.
The following quantities were used in the validity testing procedure: The procedure implemented for a quantitative error evaluation is explained below referring as an example to the results obtained for a specific dataset (APDM at medium angular rate, 3D): 1. For each condition, the relative differences δ and the absolute errors e are plotted together for each value of β (Fig. 5). 2. The β value (β opt ) corresponding to the lowest absolute error (e min ) is identified. 3. The interval of β values (referred to as optimal interval, β opt ) for which e falls within e min and e min + 0.5 deg was identified (blue-solid horizontal segment in Fig. 6). This threshold was defined considering differences smaller than 0.5 deg as not relevant, being of the same amplitude of the errors affecting the SP estimates (as stated in section II.E). 4. Finally, the β value (β * ) which corresponds to the minimum of the relative differences was identified (greendashed vertical line in Fig. 6).
It is therefore possible to verify if β * is included in β opt and to determine the absolute error e * obtained for β * .
This procedure was repeated for each considered condition. The absolute orientation error of each MIMU is also represented. Fig. 6. In this picture the optimal interval (Δβ opt ) and the β suboptimal value (β * ) were identified with the blue-solid horizontal segment and the green-dashed vertical line, respectively. In this case β * is included in Δβ opt . The absolute orientation error of each MIMU is also represented.
III. RESULTS Noise description for each sensor of each MIMU and gyroscope bias were reported in Table. II and Table. III, respectively (Appendix) [44].
To facilitate the comparison across the different experimental conditions, the values of β opt and β * were represented in Fig. 7.
The values of β opt and β * , along with the relevant errors e min and e * are reported in Table. I. The last column of Table I shows the absolute errors (e M AD ) obtained by using the default value of β = 0.1 rad/s adopted by Madgwick et al. in the original MATLAB implementation of their code. In particular, e * and e M AD are reported in bold when their values fall within the interval e min + 0.5 deg.

IV. DISCUSSION A. Different Optimal Parameter Value for Different Experimental Conditions
In the last decades, great efforts have been dedicated to present different analytic formulations for improving the accuracy of the orientation estimate based on data fusion. In this perspective, several studies have compared the filter performances in different conditions [22], [27], [36] in order to identify the most effective solutions. However, to carry out a fair and general comparative evaluation among methods, it is necessary to preliminarily chose effective values of the sensor fusion algorithms parameters. In general, parameter values are tuned based on trial-and-error by minimizing the overall differences between the gold standard orientation data and those estimated by the sensor fusion algorithm using the data collected during a specific movement by a specific hardware. This approach allows to find optimal parameters values, but it can be only applied to laboratory recordings.
The importance of a proper choice of the values of the sensor fusion algorithms parameters is widely recognized [22], [26], [2] but to the best of the authors' knowledge the literature is still lacking methods for their estimate. The only published paper reports a genetic algorithm to automatically find the optimal weight for a sensor fusion algorithm given the specific dataset and the ground truth values [25].
Our results confirmed that the selection of the optimal value of the sensor fusion parameters depends on the specific experimental conditions considered. In fact, by observing Fig. 7, there is not a common intersection interval among the β opt in all cases. This means that the value of β opt which minimized the absolute orientation errors varied as commercial product and motion characteristics (range of variation) changed. As a consequence, it was impossible to identify an optimal interval of β values suitable for all 18 conditions (3 different commercial product characteristics, 3 angular rates, and 2 different conditions of motion dimensionality). This finding is in accordance with other studies [21], [22], [28] which state that parameter values have a strong influence on the orientation accuracy at different frequency and amplitude of the performed movements. As shown in Fig. 7, under the same angular rate and dimensionality conditions, the β opt ranges highly varied depending on the specific commercial Fig. 7. The absolute error obtained considering Δβ opt , e min , was reported above each segment for each considered condition. When β * values were outside the Δβ opt interval, the amount of the absolute errors obtained considering β * , referred as e * , was reported above β * identified by the black dot. product considered due to the different noise characteristics (see Table. II and Table. III in the appendix). This is expected considering that the final orientation estimate is determined by weighting the information provided by the gyroscope, accelerometer and magnetometer at each time step, and each weight should take the specific sensor noise characteristics into account. This is quite evident by a preliminary comparison of the gyroscope bias changes observed before and after the acquisition for the different MIMUs (e.g. 0.17 deg/s for APDM vs 0.03 deg/s for Xsens and 0.12 deg/s for APDM vs 0.02 deg/s for Xsens at medium and high angular rates, respectively).
In addition to sensor noise characteristics, also the motion characteristics have an influence on the definition of the β opt (Fig. 7). In fact, when increasing the angular rate, the ratio between the body acceleration and gravity component increases and the MIMU inclination, as provided by the accelerometer readings only, is less reliable. However, as it can be noticed in Table. I, the latter statement does not necessarily implies lower values of β opt in order to reduce the accelerometer contribution. In fact, the drift on the gyroscope vertical axis can be only corrected by the magnetometer readings (in absence of ferromagnetic disturbances, as we can consider the experiments performed in this study). In this case a high value of β allows for the magnetometer contribution to influence the resulting orientation. In conclusion, since the mathematical formulation of MAD adopts the same β value for weighting both accelerometer and magnetometer signals, it is very difficult to generalize the results interpretation. A similar consideration can be drawn when comparing motion with different characteristics and involving different rotation axes (e.g. single axis movement versus 3D movements) as the ratio between the body acceleration and gravity component can differ even under the same angular rate, depending on the specific motion.
The abovementioned considerations further contribute to acknowledge the importance of selecting parameter values that enhance filter performance under specific experimental conditions.

B. Validity, Limitations, and How to Use the Rigid-Constraint Approach in Practice
The proposed "rigid-constraint" approach relies on the general assumption that the β value minimizing the relative orientation changes, is also a suboptimal choice for minimizing orientation errors under the same specific experimental conditions. In this respect, it can be noted that the "rigid-constraint" approach resulted in the identification of suboptimal β values in most cases as for 13 different conditions out of 18, the estimated β * values were included in the β opt interval (difference between e * and e min lower than 0.5 deg). In the remaining cases, the error increase with respect to β opt was lower than 2.5 deg. Conversely, when using a predefined default value for the parameter (β = 0.1 rad/s as proposed in Madgwick et al.), a general and consistent error increase between 0.9 deg up to 3.2 deg was observed across conditions. Moreover, the "rigid-constraint" approach was effective in providing the best available orientation with Xsens, whereas considering APDM and Shimmer, the suboptimal β values led to slightly higher errors. All errors obtained refer to the specific commercial product employed in this study; a further reduction of the orientation errors is likely when using hardware with better noise characteristics.
As already mentioned, some limitations must be considered when applying the "rigid-constraint" approach. In fact, when r is null, the assumption that the sources of noise affecting accelerometers and magnetometers are different is no longer valid. In fact, if the relative distance approaches zero the differences between the two sensed accelerations and magnetic fields, are negligible, regardless of the body accelerations (7) and ferromagnetic disturbances (8). In this case, by setting β to a high value (to limit the gyroscope contribution) the relative orientation difference will be small due to the similarity of the measurements, but this does not guarantee for a small absolute error, especially when the accelerometer and magnetometer readings are corrupted by important body accelerations and distorted magnetic fields, respectively. From a practical point of view, it is suggested to position the two MIMUs so that r is sufficiently large (at least a few centimeters) compatibly with the size of the rigid body.
Regarding the applicability of the approach, often in clinical gait analysis and in sport applications, sensors can be firmly attached to the subject by means of mounted rigid plastic plates using elastic straps. In light of the miniaturization progress, it is reasonable to expect that for several applications a specifically designed plastic plate can host two MIMUs during the recording and the proposed method can be applied for filter tuning, also online. An example of this situation is described in [45], where a support was designed to be attached to the foot. Alternatively, when the interest is the analysis of a given motor task such as gait clinical test or when the same sensors are used during the entire data collection, then it would be reasonable to define the suboptimal parameter values on preliminary movement data acquisitions which mimic the actual experimental conditions (i.e. sensor model, motor task, speed, etc.). A similar approach was implemented in the article by Cardarelli et al., [46]. In that case, the authors estimated the orientation of a lower back-mounted MIMU to remove the gravity for estimating the position by means of double-integration technique. This operation leads to an important position drift. To minimize this effect, they used a Weighted Fourier Linear Combiner whose parameters where tuned by minimizing the estimated position difference from the position of the marker mounted on the MIMU. In addition, the authors recognized the importance of changing the parameters according to the scenario under analysis. The main finding was that the weights obtained from a first dataset as training set were useful to improve the accuracy of the position estimates in test set with similar motion conditions.
To summarize, when using the MAD filter, and in general any sensor fusion filter for orientation estimation, it might be not appropriate to select parameters values reported in the original articles as those values usually reflect the specific experimental conditions and hardware employed for data recording. This caution should be exercised especially when different sensors and motor tasks are analyzed. When neither the preliminary ad-hoc experiment nor two MIMUs are available, results reported in Table I can provide a preliminary indication on the expected optimal range for β values when employing the same commercial products and similar type of motion. If even the commercial product employed is different from the three described in this article, it is still possible to refer to the results of Table. I by considering the product whose noise characteristics (listed in Table. II and Table. III) best reflect the one employed.

V. CONCLUSION
The proposed "rigid-constraint" approach was designed to estimate the suboptimal parameter value of the sensor fusion algorithms by exploiting the hypothesis that a set of MIMUs aligned has a null and constant relative orientation over time. The proposed approach does not require the knowledge of any orientation reference. In the majority of cases the approach selected suboptimal values which were included in the optimal intervals, whereas in the remaining the errors were acceptable (maximum difference equal to 2.5 deg and less than 1.5 deg, on average).
Another key finding of this study was the absence of a unique parameter value interval suitable for all experimental conditions. This empirical evidence is in accordance with the previous study of [22] in which the authors hypothesized the existence of a strong relationship between the parameter values and the corresponding orientation accuracy according to the experimental conditions. For this reason, the parameters tuning is crucial for every sensor fusion algorithm to obtain reliable orientation estimates.
The filter from Madgwick et al. is popular, easy to use, open source, and tuned with a single parameter. For this reason, it was chosen in this work as a paradigm to test the effectiveness of the method proposed. However, a single value for that parameter may be a limitation when long acquisitions are performed or when ferromagnetic disturbances are present. Future work will be devoted to test the validity of the proposed approach to filters requiring more than one parameter. In such case the search space would be extended, and the analysis of the most relevant parameters should be performed first to limit the increase of the computational burden. Very preliminary results [47] suggest that errors obtained by rigid-constraint approach with other filters are reasonable. This is justified by the fact that the rigid-constraint approach relies only on the orientation estimated by a filter, regardless of its specific mathematical structure.
In the following Table. II, the noise standard deviation of the sensors of each unit is reported. The evaluation was carried out on 60 seconds of static acquisition.
Firstly, it is possible to note that all the sensors embedded in different commercial products were characterized by different noise value. On the other hand, by comparing the sensor noise values of two units of the same commercial product, it may be observed that the noise values are not exactly the same. The differences, however, are generally lower than those obtained from the comparison of different commercial products [48].  Table. III the three-axial bias of each gyroscope was computed in the static trials before and after each dynamic acquisition to assess the bias changes. The bias at the beginning of the acquisition and the difference () with respect to the end were listed in the Table. III.
It is possible to assess that, in the majority of cases, the biases are not constant between the beginning and the end of the trials. In particular, the greatest differences arose for each trial of APDM with values of up to 0.17 deg/s after less than one minute. It has to be said that these measures of bias before and after the dynamic trials are not meant to be a proper characterization of the bias instability (which can be computed through the Allan deviation over a long time acquisition [49], [50]), but they can give an overview of the severity of the bias changes within the same recording. As it can be seen from the several attempts published over the years to tackle the online bias estimation (e.g. [9], [13]- [15]), a varying bias is the most problematic issue affecting the orientation estimates obtained using the gyroscope since these changes are difficult to be predicted and modelled.
Finally, the specifications of the calibrated data for each commercial product used are listed in the Table. IV.