Magnetometer calibration: hard-iron and soft-iron correction by ellipsoid fitting

A rotated three-axis magnetometer traces a shifted, tilted ellipsoid instead of a sphere, because its output mixes sensor imperfections with permanent and induced fields of the host. This note derives hard-iron and soft-iron estimation by ellipsoid fitting, the orthogonal ambiguity that magnetometer data alone cannot resolve, the effect of data coverage, and the mapping of residuals into heading error, then covers temperature, board currents and validation.

Measurement model and the physical origin of each term

A sensor rotated in a constant local field m (navigation frame) reads, at time t, according to the linear model shared across the literature [1], [7], [9]:

R_t rotates from navigation to body frame; ε_t is wideband noise.

• Hard-iron offset (b_hi): field of permanently magnetised parts (speakers, motors); it rotates with the sensor, so it is a constant body-frame vector. DC trace currents add to it [2].
• Zero-field offset (b_zb): inseparable from hard iron, so vendors often leave it uncalibrated [1], [3]. One MEMS data sheet gives ±40 μT typical at 25 °C [15], the order of the geomagnetic field itself.
• Soft-iron matrix (C_si): field induced in unmagnetised ferromagnetic parts (shields, batteries); linear only without hysteresis [1], [10].
• Scale, non-orthogonality, misalignment (C_sc, C_no, R_mb): gain differences, non-perpendicular axes [7], [8] and rotation relative to the inertial axes [9].

Only the lumped W and b are estimable; all linear time-invariant distortions collapse into one affine map [8]. The assumptions most often violated are constant W, b (rigidly attached disturbers) and constant m during calibration [9].

The ellipsoid constraint and the general quadric

Without noise, W⁻¹(h − b) = R_tm; rotations preserve norm, so at every attitude:

A is symmetric positive definite, so (2) is an ellipsoid centred at b with principal axes along the eigenvectors of A and semi-axes F/√λ_i [1]. Expanding yields the general quadric:

s > 0 is an arbitrary scale; the coefficients are defined only up to a factor. Recovering the calibration is direct [7], [9]:

The second follows from b^TMb − k = F²/s. That is 3 + 6 + 1 = 10 parameters [1]. Heading needs only F = 1 [9]; absolute work takes F from a geomagnetic model [7].

Fitting methods

Algebraic least squares

Each sample contributes a row d_i = [x², y², z², 2yz, 2xz, 2xy, 2x, 2y, 2z, 1], and the sum of squared algebraic distances is minimised over the 10-element coefficient vector p:

The constraint excludes p = 0. With C = I the solution is the right singular vector of D for its smallest singular value: cheap, but M need not be positive definite, and with little rotation the result may not be an ellipsoid [9]. Squared raw readings make D ill-conditioned; centre and scale the data first.

Ellipsoid-specific constraint (Li and Griffiths)

Li and Griffiths gave a sufficient condition on the quadratic coefficients p = [a, b, c, f, g, h, …] (their notation) [4], [5]:

With free scale this becomes the normalisation 4J − I² = 1, a quadratic constraint C in (5), and the problem reduces to a generalised eigenproblem. The result is the algebraic least-squares fit among ellipsoids whose short radius is at least half the long one, with a search step for other data [4]; above axis ratio 2 it is biased [5]. Magnetometer ellipsoids are usually near-spherical, but strong soft iron [2] warrants a check.

Bias of the algebraic fit

For a noisy sample h̃ = h + ε, ε ~ N(0, σ²I), the expected left-hand side of (2) is [7]:

The cross term is zero-mean, the last term positive: residuals are not zero-mean and the estimator is biased. More samples do not help; the algebraic fit is inconsistent, as is the geometric fit, with smaller asymptotic bias [6]. Adjusted least squares, which also estimates the noise variance, is consistent and has been applied to magnetometers [7]. With A ≈ I the relative error in the squared radius is 3σ²/F² = 2×10⁻⁴ for σ = 0.4 μT, F = 48.2 μT (derived): minor under full coverage, but with band-limited data it migrates into the centre estimate.

Geometric fitting and maximum likelihood

