{'fraction_non_alphanumeric': 0.028344999132295944, 'fraction_numerical': 0.00578469370046856, 'mean_word_length': 4.651356652500818, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Efficient Simultaneous Calibration of a Magnetometer and an Accelerometer

This paper describes a calibration algorithm to simultaneously calibrate a magnetometer and an accelerometer without any information besides the sensors readings. Using a linear sensor model and maximum likelihood cost, the algorithm is able to estimate both sensors' biases, gains, and covariances, besides sensor orientations and Earth's fields. Results show errors of less than 0.1 standard deviations in simulation, and high-quality estimates with real sensors even when the algorithm's assumptions are violated.


I. INTRODUCTION
S ENSOR reading is an important step in a robot's control loop. However, the values obtained may be inaccurate because of incorrect sensor calibration. In the field of aerial vehicles, accelerometers and magnetometers are frequently used to estimate attitude [1]. As both these types of sensors are calibrated using Earth's gravitational or magnetic field as reference, most research into calibration of these devices focuses on magnetometer calibration, which can also be used to calibrate accelerometers in stationary settings. This work focuses on batch algorithms, that is, algorithms that collect a data set and then calibrate the sensors, which can be used to provide the initial conditions for online methods that solve the calibration problem by filtering [2] or iterating [3], [4].
A common approach to calibration is to assume that the magnetic field's norm is known or unitary, and then adjust the gain and bias so that the readings match this value. [5] extended the TWOSTEP algorithm, which tries to minimize the difference between the field's norm and the sensed norm, to calibrate both the bias and gain. [6] proved that noiseless magnetometer readings, including soft-and hard-iron effects, occur on the surface of an ellipsoid manifold, and then developed an algorithm to calibrate the bias and gain.
Although both methods can be used to calibrate accelerometers as well as magnetometers, the calibration must be performed independently, so the magnetometer and accelerometer are calibrated in different frames, and the rotation between them is unknown. Furthermore, neither method is able to estimate the direction of the magnetic field, which may be unknown and is essential for attitude estimation using magnetometers [2], [7], [8]. [9] partially solves this issue by calibrating the magnetometer while also capturing a reference for the z direction, provided by some independent filter. This allows the direction of the Earth's magnetic field to be estimated, and the calibration to be performed in some previously known reference frame, while also estimating the magnetometer's bias and gain. However, a reference orientation, whether reliable or not, is not always available, thus limiting the use of algorithm. [10] tried to solve both problems by simultaneously calibrating the magnetometer and accelerometer, allowing the magnetic field, the attitude, and all the sensors' parameters to be estimated. This eliminates the need for a reference so that the calibration can be performed from sensor readings alone, that is, no external apparatus or knowledge besides the sensor readings is required. However, the problem is replete with lowquality local minima, which sometimes significantly degrade the calibration performance, and the solution has a very high computational cost, making the algorithm impractically slow.
The novelties in this paper are the improvement over [10] made possible by the use of better initial estimates of the parameters, leading to better minima, and the approximation for the rotation estimation, which also helps to avoid poor minima and reduce computational time because of its closedform solution, making it suitable for use in practice. The initial conditions are given by an adaptation of [9], where the z reference is replaced by an estimate computed from the accelerometer. To our knowledge, the algorithms described in the present work and in [10] are the only accelerometer and magnetometer calibration algorithms that are able to perform the complete calibration from collected data alone. Moreover, the algorithm proposed here is also the most statistically correct and generic solution based on linear models since all the parameters (gains, biases, Earth's fields, and orientations) must be estimated at the same time.
The algorithm is tested in simulated and real sensors. The simulations show that the algorithm correctly calibrates different sets of sensor parameters and allows comparison between the estimates and the ground truth. Several variations of the algorithm are tested, and the best-performing one is also used to calibrate the real sensors. The data capture was designed to violate the basic assumptions underlying the algorithm so that its robustness to adverse real-world conditions could be tested. It is shown that even in this imperfect setting the proposed algorithm is able to achieve high-quality calibration.
The paper is organized as follows. Sections II and III describe the sensor model used and the collection and preprocessing of the data, respectively. Based on the sensor model, a maximum likelihood cost function is described in Section IV. The initial parameters' estimates are given in Section V, and the optimization steps are described in Section VI. The results of simulated and real-world experiments are reported in Section VII and are used to evaluate the proposed method.
The concluding remarks are presented in Section VIII.

II. SENSOR MODEL
The linear sensor model with Gaussian noise, where the correct value is scaled by a matrix and added to a bias and a noise, is frequently used in the literature [5], [6] and is also used in this paper. The value v s ∈ R 3 read by the sensor s, where s = a and s = m for the accelerometer and magnetometer, respectively, is given by: where K s ∈ R 3×3 and b s ∈ R 3 are the gain matrix and bias vector associated with the sensor, v * s ∈ R 3 is the nominal value to be read, and Σ s ∈ S 3 ++ is the noise covariance, which is a positive-definite matrix. It is important to note that the linear sensor model for the magnetometer is also able to handle softand hard-iron effects, which can be embedded in the sensor's gain and bias [6].
Conversely, given the parameters and a value read by the sensor, the nominal value can be estimated as: where it is assumed that the gain matrix K s is invertible. If it is not, then its columns are not linear independent and the measured values v s must lie in a reduced subspace spanned by these columns. Therefore, only the case where K s is invertible is considered, and the nominal measurement v * s can be recovered.
The nominal value for a sensor is given by the Earth's respective field rotated by the current orientation of the sensor frame, so that v * s = Rv I s , where v I s is the field value for a given sensor in the fixed inertial frame I and R ∈ SO(3) is the rotation representing the current sensor orientation. Since the choice of I is arbitrary, one can choose it such that the gravity g only has a negative component in z and the magnetic field h has a positive component in x, such that It should be noted that the assumption that the magnetic field component h x is greater than zero can be satisfied almost everywhere except at the magnetic poles, where the magnetic field is perpendicular to the Earth's surface and the magnetometer does not provide more orientation information than the accelerometer does for the gravitational field. In this particular case, one cannot use this combination to perform the calibration, but one can do so in all other contexts, as will be shown.
As discussed in [10] and highlighted later, the sensor frame also has multiple possible definitions, since the gain matrix can be written as K s = RK ′ s , where R is a rotation matrix. Therefore, some constraints must be fixed to make the solution unique. This paper assumes that K a is triangular, which is a sufficient assumption to define a single solution. Moreover, K m is allowed to have any value, so that it can incorporate the misalignment between the accelerometer and magnetometer.
Since during the calibration the algorithm does not use information about the body frame, and this may not be available if the calibration is being performed outside the desired body, other algorithms such as [11] can be used to align the sensor and body frame.

III. DATA PREPROCESSING
Assuming the linear model introduced in Sec. II and that the sensors can be held relatively still when readings are being taken so that a set of sensor measurements have the same nominal value, the measured values can be written as: where N is the number of measurement sets, ∆ s [i] is the number of data points collected for the i-th set, µ s [i] is the mean value shared by the measurements in the i-th set, and R i is the orientation at which the i-th set was collected.
Since the measurements in each set are assumed to have sampled the same information, the mean values µ s [i], which are different for each set, and the covariance matrix Σ s , which is shared among all sets, can be estimated.
The minimum variance unbiased estimator for µ s [i] is the sample mean [12], given by: and has a distribution described by: Since the sample meanμ s [i] characterizes the values for the i-th set, it can be used instead of them as a single measurement with appropriate weight, associated with the number of measurements ∆ s [i], during optimization to decrease the computing time.
The minimum variance unbiased estimator for the noise covariance is given by: where the subtraction of N from the denominator is a generalization of the Bessel's correction [13] for N measurement sets, instead of 1 as in the original correction.
IV. COST FUNCTION Given that it is assumed that only the measured values are known, the full parameter set is given by all the sensors' and fields' parameters and all orientations. This set can be written as: Note that this is the most generic setting possible for this problem with the assumption of affine measurements, since all parameters are considered unknown and any previously known parameter can be fixed to the known value.
Using the model defined in Eq. (3), the negative loglikelihood can be used to create a cost function [14], which is given by: In this format, the calibration problem becomes a traditional minimization problem over all the parameters in Θ, where the restrictions over Σ s , g z , h x , and R i must be satisfied. A simplified cost function can be defined on a smaller parameter set Θ ′ = Θ \ {Σ a , Σ m }, where the covariances are replaced by their estimates computed using Eq. (6). In this case, the maximum likelihood cost can be simplified to , which greatly reduces the computational cost.
Since both cost functions are subject to many local minima, which can lead to poor estimates as shown in [10], a good initial estimate must be provided. Then, instead of optimizing all the parameters at once, it will be shown that optimizing subsets one at a time greatly simplifies the problem. The next two sections deal with these problems.

V. INITIAL ESTIMATE
This section is an adaptation of the work presented in [9], where a magnetometer was calibrated using a known correct reference measure of the z direction. Since it is assumed that no such z is available, this reference is replaced here by the estimated accelerometer's field value, which only has a negative component in the z direction.
Furthermore, the constraint that the magnetometer cannot have mirrored axis relative to the inertial axis is relaxed, which was not possible in the original reference. The reader is referred to the original paper for a more detailed derivation of the equations and implementation details.
Assuming that the field v I s is unitary and that the noise ǫ s can be neglected, the following approximation can be performed: which follows directly from taking the norm of both sides of Eq. (1) and replacing the single measurements by their mean given by Eq. (4), which has smaller variance and thus is a better estimate.
By separating the known terms in Eq. (9a) from the unknown terms and considering all N sample sets, the approximation can be written as: where ⊗ denotes the Kronecker product and vec denotes the vectorization operator. Assuming that A is symmetric, the nontrivial solutionη that minimizes the approximation error is given by the right eigenvector corresponding to the smallest singular value of J [9]. Since any αη is also a solution to the approximation in Eq. (10), one can find the correct α by manipulating the elements in Eq. (9), which provides the correct scale. Hence the correct value is given by: From the definitions in Eq. (9) and the estimatesη and α, the gain and bias must satisfŷ which gives the bias estimate for each sensor directly. The gainK s can not be uniquely determined, as anyK s R with RR T = I 3 is also a solution to Eq. (11). If the upper triangular solution obtained by the Cholesky decomposition is denoted byK s , the estimates can be written as: where it is assumed that onlyK m is subject to rotation as the rotation between the two gains is fixed and there would be multiple solutions if both had full unknown rotations. So far, the solution described in [9], which provides estimates for the sensor bias and gain, has been used to compute the initial estimates. This solution could be applied here since the accelerometer and magnetometer could be considered decoupled problems. However, in Eq. (12),K a can be estimated directly whileK m depends on the unknown rotation R. Since only the magnetometer was calibrated in [9], this was not an issue as the ground-truth z direction was available. Moreover, the magnetic field component h z could also be estimated directly. To solve this problem, the ground truth z was replaced by the estimate of the nominal value v * a from the accelerometer, since this is already calibrated. Therefore, estimates for the rotation R between the accelerometer and the magnetometer and the component h z of the magnetic field can be found simultaneously by solving the following optimization problem: where ∆ a and ∆ m are used to weight the number of samples in each data set. As stated earlier, this equation is similar to the one in [9], but the ground truth z is replaced by the estimates Test showed that, unlike in [9], this problem is subject to a few low-quality local minima because of the replacement of precise knowledge of the z direction by the accelerometer estimate and that the initial condition R = I 3 and h z = 0 sometimes leads to poor estimates. Hence, an optimization with random restarts was performed and the initial conditions were given by a random rotation as described in [15] and h z ∼ U([−1, 1]). Tests showed that at most 100 iterations were needed to find good initial estimates.
Since it was assumed that the fields were unitary, it follows thatĥ As discussed in [10], having the fields' and gains' scales as variables leads to multiple solutions, since the gain can be multiplied by a constant while the field is divided by the same value. Therefore, in this paper the value g z = −1 and h x = 1 are adopted, and the magnetometer estimates must be adjusted as:K ′ m =K mĥx ,ĥ ′ z =ĥ z /ĥ x , while the accelerometer does not have to be adjusted because it already has a unitary field.
Finally, the rotation R i associated with each sample set is given by the eigenvector corresponding to the minimal eigenvalue of B i [16], where B i is defined as: and the weights w a and w m introduce uncertainty weighting. All initial estimates therefore have closed-form efficient solutions except for R and h z . However, Eq. (13) has few decision variables and can be efficiently optimized by random restarts and gradient descent.

VI. OPTIMIZATION STEPS
As the cost functions in Eqs. (7) and (8) have many joint parameters, the parameter set is split into four disjoint subsets that can be solved efficiently one at a time. The optimization starts with the estimates provided in Sec. V and then iterates each subset optimization while keeping the other parameters fixed until convergence is achieved. Convergence is defined by J k−1 − J k < γ, where J k is the cost after the k-th iteration and γ is the stop parameter. Since the true values for any parameter in Θ are not available, the hat notation is dropped where convenient to simplify the expressions.
This section is an extension of the optimization presented in [10]. The main differences are, firstly, the approximation for the rotation estimation described in Sec. VI-A. This uses Eq. (14) to compute an efficient approximation to the real orientation during the optimization, which is not present in [10], leading to a significant slow-down of the algorithm and worse results. Secondly, the format of the gain matrix K a computed in Sec. VI-C is assumed to be triangular here because its initial condition is given by a Cholesky decomposition while in [10] it is assumed to be symmetric. The main difference between the method proposed here and the one described in [9], apart from the fact that two sensors are used and no ground-truth is available in the former, is that the optimization is performed in steps with efficient solutions, while in [9] the initial estimate is computed and then all parameters are optimized at the same time, preventing the inherent structure of the problem from being exploited.

A. Rotation Estimation
From Eq. (8), it is clear that the cost for the i-th rotation is given by: . As a closed-form solution could not be found, the problem is optimized directly through gradient descent using a quaternion to represent the rotation. However, the approximation used to compute the initial estimates, which is given by the eigenvector associated with the minimum eigenvalue of B i from Eq. (14), can be used to speed-up the initial steps of the optimization, when the estimates are still far from the true values. The performance of this approximation approximation is evaluated in the experimental section.

B. Bias and Field Estimation
Since the sensor readings are assumed linear in the biases and fields and have Gaussian noise, the problem can be posed as a generalized least squares (GLS) problem [17] given by finding β to approximate y = Xβ + ǫ, E[ǫ|X] = 0, Var[ǫ|X] = Ω, whose solution is given by: Writingμ a [i] in the GLS format using Eqs. (3) and (5) gives   μ a [1] . . .
where z is the unitary vector in the z direction. As Sec. V established that g z = −1, the gain matrix is adjusted to ensure that this equality is still valid. The estimate is therefore adjusted so thatK ′ a = −K aĝz , whereĝ z is the estimated gravity.
The magnetometer estimateμ m [i] can be written in a similar fashion, allowingĥ x ,ĥ z andb m to be estimated. As h x = 1, it follows thatK ′ m =K mĥx andĥ ′ z =ĥ z /ĥ x , which maintains the scale.

C. Gain Estimation
The gains can also be written in the GLS format using Eqs. (3) and (5). Assuming that K a is upper triangular, which is valid for the initial estimate since it comes from a Cholesky decomposition, one has that: where (v) i represents the i-th component of the vector v. Again, a similar formulation can be used for the magnetometer, but the gain K m is considered full, as it is the product of a rotation with an upper triangular matrix.

D. Covariance Estimation
The maximum likelihood covariance estimator is similar that given in Eq. (6), but the referenceμ s [i] is replaced by an estimate of µ s [i] using parameters in Θ, and Bessel's correction is not used, since it lowers the likelihood of the data in order to provide an unbiased estimate. The new estimate is then given by:

VII. EXPERIMENTS
Two experiments are performed to validate the proposed algorithm: calibration of simulated sensors and calibration of real sensors. The simulation is used to evaluate the effects of changing the stop condition, the number of sampled intervals, and the algorithm hypothesis. The performance is analyzed in terms of errors in relation to the ground-truth values, which are available in the simulation, and the total time required for the simulation. Using a Monte Carlo approach, many combinations of sensors can be simulated to evaluate the robustness of the proposed algorithm.
The performance with test data, which is not used for training, can be evaluated to test the generalization of the algorithm [14].
In the real experiment, a magnetometer and an accelerometer had their readings measured while being held by hand inside a building, which slightly violates the hypothesis that all samples measure the same value. It will be shown that the algorithm is able to fit the sampled data correctly.
It is important to note that only the algorithm derived in this paper is used for calibration since the authors do not know of any other work, apart from their own previous contribution [10] upon which this paper is based, that describes an algorithm that can perform this kind of calibration without external references. If an algorithm such as the one proposed in [9] were used, additional information would have to be provided, making comparison of the results impossible since the calibration would be performed on different data. Thus, no comparison with any previous work is performed.

A. Simulated calibration
To evaluate the calibration algorithm, 100 Monte Carlo simulations with new sensor parameters on each run were performed, allowing evaluation of both common and possible exceptional calibration behavior. This analysis is important to show that the method described handles different sensor characteristics well. Table I  The rotations for each set were randomly created using the method described in [15], and the stop condition was set to γ = 10 −4 .
When varying the stop condition γ, the same sample sets were used for each value tested, so that the performance can be viewed as the same algorithm stopping at different times. This avoids differences caused by random fluctuations and allows better comparison. For a varying number of sets, runs with N + 1 sample sets merged a new set to the previous N already used, simulating collection of an increasing number of sets and also preventing random fluctuations from significantly affecting the comparison. Alternatively, the training with N sets can be viewed as training on a subset of the data set and the covariances are given by where α s = 10 ǫ , ǫ ∼ U([−2, −4]), which guarantees that the covariance matrices are positive definite, but allows high condition numbers. Instead of using the nominal value K m , the magnetometer is rotated by R ′ and mirrored by M , so that the resulting gain matrix is given by K ′ m = M R ′ K m , showing that the algorithm described is able to handle any configuration between the accelerometer and magnetometer, unlike [9]. The rotation R ′ is chosen randomly [15], and M = diag(δ 1 , δ 2 , δ 3 ), δ i ∼ U ({−1, 1}), is a random diagonal matrix that indicates the mirroring of each component.
After calibration, the estimated value of Θ was used to rebuild estimatesμ s [i] for µ s [i], as in Eq. (3). To compare the quality of the reconstructions, the following measurement was used:  where Σ s is the real sensor covariance. This measurement is derived from the Mahalanobis distance [18] and represents the average error in standard deviations.
The rotation optimization described in Sec. VI-A can be performed in one of two ways: always using gradient descent in Eq. (15) or using the approximation in Eq. (14) until J k > J k−1 and then changing to gradient descent in subsequent iterations. As the approximation does not minimize the cost directly, eventually the change to gradient descent will be made. Another variation involves re-estimating the covariance using Eq. (16) or only using the initial estimate given by Eq. (6). As some approaches consider the covariance matrix to be diagonal [6], this restriction is also evaluated.
Hence there are four algorithms to compare: not refitting covariance and optimizing the rotation directly using gradient descent (NCDR); not refitting covariance and approximating the rotation (NCAR); and fitting full covariance and diagonal covariance and approximating the rotation (FCAR and DCAR, respectively). The results obtained while readjusting the covariance with direct rotation optimization are not reported as these were considerably worse.
For each sample set, another set with different samples and rotations for the same sensor parameters was created to evaluate the generalization of the fitted parameters, that is, if the parameters are able to fit the data not present during calibration well. These new samples follow the description in Sec. III and provide new valuesμ ′ [i], which are used to compute the new rotations. The rotations were first approximated and then optimized directly to reduce the cost function, as described in Sec. VI-A. Figure 1 shows the reconstruction error for different stop conditions. For γ < 10 −4 , the error does not change signifi- cantly, justifying the choice of this value as a reference when evaluating different number of samples. Except for NCDR, all the algorithms have similar behavior for lower values of γ, while adjusting the covariance was beneficial for larger stop conditions. However, as larger stop conditions can only be justified if the computational cost is very high, it will be shown that this difference is not relevant. NCDR has slightly worse magnetometer calibration because of the failure of the gradient descent method to escape local minima.
It is interesting to note the apparent magnetometer overfitting for γ < 10 −2 , as the reconstruction error for the training set increases even though the cost function decreases. Nonetheless, Fig. 2 shows that the reconstruction error for the new data not used during training generally decreases. This suggests that, although the reconstruction worsens, the parameters improve the generalization. The significantly higher error for one case with NCDR indicates that it is not as robust as the other methods using approximations. For all the algorithms, about 75% of all cases have errors of less than 0.1 standard deviations from the correct value for both sensors, showing that the proposed method is able to fit the parameters well. Figure 3 shows that all the algorithms eventually converge to solutions that are close to each other as the number of samples increases, with errors on every run lower than 0.1 standard deviations for N > 25 for all methods, except for one run with NCDR. However, a diminishing return on the number of samples is observed for N > 20 as the error decreases very little if more data is collected. The NCAR algorithm appears to be the most robust, with a lower error spread for a smaller number of samples, so that less user data is needed to achieve acceptable performance for most situations. The NCDR algorithm usually has the largest errors, which, as in the previous case, are the result of poor local minima. Both algorithms with readjusted covariance have similar performance, indicating that the diagonal simplification is a good one. Figure 4 shows slightly better performance overall, indicating that the estimates obtained are able to generalize and fit new unseen samples well.
While the NCAR, FCAR and DCAR algorithms had similar performances for reasonable stop conditions and numbers of samples, the NCDR algorithm not only had the worst performance for all the conditions simulated but also required more computational time, as shown in Fig. 5a. It should be noted that, although DCAR has a simpler covariance to fit than FCAR, it has a higher computational cost because of implementation details. Hence, the calibration time should be about the same for these algorithms in a final implementation. As stated earlier, using higher values of γ and covariance reestimation would be justified if the computational cost was high. However, the cost for NCAR is significantly lower, with γ = 10 −2 having the same cost as γ = 10 for FCAR. Therefore re-estimating the covariance using the rotation approximation does not improve performance and instead slows the processing considerably.
As shown in Figs. 1 and 2, γ values of less than 10 −4 do not improve the performance significantly but cause the computing time for the approximation methods to increase. This occurs because rotation optimization is still taking place, which, although having little effect on the calibration, is time consuming. This also explains why the time taken by the NCDR algorithm stabilizes.
The time for NCAR is almost constant for all the numbers of samples tested, as shown in Fig. 5b. This occurs despite the fact that the cost function in Eq. (8) is a sum that increases as the number of sample sets increases and the stop condition compares two iterations of the algorithm. As more terms are used in the sum, the difference J k−1 − J k would also be expected to increase, so that it takes longer to satisfy the stop condition. Moreover, the optimization steps also increase in complexity, requiring more time for each step in Sec. VI. Nonetheless, the presence of more samples increases the speed at which the estimates converge, requiring less iterations to achieve smaller values of J.

B. Real calibration
Once it has been shown that the proposed algorithm works correctly with different simulated sensors, calibration of real sensors must be tested. The main difference between the real and simulation experiments is that, although the simulation allowed a wide variety of sensor parameters to be tested, it also simulated a perfectly still system during the capture of each dataset, which was an assumption made in Sec. III. Since the NCAR algorithm showed the best performance across all simulations, it was chosen to calibrate the real system.
In real-world applications the system may not be perfectly still, and the robustness of the algorithm in applications where this is the case must be tested. Another assumption was that the magnetic field is constant throughout the measurements, since it is assumed that it can be written as in Eq. (2). This also may not be valid in the real world because of the influence of power cables and any metal structures nearby.
To test the algorithm's robustness under these conditions, a RoboVero board was held by hand to introduce vibration and measurements were taken inside a large room so that magnetic disturbances were present but not overwhelming. Both the accelerometer and the magnetometer were configured to have the same sample rate so that both sensors took the same number of measurements for each set. The set length was determined by using time slices in which the norms of the values read were kept approximately constant. This also introduces errors as the underlying values can vary without changing the norm considerably.
The results are shown in Fig. 6, where the darker lines are the estimated mean sensor readings for each measurement set. These are computed from the estimated parameters using Eq. (3), with length equal to the number of measurements ∆ s [i] in each set. The light lines are the real sampled values. Fig. 6a shows that the estimated values for the accelerometer were very close to the real ones despite the sensors not being held still during each set. This confirms the robustness of the algorithm when the first assumption is violated and supports its use when the sensors cannot be kept stationary. Fig. 6b shows that the magnetometer is also calibrated correctly but that its measurements vary a lot more than the accelerometer measurements. As a result, some measurement sets have slight offsets, as is the case of the sets around 20000 and 26000 samples. These errors are associated with distortion of the magnetic field by the building's metallic structure. However, the small size of these offsets show that this does not affect significantly the calibration, confirming the algorithm's robustness in the presence of this kind of distortion.
It will be recalled that the calibration was performed using only the measured values shown in Fig. 6 without any external references or sensors. Therefore, the proposed algorithm's robustness to violation of the assumptions about the only measurements available, together with the high-quality calibration that it provides, should make it a good candidate for any calibration of these kinds of sensors.

VIII. CONCLUSION
This paper has described an algorithm to calibrate the parameters of an accelerometer and a magnetometer using only sensor measurements without any external information. The algorithm is designed to estimate the gain, bias, and covariance of each sensor, as well as the orientation of each measurement and the direction of the Earth's magnetic field. Hence, the calibration can be performed almost anywhere by anyone as only the sensor readings are required. This is the most generic setting for this kind of problem, since all the parameters are calibrated and any external information can be considered a constraint on the parameters and does not make the calibration harder.
A comparison was made of the base method and its variants, which use an approximation of the cost function being minimized. Simulated results show that the simplest variant is also the fastest and most robust, with the smallest worst-case errors. All the algorithms are able to achieve an error of less than 0.1 standard deviations for the data used to calibrate the devices, indicating that the algorithms were correctly trained, and also for new data not used in the calibration, indicating that the parameters obtained are suitable for use when taking sensor readings and that the learning algorithm was able to generalize. Test with real sensors and adverse conditions that violate the assumptions underlying the derivation of the algorithms showed that high-quality calibration can be achieved in noncontrolled settings, confirming the robustness of the algorithm.
Future studies should investigate use of the estimated sensor rotations to compute the rotation between body and sensor frame, as described in [11], since the sensor frame used was specific to this study and may be different from the desired frame. Another possibility for further study is to integrate external references, such as vision systems [19], so that other sources of information can be used during calibration to reduce the errors even further.