With isotropic noise in the sensor frame, the maximum-likelihood (ML) estimator minimises orthogonal distances to the ellipsoid. Vasconcelos et al. show that a sub-optimal formulation, with the noise defined in the sphere space, can be searched over the parameters alone [8]:

T is a general 3×3 matrix, solved iteratively (Gauss–Newton, Newton) from an algebraic or two-step start [8], [10]; geometric fitting costs more and depends on initialisation [5]. Kok and Schön pose calibration as ML orientation estimation, optimising the likelihood of EKF prediction errors jointly over magnetometer parameters, gyroscope bias and noise covariances, initialised by ellipsoid fitting and alignment [9].

Reduced models: 4 and 7 parameters

If soft iron is weak, set W = I. Expanding (2) gives h^Th = 2b^Th + (F² − b^Tb), which is linear in the unknowns [1]:

Rows of X are [h_x, h_y, h_z, 1], Y holds ‖h_i‖²; X^TX and X^TY accumulate in constant memory [1]. The 7-parameter model adds a diagonal gain. Smaller models cut variance but leave unmodelled soft iron as bias — a deliberate trade when coverage is limited.

Obtaining the correction matrix

Two common choices of W⁻¹ satisfy W^−TW⁻¹ = A:

The symmetric square root shrinks the ellipsoid onto a sphere along its principal axes without adding rotation [1], [7]. Kok and Schön start from a Cholesky factor and resolve the remaining rotation separately [9]. The two differ only by an orthogonal matrix: W⁻¹_ch = QW⁻¹_sym.

Identifiability: the orthogonal ambiguity and alignment

Since ‖UW⁻¹(h − b)‖ = ‖W⁻¹(h − b)‖ for any orthogonal U, any estimator that uses magnetometer data alone determines W⁻¹ only up to a left orthogonal factor:

The nine elements of W appear as only six coefficients [1]; in W = R_LS_LV_L the alignment V_L is unobservable [8]; the rotation of the calibrated sphere is unknown [9].

What the symmetric choice does and does not fix

By polar decomposition the true correction is W⁻¹ = QP, with P = A^1/2 symmetric positive definite and Q orthogonal. The symmetric estimate returns P, so the corrected vector is

— the true field, still rotated by Q^T. Q = I only if the true W⁻¹ is symmetric positive definite, which misalignment and non-orthogonality generally are not. The NXP claim that symmetry prevents spurious rotation [1] holds only in that no rotation beyond the principal axes is added; physical misalignment is not recovered.

A residual rotation θ_z about the vertical is a constant heading offset. A rotation θ_h about a horizontal axis leaks Zθ_h of the vertical component Z = F sin I into the horizontal plane; with H = F cos I the heading error is sinusoidal with peak tan I · θ_h (derived), like tilt error [3]. At I ≈ 58°, tan I ≈ 1.6: 1° of unobservable horizontal rotation yields about 1.6° of heading error behind a flawless fit.

Resolving the ambiguity

Kok and Schön use the rotation invariance of inner products: the vertical v_t from accelerometer and gyroscope, dotted with the calibrated field, must equal the constant vertical component m_z [9]:

Varied attitudes determine R_D; joint ML then refines everything. Over 24 orientations of a right-angled block, mean heading error fell from 1.28° (initial) to 0.76° (ML), maximum from 4.36° to 2.48° [9]. Vasconcelos et al. align against external reference vectors as an orthogonal Procrustes problem: two independent vectors fix SO(3), three are needed for O(3) [8].

Data coverage, noise and conditioning

The fit's information is the curvature the samples span; samples at one attitude differ only by noise, so roll and pitch must vary clearly [1]. Gebre-Egziabher et al. showed that a non-aerobatic aircraft traces only a narrow strip; with a 10° strip and high noise their solver diverged about as often as it converged, a two-step initialiser removed divergence, and z-offset uncertainty was 3–5 times that of x [10]. Vasconcelos et al. report good ML accuracy on ring data (free yaw, ±20° pitch), degrading somewhat as curvature shrinks [8].

A simple Monte Carlo example separates these effects (derived): F = 48.19 μT, I = 58.19° (Ankara, WMM2025 [12]); symmetric soft iron with diagonal (1.08, 0.95, 1.02), off-diagonals 0.015–0.035; b = (25, −12, 40) μT; white noise σ per axis; 300 samples; 10-parameter algebraic fit on normalised data, symmetric root. “±φ”: free yaw, roll and pitch uniform within ±φ. Heading error uses ideal tilt compensation at 72 headings, “level” at zero tilt, “tilted” at roll/pitch ∈ {−40°, 0°, +40°}; medians of 100 trials.

Table 1. Effect of coverage and noise (illustrative simulation, 10-parameter algebraic fit)

At ±20° the centre is off by about 15 μT, yet level heading error is only 0.65°: the fit is right where the data are and wrong elsewhere, exceeding 20° at 40° tilt. With 3000 and 10 000 samples the centre error did not shrink; bias (7) and narrow-band conditioning are not cured by sample count. Noise and coverage interact: 2 μT is tolerable with full coverage, but at ±60° the centre error reaches 3.6 μT. A 1° rotation added to the true W gave maximum level heading errors of 1.67° about a horizontal (x) axis and 1.17° about the vertical (z), matching tan I ≈ 1.61 and the constant-offset prediction, with residuals at noise level.

Heading error budget

After tilt compensation, for the horizontal field h = H[cos ψ, sin ψ] and a small perturbation δh, δψ = (h × δh)/H²; amplitudes are independent of sign convention.

Residual hard iron. With δh = Δb:

This is the first harmonic, matching the NXP result [3]. Rule of thumb: at H = 10 μT, 3° (0.05 rad) requires field error below 0.5 μT [2].

Residual soft iron and gain. With residual matrix I + E, δh = Eh:

Derivation: apply the double-angle identities to h × Eh = H²(E₁₀cos²ψ + (E₁₁ − E₀₀) sin ψ cos ψ − E₀₁sin²ψ). The symmetric part gives a second harmonic, equal at ψ and ψ + π as NXP notes [3]; the antisymmetric part is a constant — the vertical-axis rotation the fit cannot see. (14) and (15) reproduce the classical compass-swinging model [10]:

Tilt and horizontal-axis rotation. For a tilt error δθ, or a residual rotation about a horizontal axis, δψ_max ≈ tan I · δθ [3]; the small-angle form breaks down towards the magnetic poles.

Noise. With noise component σ_⊥ perpendicular to the field in the horizontal plane, σ_ψ ≈ σ_⊥/H; averaging N independent samples reduces it by √N [3].

Reference model. The WMM2025 1σ uncertainty of declination, needed for true north, is √(0.26² + (5417/H)²) degrees with H in nT; regions with H < 2000 nT are blackout zones and 2000–6000 nT caution zones [11].

Table 2 is an example budget for Ankara (39.93° N, 32.86° E, 2026-09-25; WMM2025: F = 48.19 μT, H = 25.40 μT, I = 58.19° [12]); inputs are illustrative except the data-sheet row.

Table 2. Example heading error budget, H = 25.40 μT, I = 58.19°

Root-sum-square of the “good calibration” items (0.1 μT, 1%, both 0.5° rotations, 0.2° tilt, averaged noise, WMM) is about 1.1° (derived); RSS assumes independence and mixes peaks with 1σ, so it indicates magnitude only. Alignment and board currents dominate, and neither shows in the fit residual. The same 0.5 μT gives 2.65° at 69.7° N, where H = 10.8 μT (WMM2025 [12]).

Practical problems

Temperature drift

On Ørsted's fluxgate, sensitivities and offsets varied linearly with temperature, non-orthogonalities did not; the in-flight model has linear terms in sensor and electronics temperature and in time [13]. A consumer MEMS data sheet gives a sensitivity temperature coefficient of ±0.01 %/K typical, ±0.03 %/K maximum, and a typical zero-field offset of ±2 μT over −40…+85 °C after calibration with the vendor's software [15]. Common gain change cancels in heading, since only differences enter (15); offset drift enters (14) directly. The 4.5° in Table 2 is that whole typical residual; the data sheet does not separate the thermal share. Remedies: a temperature-dependent offset model from climatic-chamber characterisation, in-use recalibration, or both.

Board currents and dynamic disturbers

For a long straight trace, Ampère's law gives

so 1 A gives 200 μT at 1 mm and 20 μT at 1 cm [2]. Constant currents join the hard-iron offset; load-dependent ones are the problem, and a supply trace following processor load, backlight or RF amplifier near the sensor is called the hardest case [2]. A 50 mA change at 10 mm is 1.0 μT, i.e. 2.3° at H = 25.4 μT (derived). With sampling up to about 100 Hz, fast currents are unseen while slow ones modulate heading [2]. A steel shield creates strong soft iron and screens the geomagnetic field; any shield should be low-permeability and at least 7 mm away [2].

Hard and soft iron that change after assembly

Calibration holds only in the final mechanical configuration. In Renaudin et al., a taped screw moved the combined offset from (−4.64, −0.24, 2.95) μT to (4.29, 44.15, −27.86) μT and changed off-diagonal terms; calibration cut the maximum heading difference from 239° to about 5° with the screw and from 18° to ±2° otherwise [7]. Ørsted's pre-flight/in-flight offset difference was attributed to tests not run in flight configuration (flight battery, live solar panels) [13]. A hinged part open or closed is another configuration [2]; strong-field exposure can change magnetisation, and hysteresis breaks linearity [1].

Repeatability in production

Either calibrate every unit fully, or use a per-design soft-iron matrix with per-unit offsets; the latter assumes little unit-to-unit soft-iron variation, which must be verified statistically on a lot. Station cleanliness (fixtures, motors, steel benches) is part of the calibration; a constant indoor field takes care [9].

Validation infrastructure and method

The ellipsoid residual cannot see the orthogonal ambiguity or out-of-coverage error; validation needs an independent reference.

• Helmholtz coils: nulling the geomagnetic field measures the zero-field offset; known vectors give scale and cross-axis gain [7]. Verify coil constants with a reference fluxgate; monitor coil current and temperature.
• Scalar reference: Merayo et al. estimated three offsets, scales and non-orthogonality angles using only a scalar reference in the Earth's field, reporting 0.5 nT absolute precision [14]. On Ørsted the calibrated vector instrument agreed with the scalar Overhauser magnetometer to 0.33 nT rms; the rotation came separately from a star imager, to better than 4 arcsec [13].
• Non-magnetic turntable: three-axis platforms of aluminium, brass, plastic and glass with piezo motors have been used [7]; steel has relative permeability ≈ 100, aluminium and brass ≈ 1 [2].
• Heading reference: GNSS/INS [7], navigation-grade INS [10] or 90° steps of a right-angled block [9]; the block misses the constant A in (16) but catches the harmonics.
• Climatic chamber: fans, heaters and metal structure can distort the test-volume field; map it and wait for thermal equilibrium.

Metrics: heading error against the reference decomposed into the harmonics of (16); tests at tilts outside the calibration envelope; residuals on held-out data. Kok and Schön check that normalised EKF residuals resemble N(0, 1) on estimation and validation sets [9]; in Gebre-Egziabher et al. the estimator covariance bounded actual errors, a conservative metric [10].

Common mistakes

• Min/max offset and scale: axis-extreme midpoints ignore cross-axis coupling [7] and rely on the noisiest samples. On the Table 1 full-sphere data, per-axis min/max (offset plus diagonal scale) found the centre within about 0.5 μT but left 5.7° maximum heading error from off-diagonal soft iron; the 10-parameter fit gave 0.2° (median of 100 trials, derived).
• Assuming a sphere: fine for weak soft iron [1]; otherwise the offset is right but a second harmonic remains — 0.3 μT centre error but 7.3° heading error in the same simulation.
• Judging the fit by one number: residual RMS or normalised fit error [1] is blind to the orthogonal ambiguity and flatters narrow-coverage fits (Table 1, ±20°).
• Treating the symmetric correction as physical truth: by (12) misalignment remains; 1° horizontal is about 1.6° of heading error.
• Calibrating in the lab, using in the field: a different assembly means different parameters [7], [13].
• Ignoring saturation: exceeding the measurement range, offsets included, invalidates the linear model [15].

References

1. T. Ozyagcilar, “Calibrating an eCompass in the Presence of Hard- and Soft-Iron Interference,” Freescale/NXP Application Note AN4246, Rev. 4.0, 2015. https://www.nxp.com/docs/en/application-note/AN4246.pdf
2. T. Ozyagcilar, “Layout Recommendations for PCBs Using a Magnetometer Sensor,” Freescale/NXP Application Note AN4247, Rev. 4.0, 2015. https://www.nxp.com/docs/en/application-note/AN4247.pdf
3. T. Ozyagcilar, “Accuracy of Angle Estimation in eCompass and 3D Pointer Applications,” Freescale/NXP Application Note AN4249, Rev. 1.0, 2015. https://www.nxp.com/docs/en/application-note/AN4249.pdf
4. Q. Li, J. G. Griffiths, “Least squares ellipsoid specific fitting,” Geometric Modeling and Processing 2004, Proceedings, IEEE, pp. 335–340, 2004. https://doi.org/10.1109/GMAP.2004.1290055
5. M. Zhao, X. Jia, L. Ma, Y. Shi, J. Jiang, Q. Li, D.-M. Yan, T. Huang, “A Bayesian Approach Toward Robust Multidimensional Ellipsoid-Specific Fitting,” arXiv:2407.19269, 2024. https://arxiv.org/abs/2407.19269
6. K. Usevich, I. Markovsky, “Adjusted least squares fitting of algebraic hypersurfaces,” arXiv:1412.2291, 2015. https://arxiv.org/abs/1412.2291
7. V. Renaudin, M. H. Afzal, G. Lachapelle, “Complete Triaxis Magnetometer Calibration in the Magnetic Domain,” Journal of Sensors, vol. 2010, art. 967245, 2010. https://doi.org/10.1155/2010/967245
8. J. F. Vasconcelos, G. Elkaim, C. Silvestre, P. Oliveira, B. Cardeira, “Geometric Approach to Strapdown Magnetometer Calibration in Sensor Frame,” IEEE Transactions on Aerospace and Electronic Systems, vol. 47, no. 2, pp. 1293–1306, 2011. https://doi.org/10.1109/TAES.2011.5751259
9. M. Kok, T. B. Schön, “Magnetometer Calibration Using Inertial Sensors,” IEEE Sensors Journal, vol. 16, no. 14, pp. 5679–5689, 2016. https://doi.org/10.1109/JSEN.2016.2569160
10. D. Gebre-Egziabher, G. H. Elkaim, J. D. Powell, B. W. Parkinson, “Calibration of Strapdown Magnetometers in Magnetic Field Domain,” Journal of Aerospace Engineering, vol. 19, no. 2, pp. 87–102, 2006. https://doi.org/10.1061/(ASCE)0893-1321(2006)19:2(87)
11. NOAA NCEI, “World Magnetic Model: Accuracy, Limitations, and Error Model (WMM2025),” accessed 2026-09-25. https://www.ncei.noaa.gov/products/world-magnetic-model/accuracy-limitations-error-model
12. A. Chulliat et al., “The US/UK World Magnetic Model for 2025–2030: Technical Report,” NOAA National Centers for Environmental Information, 2025. https://doi.org/10.25923/prbc-s316
13. N. Olsen, L. Tøffner-Clausen, T. J. Sabaka, P. Brauer, J. M. G. Merayo, J. L. Jørgensen, J.-M. Léger, O. V. Nielsen, F. Primdahl, T. Risbo, “Calibration of the Ørsted vector magnetometer,” Earth, Planets and Space, vol. 55, no. 1, pp. 11–18, 2003. https://doi.org/10.1186/BF03352458
14. J. M. G. Merayo, P. Brauer, F. Primdahl, J. R. Petersen, O. V. Nielsen, “Scalar calibration of vector magnetometers,” Measurement Science and Technology, vol. 11, no. 2, pp. 120–132, 2000. https://doi.org/10.1088/0957-0233/11/2/304
15. Bosch Sensortec, “BMM150 Geomagnetic Sensor — Data sheet,” BST-BMM150-DS001-05, Rev. 1.4, 2020. https://www.bosch-sensortec.com/media/boschsensortec/downloads/datasheets/bst-bmm150-ds001.pdf