Calibration of 3D field sensors

0Associated
Cases 
0Associated
Defendants 
0Accused
Products 
42Forward
Citations 
0
Petitions 
2
Assignments
First Claim
1. An apparatusbased method for calibration of a threedimensional (3D) field sensor, comprising:
 accessing a plurality of samples, the samples from the 3D field sensor, each sample representing a magnitude and an orientation of a 3D field sensed by the 3D field sensor;
using the plurality of samples, determining a plurality of parameters of an ellipsoid, wherein the determining the plurality of parameters is performed so that the ellipsoid fits the plurality of samples, wherein the following is used to represent the ellipsoid;
x^{2}+y^{2}+z^{2}−U(x^{2}+y^{2}−2z^{2})−V(x^{2}−2y^{2}+z^{2})−4Mxy−2Nxz−2Pyz−Qx−Ry−Sz−T=0, and wherein determining a plurality of parameters determines the parameters {right arrow over (s)}_{LS}=(U,V,M,N,P,Q,R,S,T)^{T};
determining a transformation that transforms the ellipsoid into a sphere;
applying the transformation to a sample to create a transformed sample; and
using the transformed sample to determine at least an orientation of the apparatus.
2 Assignments
0 Petitions
Accused Products
Abstract
An apparatusbased method is disclosed for calibration of a 3D field sensor. The method includes accessing a plurality of samples, where the samples are from the 3D field sensor. Each sample represents a magnitude and an orientation of a threedimensional field sensed by the 3D field sensor. Using the plurality of samples, a plurality of parameters is determined of an ellipsoid. The determination of the plurality of parameters is performed so that the ellipsoid fits the plurality of samples. A transformation is determined that transforms the ellipsoid into a sphere. The transformation is applied to a sample to create a transformed sample. Apparatus and signal bearing media are also disclosed.
49 Citations
View as Search Results
Calibrating sensor unit orientation for use in a vehicle monitoring system  
Patent #
US 9,958,473 B1
Filed 05/01/2017

Current Assignee
SmartDrive Systems Inc.

Sponsoring Entity
SmartDrive Systems Inc.

GYROSCOPE CALIBRATION  
Patent #
US 20150260545A1
Filed 05/27/2015

Current Assignee
Google LLC

Sponsoring Entity
Google LLC

Interfacing application programs and motion sensors of a device  
Patent #
US 9,811,174 B2
Filed 04/28/2015

Current Assignee
Invensense Incorporated

Sponsoring Entity
Invensense Incorporated

Systems and methods for efficiently generating magnetic resonance images from incomplete data  
Patent #
US 10,048,338 B2
Filed 04/23/2015

Current Assignee
Mayo Foundation For Medical Education And Research

Sponsoring Entity
Mayo Foundation For Medical Education And Research

Systems and Methods for Efficiently Generating Magnetic Resonance Images from Incomplete Data  
Patent #
US 20150309136A1
Filed 04/23/2015

Current Assignee
Mayo Foundation For Medical Education And Research

Sponsoring Entity
Mayo Foundation For Medical Education And Research

Integrated motion processing unit (MPU) with MEMS inertial sensing and embedded digital electronics  
Patent #
US 10,288,427 B2
Filed 03/19/2015

Current Assignee
Invensense Incorporated

Sponsoring Entity
Invensense Incorporated

MEMS rotation sensor with integrated electronics  
Patent #
US 9,846,175 B2
Filed 01/27/2015

Current Assignee
Invensense Incorporated

Sponsoring Entity
Invensense Incorporated

Interfacing application programs and motion sensors of a device  
Patent #
US 9,342,154 B2
Filed 01/06/2015

Current Assignee
Invensense Incorporated

Sponsoring Entity
Invensense Incorporated

DYNAMIC MAGNETOMETER CALIBRATION  
Patent #
US 20150018099A1
Filed 09/29/2014

Current Assignee
Microsoft Technology Licensing LLC

Sponsoring Entity
Microsoft Technology Licensing LLC

Dynamic magnetometer calibration  
Patent #
US 9,250,300 B2
Filed 09/29/2014

Current Assignee
Microsoft Technology Licensing LLC

Sponsoring Entity
Microsoft Technology Licensing LLC

PORTABLE DEVICE AND METHOD FOR CONTROLLING THE SAME  
Patent #
US 20150070395A1
Filed 09/12/2014

Current Assignee
KIM Seungil

Sponsoring Entity
KIM Seungil

Physical quantity data correcting device and physical quantity data correcting method  
Patent #
US 9,897,446 B2
Filed 03/05/2014

Current Assignee
Asahi Kasei Microdevices Corporation

Sponsoring Entity
Asahi Kasei Microdevices Corporation

Electromagnetic position and orientation sensing system  
Patent #
US 8,723,509 B2
Filed 05/21/2013

Current Assignee
Brown University

Sponsoring Entity
Brown University

Magnetization device and method  
Patent #
US 9,257,220 B2
Filed 03/05/2013

Current Assignee
eZono AG

Sponsoring Entity
eZono AG

Magnetic position detection system  
Patent #
US 9,459,087 B2
Filed 03/05/2013

Current Assignee
eZono AG

Sponsoring Entity
eZono AG

Distributed systems and methods to measure and process sport motions  
Patent #
US 9,060,682 B2
Filed 10/25/2012

Current Assignee
AlpineReplay Inc.

Sponsoring Entity
AlpineReplay Inc.

Inuse automatic calibration methodology for sensors in mobile devices  
Patent #
US 9,683,865 B2
Filed 10/12/2012

Current Assignee
Invensense Incorporated

Sponsoring Entity
Invensense Incorporated

Magnetometer bias and anomaly detector  
Patent #
US 8,577,640 B2
Filed 08/10/2012

Current Assignee
Invensense Incorporated

Sponsoring Entity
Invensense Incorporated

Integrated motion processing unit (MPU) with MEMS inertial sensing and embedded digital electronics  
Patent #
US 8,997,564 B2
Filed 06/08/2012

Current Assignee
Invensense Incorporated

Sponsoring Entity
Invensense Incorporated

Wireless controller  
Patent #
US 9,849,376 B2
Filed 05/02/2012

Current Assignee
Microsoft Technology Licensing LLC

Sponsoring Entity
Microsoft Technology Licensing LLC

Low inertia frame for detecting coriolis acceleration  
Patent #
US 8,539,835 B2
Filed 03/22/2012

Current Assignee
Invensense Incorporated

Sponsoring Entity
Invensense Incorporated

Device and method of gyro sensor calibration  
Patent #
US 9,146,134 B2
Filed 11/08/2011

Current Assignee
Qualcomm Inc.

Sponsoring Entity
AlpineReplay Inc.

Imaging probe and method of obtaining position and/or orientation information  
Patent #
US 9,597,008 B2
Filed 09/06/2011

Current Assignee
eZono AG

Sponsoring Entity
eZono AG

Reducing power consumption or error of digital compass  
Patent #
US 9,541,393 B2
Filed 06/30/2011

Current Assignee
Qualcomm Inc.

Sponsoring Entity
Qualcomm Inc.

Vertically integrated 3axis MEMS angular accelerometer with integrated electronics  
Patent #
US 8,960,002 B2
Filed 04/28/2011

Current Assignee
Invensense Incorporated

Sponsoring Entity
Invensense Incorporated

MEASUREMENT METHOD FOR A COMPONENT OF THE GRAVITY VECTOR  
Patent #
US 20110196636A1
Filed 01/28/2011

Current Assignee
Baker Hughes Incorporated

Sponsoring Entity
Baker Hughes Incorporated

CALIBRATING MULTIDIMENSIONAL SENSOR FOR OFFSET, SENSITIVITY, AND NONORTHOGONALITY  
Patent #
US 20110106477A1
Filed 11/04/2009

Current Assignee
Qualcomm Inc.

Sponsoring Entity
Qualcomm Inc.

Calibrating multidimensional sensor for offset, sensitivity, and nonorthogonality  
Patent #
US 8,645,093 B2
Filed 11/04/2009

Current Assignee
Qualcomm Inc.

Sponsoring Entity
Qualcomm Inc.

Accurate magnetic compass in mobile electronic device  
Patent #
US 8,880,373 B2
Filed 11/04/2009

Current Assignee
Qualcomm Inc.

Sponsoring Entity
Qualcomm Inc.

SYSTEMS AND METHODS FOR CALIBRATION OF GYROSCOPES AND A MAGNETIC COMPASS  
Patent #
US 20110066392A1
Filed 09/14/2009

Current Assignee
Honeywell International Inc.

Sponsoring Entity
Honeywell International Inc.

Systems and methods for calibration of gyroscopes and a magnetic compass  
Patent #
US 8,108,171 B2
Filed 09/14/2009

Current Assignee
Honeywell International Inc.

Sponsoring Entity
Honeywell International Inc.

Systems and methods for gyroscope calibration  
Patent #
US 8,548,766 B2
Filed 09/14/2009

Current Assignee
Honeywell International Inc.

Sponsoring Entity
Honeywell International Inc.

STORAGE MEDIUM STORING INFORMATION PROCESSING PROGRAM, AND INFORMATION PROCESSING APPARATUS  
Patent #
US 20100286940A1
Filed 07/28/2009

Current Assignee
Nintendo Company Limited

Sponsoring Entity
Nintendo Company Limited

Storage medium storing information processing program, and information processing apparatus  
Patent #
US 9,050,525 B2
Filed 07/28/2009

Current Assignee
Nintendo Company Limited

Sponsoring Entity
Nintendo Company Limited

Controlling and accessing content using motion processing on mobile devices  
Patent #
US 8,462,109 B2
Filed 06/16/2009

Current Assignee
Invensense Incorporated

Sponsoring Entity
Invensense Incorporated

SYSTEM FOR MEASURING A MAGNETIC FILED COMPRISING A THREEAXIS SENSOR FOR MEASURING A MAGNETIC FILED THAT IS ABLE TO MOVE TOGETHER WITH A CARRIER THAT DISRUPTS THE MEASUREMENTS, AND ASSOCIATED METHOD  
Patent #
US 20110153233A1
Filed 05/14/2009

Current Assignee
Movea, Commissariat a LEnergie Atomique

Sponsoring Entity
Movea, Commissariat a LEnergie Atomique

System for measuring a magnetic field comprising a threeaxis sensor for measuring a magnetic field that is able to move together with a carrier that disrupts the measurements, and associated method  
Patent #
US 9,459,098 B2
Filed 05/14/2009

Current Assignee
Movea, Commissariat a LEnergie Atomique

Sponsoring Entity
Movea, Commissariat a LEnergie Atomique

DERIVING TILTCORRECTED SEISMIC DATA IN A MULTIAXIS SEISMIC SENSOR MODULE  
Patent #
US 20100195438A1
Filed 02/05/2009

Current Assignee
WesternGeco LLC

Sponsoring Entity
WesternGeco LLC

Deriving tiltcorrected seismic data in a multiaxis seismic sensor module  
Patent #
US 8,199,611 B2
Filed 02/05/2009

Current Assignee
WesternGeco LLC

Sponsoring Entity
WesternGeco LLC

MOBILE DEVICES WITH MOTION GESTURE RECOGNITION  
Patent #
US 20090265671A1
Filed 10/15/2008

Current Assignee
Invensense Incorporated

Sponsoring Entity
Invensense Incorporated

Wafer scale chip scale packaging of vertically integrated MEMS sensors with electronics  
Patent #
US 8,508,039 B1
Filed 05/08/2008

Current Assignee
Invensense Incorporated

Sponsoring Entity
Invensense Incorporated

Interfacing application programs and motion sensors of a device  
Patent #
US 8,952,832 B2
Filed 04/21/2008

Current Assignee
Invensense Incorporated

Sponsoring Entity
Invensense Incorporated

Downhole calibration system for directional sensors  
Patent #
US 20040149004A1
Filed 02/04/2003

Current Assignee
Weatherford Technology Holdings LLC

Sponsoring Entity


Device for calibration of magnetic sensors in three dimensions  
Patent #
US 20060164073A1
Filed 07/01/2002

Current Assignee
EUROPEAN ORGANISATION FOR NUCLEAR RESEARCH

Sponsoring Entity
EUROPEAN ORGANISATION FOR NUCLEAR RESEARCH

Three dimensional magnetic field sensor  
Patent #
US 6,278,271 B1
Filed 03/29/1999

Current Assignee
MELEXIS TESSENDERLO NV

Sponsoring Entity
SENTRON AG

Navigation system and method  
Patent #
US 6,014,610 A
Filed 01/30/1998

Current Assignee
Continental Automotive Systems US Incorporated

Sponsoring Entity
GREENFIELD ENTERPRISES INC. DBA GREENFIELD ASS.

Threeaxis fiber optic vector magnetometer  
Patent #
US 5,243,403 A
Filed 09/30/1991

Current Assignee
UNITED STATES OF AMERICA THE AS REPRESENTED BY THE SECRETARY OF THE NAVY

Sponsoring Entity
UNITED STATES OF AMERICA THE AS REPRESENTED BY THE SECRETARY OF THE NAVY

Calibration of magnetic compasses  
Patent #
US 5,117,375 A
Filed 10/19/1989

Current Assignee
BAE Systems plc

Sponsoring Entity
British Aerospace Inc.

Method for determining route angles  
Patent #
US 4,870,602 A
Filed 12/19/1986

Current Assignee
LITEF GMBH

Sponsoring Entity
LITEF GMBH

41 Claims
 1. An apparatusbased method for calibration of a threedimensional (3D) field sensor, comprising:
 accessing a plurality of samples, the samples from the 3D field sensor, each sample representing a magnitude and an orientation of a 3D field sensed by the 3D field sensor;
using the plurality of samples, determining a plurality of parameters of an ellipsoid, wherein the determining the plurality of parameters is performed so that the ellipsoid fits the plurality of samples, wherein the following is used to represent the ellipsoid;
x^{2}+y^{2}+z^{2}−U(x^{2}+y^{2}−2z^{2})−V(x^{2}−2y^{2}+z^{2})−4Mxy−2Nxz−2Pyz−Qx−Ry−Sz−T=0, and wherein determining a plurality of parameters determines the parameters {right arrow over (s)}_{LS}=(U,V,M,N,P,Q,R,S,T)^{T};
determining a transformation that transforms the ellipsoid into a sphere;
applying the transformation to a sample to create a transformed sample; and
using the transformed sample to determine at least an orientation of the apparatus.  View Dependent Claims (2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12)
 accessing a plurality of samples, the samples from the 3D field sensor, each sample representing a magnitude and an orientation of a 3D field sensed by the 3D field sensor;
 13. (canceled)
 22. A computerreadable medium tangibly embodying a program of machinereadable instructions executable by a processor to perform operations for calibration of a threedimensional (3D) field sensor, the operations comprising:
 accessing a plurality of samples, the samples from the 3D field sensor, each sample representing a magnitude and an orientation of a threedimensional field sensed by the 3D field sensor;
using the plurality of samples, determining a plurality of parameters of an ellipsoid, wherein the determining of the plurality of parameters is performed so that the ellipsoid fits the plurality of samples, wherein the following is used to represent the ellipsoid;
x^{2}y^{2}+z^{2}−U(x^{2}+y^{2}−2z^{2})−V(x^{2}−2y^{2}+z^{2})−4Mxy−2Nxz−2Pyz−Qx−Ry−Sz−T=0, and wherein determining a plurality of parameters determines the parameters {right arrow over (s)}_{LS}=(U,V,M,N,P,Q,R,S,T)^{T};
determining a transformation that transforms the ellipsoid into a sphere;
applying the transformation to a sample to create a transformed sample; and
using the transformed sample to determine at least an orientation of the apparatus.  View Dependent Claims (23)
 accessing a plurality of samples, the samples from the 3D field sensor, each sample representing a magnitude and an orientation of a threedimensional field sensed by the 3D field sensor;
 24. An apparatus, comprising:
 means for measuring a threedimensional (3D) field and for producing samples representing a magnitude and an orientation of the threedimensional field;
means for accessing a plurality of samples from the means for measuring;
means, responsive to the plurality of samples, for determining a plurality of parameters of an ellipsoid, wherein the determining of the plurality of parameters is performed so that the ellipsoid fits the plurality of samples, wherein the following is used to represent the ellipsoid;
x^{2}+y^{2}+z^{2}−U(x^{2}+y^{2}−2z^{2})−V(x^{2}−2y^{2}+z^{2})−4Mxy−2Nxz−2Pyz−Qx−Ry−Sz−T=0, and wherein determining a plurality of parameters determines the parameters {right arrow over (s)}_{LS}=(U,V,M,N,P,Q,R,S,T)^{T};
means for determining a transformation that transforms the ellipsoid into a sphere;
means for applying the transformation to a sample to create a transformed sample; and
means for using the transformed sample to determine at least an orientation of the apparatus.  View Dependent Claims (25)
 means for measuring a threedimensional (3D) field and for producing samples representing a magnitude and an orientation of the threedimensional field;
 26. An apparatus, comprising:
 a threedimensional (3D) field sensor, the 3D field sensor adapted to produce field samples, each field sample representing a magnitude and an orientation of a 3D field sensed by the 3D field sensor;
a calibration module adapted to access a plurality of the field samples from the 3D field sensor, the calibration module adapted to use the plurality of samples to determine a plurality of parameters of an ellipsoid, wherein the determining the plurality of parameters operation is performed so that the ellipsoid fits the plurality of field samples, wherein the following is used to represent the ellipsoid;
x^{2}+y^{2}+z^{2}−U(x^{2}+y^{2}−2z^{2})−V(x^{2}−2y^{2}+z^{2})−4Mxy−2Nxz−2Pyz−Qx−Ry−Sz−T=0, wherein the calibration module is adapted to determines the parameters {right arrow over (s)}_{LS}=(U,V,M,N,P,Q,R,S,T)^{T}, and wherein the calibration module is further adapted to determine a transformation that transforms the ellipsoid into a sphere and to apply the transformation to a field sample to create a transformed field sample; and
an application configured to use the transformed sample to determine at least an orientation of the apparatus.  View Dependent Claims (27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37)
 a threedimensional (3D) field sensor, the 3D field sensor adapted to produce field samples, each field sample representing a magnitude and an orientation of a 3D field sensed by the 3D field sensor;
 38. A method comprising:
 accessing a plurality of samples, the samples from the 3D field sensor, each sample representing a magnitude and an orientation of a 3D field sensed by the 3D field sensor;
using the plurality of samples, determining a plurality of parameters of an ellipsoid, wherein the determining the plurality of parameters is performed so that the ellipsoid fits the plurality of samples;
determining a transformation that transforms the ellipsoid into a sphere, wherein determining a transformation comprises;
determining a matrix, A, a bias, {right arrow over (b)}, and a scaling factor, c in a representation of an ellipsoid of ({right arrow over (x)}−{right arrow over (b)})^{T}A({right arrow over (x)}−{right arrow over (b)})=c;
determining a matrix, D, where D^{T}D=1/cA, where D is selected such that Dˆ=λˆ, where λ is an arbitrary constant and ˆ=(1,0,0)^{T}, and where the transformation comprises the matrix D;
applying the transformation to a sample to create a transformed sample; and
. using the transformed sample to determine at least an orientation of the apparatus.  View Dependent Claims (39)
 accessing a plurality of samples, the samples from the 3D field sensor, each sample representing a magnitude and an orientation of a 3D field sensed by the 3D field sensor;
 40. An apparatus, comprising:
 a threedimensional (3D) field sensor, the 3D field sensor configured to produce field samples, each field sample representing a magnitude and an orientation of a 3D field sensed by the 3D field sensor; and
a calibration module configured to access a plurality of the field samples from the 3D field sensor, the calibration module adapted to use the plurality of samples to determine a plurality of parameters of an ellipsoid, wherein the determining the plurality of parameters operation is performed so that the ellipsoid fits the plurality of field samples, wherein the calibration module is further configured determine a transformation that transforms the ellipsoid into a sphere by determining a matrix, A, a bias, {right arrow over (b)}, and a scaling factor, c in a representation of an ellipsoid of ({right arrow over (x)}−{right arrow over (b)})^{T}A({right arrow over (x)}−{right arrow over (b)})=c, and by determining a matrix, D, where D^{T}D=1/cA, where D is selected such that Dˆ=λˆ, where λ is an arbitrary constant and ˆ=(1,0,0)^{T}, and where the transformation comprises the matrix D, and wherein the calibration module is further configured to apply the transformation to a field sample to create a transformed field sample; and
an application configured to use the transformed sample to determine at least an orientation of the apparatus.  View Dependent Claims (41)
 a threedimensional (3D) field sensor, the 3D field sensor configured to produce field samples, each field sample representing a magnitude and an orientation of a 3D field sensed by the 3D field sensor; and
1 Specification
This application is a continuation in part of U.S. patent application Ser. No. 11/219,576, filed on Sep. 2, 2005, which is herein incorporated by reference.
TECHNICAL FIELDThis invention relates generally to sensors that measure threedimensional (3D) fields and, more specifically, relates to calibration of such sensors.
BACKGROUND OF THE INVENTIONA 3D field sensor is a sensor that measures direction and magnitude of some field in threedimensional (3D) space. For example three orthogonal magnetometer elements can measure the direction and magnitude of the surrounding magnetic field. Similarly, three orthogonal accelerometer elements can sense the direction and magnitude of acceleration or more importantly the gravitational field, which is typically indistinguishable from acceleration.
As already suggested above, in almost all cases, a 3D field sensor includes at some level three singledimensional sensors oriented appropriately relative to each other. The basic assumption in reconstructing the direction and magnitude of the measured field is that these elements are orthogonal to each other and have similar scales and bias of zero.
In practice, however, these assumptions are rarely valid. Magnetometers are especially tricky, because external magnetic components attached to a device in which the magnetometer resides can result in scale and offset errors and can also effectively introduce linear dependence among the measuring axes.
To counteract these problems, calibration of the 3D field sensor is needed. Usually this process is lengthy and generally involves complicated steps for the user. For example, to find the scale of an axis, it is usually required to find the maximum and minimum value that the field can produce on that axis. This process also gives the offset, but cannot find possible linear dependence among different axes. To find possible linear dependence among different axes, it is usually required to indicate certain reference directions, which is also time consuming for the user.
As another example, a traditional way of calibrating an electronic compass is by comparing measured headings to known headings and building a lookup table. To have an accurate heading, this process must be performed whenever the device is moved to a new location. Unfortunately, 3D magnetometers are relatively new, and most calibration methods are for 2D sensors. Automatic calibration methods have been based on Kalman filtering, in for example: Bruce Hoff, Ronald Azuma: “Autocalibration of an Electronic Compass in an Outdoor Augmented Reality System”, Proceedings of IEEE and ACM Int'l Symposium on Augmented Reality 2000 (2000). This process still takes user input regarding known headings and is typically impractical for 3D field sensors.
It would therefore be desirable to provide techniques that can calibrate 3D field sensors with little if any user input.
BRIEF SUMMARY OF THE INVENTIONThe present invention provides techniques that calibrate 3D field sensors.
In an exemplary embodiment of the invention, an apparatusbased method is disclosed for calibration of a 3D field sensor. The method comprises accessing a plurality of samples, where the samples are from the 3D field sensor. Each sample represents a magnitude and an orientation of a threedimensional field sensed by the 3D field sensor. Using the plurality of samples, a plurality of parameters is determined of an ellipsoid. The determination of the parameters is performed so that the ellipsoid fits the plurality of samples. A transformation is determined that transforms the ellipsoid into a sphere. The transformation is applied to a sample to create a transformed sample.
In another exemplary embodiment, a signal bearing medium is disclosed that tangibly embodies a program of machinereadable instructions executable by a processor to perform operations for calibration of a threedimensional (3D) field sensor. The operations include accessing a plurality of samples, where the samples from the 3D field sensor. Each sample represents a magnitude and an orientation of a threedimensional field sensed by the 3D field sensor. Using the plurality of samples, a determination is made of a plurality of parameters of an ellipsoid. The determination of the plurality of parameters is performed so that the ellipsoid fits the plurality of samples. A transformation is determined that transforms the ellipsoid into a sphere. The transformation is applied to a sample to create a transformed sample.
In another exemplary embodiment, an apparatus is disclosed that includes a means for measuring a threedimensional (3D) field and for producing samples representing a magnitude and an orientation of the threedimensional field. The apparatus also includes means for accessing a plurality of samples from the means for measuring. The apparatus further includes means, responsive to the plurality of samples, for determining a plurality of parameters of an ellipsoid, wherein the determining of the plurality of parameters is performed so that the ellipsoid fits the plurality of samples. The apparatus additionally has means for determining a transformation that transforms the ellipsoid into a sphere and has means for applying the transformation to a sample to create a transformed sample.
In a further exemplary embodiment, another apparatus is disclosed that includes a threedimensional (3D) field sensor. The 3D field sensor is adapted to produce field samples, where each field sample represents a magnitude and an orientation of a 3D field sensed by the 3D field sensor. The apparatus includes a calibration module adapted to access a plurality of the field samples from the 3D field sensor. The calibration module is adapted to use the plurality of samples to determine a plurality of parameters of an ellipsoid, wherein the determining the plurality of parameters operation is performed so that the ellipsoid fits the plurality of field samples. The calibration module is further adapted to determine a transformation that transforms the ellipsoid into a sphere and to apply the transformation to a field sample to create a transformed field sample.
BRIEF DESCRIPTION OF THE DRAWINGSThe foregoing and other aspects of embodiments of this invention are made more evident in the following Detailed Description of Exemplary Embodiments, when read in conjunction with the attached Drawing Figures, wherein:
FIG. 1 is a block diagram of an exemplary apparatus for calibrating a 3D field sensor;
FIG. 2 is an exemplary representation of a transformation from an ellipsoid to a sphere;
FIG. 3 is a flowchart of an exemplary method for calibrating a 3D field sensor;
FIG. 4 is a block diagram of an exemplary apparatus for calibrating a 3D field sensor;
FIG. 5 is a flowchart of an exemplary method for filtering samples from a 3D field sensor; and
FIG. 6 is a block diagram of an apparatus including a 3D field sensor and modules for calibration of the 3D field sensor.
DETAILED DESCRIPTION OF EXEMPLARY EMBODIMENTSFor ease of reference, the present disclosure is divided into a number of different sections.
1. Introduction
The present disclosure is primarily directed to calibrations of 3D field sensors such as threeaxis magnetometers (TAM), but in principle the techniques herein can be applied to any threeaxis sensor that is used to sense a linear field. Examples of such sensors could be a threeaxis accelerometer sensing the earth's gravity or a sensor sensing an electric field in three dimensions.
Of course, there is any number of straightforward ways of calibrating such a sensor with some kind of reference signal(s). Our setting, however, is different. A TAM can be used in conjunction with a threeaxis accelerometer to obtain the orientation of the device in reference to the coordinate system set by the earth's gravitational and magnetic fields. In astronautical sciences, this is usually called the attitude of the spacecraft. We, however, are interested in the orientation of apparatus including hand held devices, such as mobile phones or game controllers. Such devices typically use inexpensive magnetometers and as a result the calibration is likely to drift in time. We therefore need an adaptive calibration that is invisible to the user.
A similar problem is faced with space crafts that are interested in measuring the earth's magnetic field. The calibration is likely to change due to the harsh conditions of take off and the different atmosphere in orbit. Fully known reference signals are not available in orbit either. There is a wealth of literature on the calibration of TAM for satellite systems.
In tests, we have observed that the calibration of a device is dependent on the location of the device. The calibration can be subject to significant changes when the device is moved as little as few meters in some direction. It is not surprising that the direction and magnitude of the magnetic field can change in such a short range, especially within buildings. Paramagnetic elements in buildings can strengthen the earth's magnetic field in one direction and electrical wiring may induce new magnetic fields. This should not, however, affect the calibration of the device sensing the field. One explanation might be that the paramagnetic components inside the device, near the magnetometers, affect the calibration differently depending on the strength of the observed field. Whatever the reason, this is observed behavior and one more reason for easy adaptive calibration.
Nevertheless, exemplary calibration herein is based on the assumption that the magnetic field is constant for sufficiently long time windows. This is in contrast with the problem for satellite systems, where the spacecraft travels in orbit and the magnetic field changes significantly in the course of one revolution around the earth. In the satellite systems case, some model of the earth's magnetic field is usually used. There is a significant difference between calibration for satellite systems and calibration for handheld devices. Obviously, when the device is rotated within the magnetic field, the sensed field is rotated in the coordinate system of the device. When the actual field is assumed to have a constant direction, then the relative angles of consecutive measurements contain information. If, however, the direction of the outside field may have changed as well, then this information may be lost. We will show that with an assumption of constant field, it is easy to determine a decoupling matrix for the calibration as well.
In the present disclosure, we discuss the calibration of a sensor that measures any 3D field, such as magnetic, gravitational or electric field. For instance, turning to FIG. 1, an exemplary apparatus 100 is shown for calibrating a 3D field sensor 110. The 3D field sensor 110 includes three singledimensional sensors 1201, 1202, and 1203 (s<sub>1</sub>, s<sub>2</sub>, and s<sub>3</sub>, respectively), each of which is parallel to an axis {right arrow over (c)}<sub>1</sub>, {right arrow over (c)}<sub>2</sub>, or {right arrow over (c)}<sub>3</sub>, which are orthogonal axes. Each sensor 120 measures an axial component of the field {right arrow over (F)}, which has a magnitude, m, and a direction, d. In the example of FIG. 1, the direction, d, is related by the angles θ, in the {right arrow over (c)}<sub>1</sub>, {right arrow over (c)}<sub>2 </sub>plane, and φ, an angle relative to the {right arrow over (c)}<sub>3 </sub>axis.
The 3D field sensor 110 produces 3D field samples 130, each of which could comprise spherical values (m, θ, φ) or field values (x,y,z) for each sensor 1201, 1202, and 1203, respectively. From field values (x,y,z), the magnitude and direction of the field can be determined. Note the 3D field sensor 110 is merely exemplary and used for exposition. The calibration module 140 operates on the 3D field samples 130 to create transformed 3D samples 150. The transformed 3D samples 150 should reduce or eliminate scale, offset, and other linear errors from the 3D field sensor 110.
When a sensor is measuring a constant 3D field (e.g., field {right arrow over (F)}), it can be assumed that all measurements (e.g., as defined by 3D field samples 130) from the 3D field sensor 110 should fall on the unit sphere. For instance, when a 3D field sensor 110 is rotated in a constant field, the measurements from the sensor 110 should lie on a spherical surface. This is shown graphically in FIG. 2 as unit sphere 220. The vectors {right arrow over (v)}<sub>1</sub>, {right arrow over (v)}<sub>2</sub>, {right arrow over (v)}<sub>3 </sub>should all be the same length in this example, and therefore measurements from 3D field sensor 110 should fall on the unit sphere 220. But misaligned sensor axes (e.g., axes {right arrow over (c)}<sub>1</sub>, {right arrow over (c)}<sub>2</sub>, and {right arrow over (c)}<sub>3</sub>), biases and scale errors will transform that sphere 220 into a quadratic surface, illustrated as exemplary ellipsoid 210, where one or more of the vectors {right arrow over (v)}′<sub>1</sub>,{right arrow over (v)}′<sub>2</sub>, {right arrow over (v)}′<sub>3 </sub>are not equivalent to their corresponding vectors {right arrow over (v)}<sub>1</sub>, {right arrow over (v)}<sub>2</sub>, {right arrow over (v)}<sub>3</sub>. Thus, one can fit a quadratic surface (e.g., ellipsoid 210) to the data (not shown in FIG. 2) from 3D field sensors 110 and find an inverse transformation that will turn that quadratic surface into a sphere. This process will then also make the measured data (e.g., the 3D field samples 130) spherical and thus this process has removed any offset, scale, and linear dependence of the 3D field sensor 110.
A model of the distortion herein is therefore is a quadratic surface such as ellipsoid 210, and the calibration module 140 finds the transformation 230 that brings that ellipsoid 210 back to the unit sphere 220. This process may be performed adaptively and computationally effectively—in real time, as new samples are collected. In particular, disclosed below is an exemplary method that can find a best fitting distortion (e.g., ellipsoid 210) recursively, so that relatively little computation is needed per new sample added to a set of samples used for calibration. This means that the calibration can be done online, during normal operation, hidden from the user, and the calibration can adapt to changes automatically without need for explicit recalibration.
The exemplary methods presented herein can determine the distortion (e.g., ellipsoid 210) from very little data. In particular, a reasonable calibration can be obtained without ever actually rotating a device a full turn. This is advantageous for example in mobile phones, where the orientation of the device is such that the screen is always facing the user and rarely away from the user. Also, the present invention has removed any explicit calibration step (unless such a step is desired), as the calibration can be updated from almost random data during operation.
In addition to offset and scale errors, 3D field sensor elements often suffer from nonlinearity. The disclosed methods herein are generally unable to correct for nonlinearities and may work incorrectly if measurements from a 3D field sensor 110 are highly nonlinear.
A short introduction will now be given to FIG. 3, which is a method 300 performed by the calibration module 140 for calibrating the 3D field sensor 110. More detailed description is given below. Briefly, the calibration module 140 accesses 3D field samples 130, which could be accessed for example from a memory (e.g., as shown in FIG. 6) or from the 3D field sensor 110. This occurs in step 310. The samples could be in any form suitable for describing a 3D field, but will typically be field values (x,y,z) for each sensor 1201, 1202, and 1203, respectively. In step 320, the calibration module 140 determines parameters for an ellipsoid (e.g., ellipsoid 210 of FIG. 2) that fits the 3D field samples 130. In step 330, a transformation is determined that transforms the ellipsoid into a sphere (e.g., unit sphere 220). In step 340, the transformation is applied to 3D field samples 130 to create transformed 3D samples 150. The transformed 3D samples 150 are therefore calibrated.
2. Distortion Model
It is assumed that a sensor measurement process can be modeled as the projection of the real field vector (e.g., field {right arrow over (F)}) on to the sensor axes and then adding an arbitrary offset. In formulas one can say that if {right arrow over (c)}<sub>1</sub>, {right arrow over (c)}<sub>2</sub>, and {right arrow over (c)}<sub>3 </sub>are the axes of the sensor, where their lengths represent a possible scale error, then the measurement process of the unit field {right arrow over (u)} is represented as<FORM>{right arrow over (m)}=C{right arrow over (u)}+{right arrow over (b)}+{right arrow over (ε)}, (1) </FORM>
where the rows of C are the {right arrow over (c)}<sub>i</sub>, {right arrow over (b)} is an arbitrary bias, and {right arrow over (ε)} is random measurement noise. Given {right arrow over (m)}, we are interested in knowing {right arrow over (u)}, so we need to know {right arrow over (b)} and D=C<sup>−1</sup>.<maths id="MATHUS00001" num="1"><math overflow="scroll"><mtable><mtr><mtd><mrow><mrow><mrow><mrow><mi>Since</mi><mo></mo><mstyle><mtext> </mtext></mstyle><mo></mo><mrow><mo></mo><mover><mi>u</mi><mo>></mo></mover><mo></mo></mrow></mrow><mo>=</mo><mn>1</mn></mrow><mo>,</mo><mrow><mi>we</mi><mo></mo><mstyle><mtext> </mtext></mstyle><mo></mo><mi>have</mi><mo></mo><mstyle><mtext> </mtext></mstyle><mo></mo><mi>from</mi><mo></mo><mstyle><mtext> </mtext></mstyle><mo></mo><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow><mo></mo><mstyle><mtext> </mtext></mstyle><mo></mo><mi>that</mi></mrow></mrow><mo></mo><mstyle><mtext></mtext></mstyle><mo></mo><mtable><mtr><mtd><mrow><msup><mrow><mo></mo><mrow><mi>D</mi><mo></mo><mrow><mo>(</mo><mrow><mover><mi>m</mi><mo>></mo></mover><mo></mo><mover><mi>b</mi><mo>></mo></mover></mrow><mo>)</mo></mrow></mrow><mo></mo></mrow><mn>2</mn></msup><mo>=</mo><msup><mrow><mo></mo><mrow><mover><mi>u</mi><mo>></mo></mover><mo></mo><mrow><mi>D</mi><mo></mo><mover><mi>ɛ</mi><mo>></mo></mover></mrow></mrow><mo></mo></mrow><mn>2</mn></msup></mrow></mtd></mtr><mtr><mtd><mrow><mo>=</mo><mrow><msup><mrow><mo></mo><mover><mi>u</mi><mo>></mo></mover><mo></mo></mrow><mn>2</mn></msup><mo>+</mo><mrow><mn>2</mn><mo></mo><msup><mover><mi>u</mi><mo>></mo></mover><mi>T</mi></msup><mo></mo><mi>D</mi><mo></mo><mover><mi>ɛ</mi><mo>></mo></mover></mrow><mo>+</mo><mrow><msup><mover><mi>ɛ</mi><mo>></mo></mover><mi>T</mi></msup><mo></mo><msup><mi>D</mi><mi>T</mi></msup><mo></mo><mi>D</mi><mo></mo><mrow><mover><mi>ɛ</mi><mo>></mo></mover><mo>.</mo></mrow></mrow></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mo>=</mo><mrow><mn>1</mn><mo>+</mo><mrow><msup><mi>ɛ</mi><mi>′</mi></msup><mo>.</mo></mrow></mrow></mrow></mtd></mtr></mtable></mrow></mtd><mtd><mrow><mo>(</mo><mn>2</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>
The points {right arrow over (m)} that satisfy (2) for ε′=0 form an ellipsoid, so it is beneficial to look for ellipsoid methods to solve for {right arrow over (m)}.
While ε′ is a random variable, ε′ has some unwanted properties. Namely, if the random variables {right arrow over (u)} and {right arrow over (ε)} are independent and zero mean, then E[ε′]=2E[{right arrow over (u)}<sup>T</sup>D{right arrow over (ε)}]+E└∥D{right arrow over (ε)}∥<sup>2</sup>┘≧0. So a least squares estimate may not be unbiased. Further ε′ depends on the very parameters we are trying to estimate. On the other hand, it is not uncommon in engineering to make assumptions about errors and in the limit as ε goes to zero this method provides the correction parameters exactly, so we ignore this problem and move on.
Notice that the resulting calibration parameters D and b only give a sort of relative calibration. With this correction, a rotation of α degrees about some axis results in a α degree turn in the corrected measurements. However, if D satisfies (2) then for any matrix R such that R<sup>T</sup>R=I, we have that RD also satisfies (2). So the actual orientation of the device relative to some reference coordinate system may actually be further lost as a result of the calibration. This problem, however, is addressed in section 5.
3. Fitting an Ellipsoid to Data
In a paper by D. A. Turner, I. J. Anderson, and J. C. Mason, entitled An algorithm for fitting an ellipsoid to data, Technical Report RR9803, School of Computing and Mathematics, University of Huddersfield (1999) (called “Turner” herein), the authors describe a method of fitting an ellipsoid to data using linear regression. The fitted ellipsoid in Turner is not the best fitting ellipsoid in terms of geometrical distance of the points to the surface. In fact, the authors in Turner note that a major drawback of the algorithm is the difficulty in interpreting the residuals that are minimized.
The authors in Turner give justification for the following quadratic form representing an ellipsoid<FORM>x<sup>2</sup>+y<sup>2</sup>+z<sup>2</sup>−U(x<sup>2</sup>+y<sup>2</sup>−2z<sup>2</sup>)−V(x<sup>2</sup>−2y<sup>2</sup>+z<sup>2</sup>)−4Mxy−2Nxz−2Pyz−Qx−Ry−Sz−T=0. (3) </FORM>
A suitable ellipsoid fitting the data is then found in Turner by finding the least squares solution to the linear equation<FORM>Λ{right arrow over (s)}<sub>LS</sub>={right arrow over (e)}, (4) </FORM>
where {right arrow over (s)}<sub>LS</sub>=(U,V,M,N,P,Q,R,S,T)<sup>T</sup>, the ith element of {right arrow over (e)} is x<sub>i</sub><sup>2</sup>+y<sub>i</sub><sup>2</sup>+z<sub>i</sub><sup>2 </sup>and the ith row of the m×9 matrix Λ is<FORM>x<sub>i</sub><sup>2</sup>+y<sub>i</sub><sup>2</sup>−2z<sub>i</sub><sup>2</sup>,x<sub>i</sub><sup>2</sup>−2y<sub>i</sub><sup>2</sup>+z<sub>i</sub><sup>2</sup>,4x<sub>i</sub>y<sub>i</sub>,2x<sub>i</sub>z<sub>i</sub>,2y<sub>i</sub>z<sub>i</sub>,x<sub>i</sub>,y<sub>i</sub>,z<sub>i</sub>,1. (5) </FORM>
Now an arbitrary ellipsoid can also be written as<FORM>({right arrow over (x)}−{right arrow over (b)})<sup>T</sup>A({right arrow over (x)}−{right arrow over (b)})=c, (6) </FORM>
where A is a symmetric positive definite matrix (and {right arrow over (x)}=(x,y,z), e.g., 3D field samples 130 from a 3D field sensor 110) and {right arrow over (b)} is a bias and c is a scaling factor. Expand this to get<FORM>{right arrow over (x)}<sup>T</sup>A{right arrow over (x)}−2{right arrow over (b)}<sup>T</sup>A{right arrow over (x)}+{right arrow over (b)}<sup>T</sup>A{right arrow over (b)}−c=0. (7) </FORM>
Expand further in terms of the elements of {right arrow over (x)} and {right arrow over (b)}, equate the left hand side (LHS) to the LHS of Eq. (3) and collect terms to get<maths id="MATHUS00002" num="2"><math overflow="scroll"><mtable><mtr><mtd><mrow><mrow><mi>A</mi><mo>=</mo><mrow><mo>(</mo><mtable><mtr><mtd><mrow><mn>1</mn><mo></mo><mi>U</mi><mo></mo><mi>V</mi></mrow></mtd><mtd><mrow><mrow><mo></mo><mn>2</mn></mrow><mo></mo><mi>M</mi></mrow></mtd><mtd><mrow><mo></mo><mi>N</mi></mrow></mtd></mtr><mtr><mtd><mrow><mrow><mo></mo><mn>2</mn></mrow><mo></mo><mi>M</mi></mrow></mtd><mtd><mrow><mn>1</mn><mo></mo><mi>U</mi><mo>+</mo><mrow><mn>2</mn><mo></mo><mi>V</mi></mrow></mrow></mtd><mtd><mrow><mo></mo><mi>P</mi></mrow></mtd></mtr><mtr><mtd><mrow><mo></mo><mi>N</mi></mrow></mtd><mtd><mrow><mo></mo><mi>P</mi></mrow></mtd><mtd><mrow><mn>1</mn><mo>+</mo><mrow><mn>2</mn><mo></mo><mi>U</mi></mrow><mo></mo><mi>V</mi></mrow></mtd></mtr></mtable><mo>)</mo></mrow></mrow><mo>,</mo></mrow></mtd><mtd><mrow><mo>(</mo><mn>8</mn><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mrow><mrow><mover><mi>b</mi><mo>></mo></mover><mo>=</mo><mrow><mfrac><mn>1</mn><mn>2</mn></mfrac><mo></mo><mrow><msup><mi>A</mi><mrow><mo></mo><mn>1</mn></mrow></msup><mo></mo><mrow><mo>(</mo><mtable><mtr><mtd><mi>Q</mi></mtd></mtr><mtr><mtd><mi>R</mi></mtd></mtr><mtr><mtd><mi>S</mi></mtd></mtr></mtable><mo>)</mo></mrow></mrow></mrow></mrow><mo>,</mo></mrow></mtd><mtd><mrow><mo>(</mo><mn>9</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>
and the scaling factor is<FORM>c=T+{right arrow over (b)}<sup>T</sup>A{right arrow over (b)}. (10) </FORM>
Observe that equation (6) has more parameters than strictly necessary to define an arbitrary ellipsoid. From Eq. (8), we can see that in Eq. (4) we have restricted ourselves to those ellipsoids for which tr(A)=3, where tr(A) means trace of the matrix A. The equations (8)(10) are invertible when matrix A satisfies this condition, so the least squares problem of Eq. (4) is then equivalent to<maths id="MATHUS00003" num="3"><math overflow="scroll"><mtable><mtr><mtd><mrow><mrow><mi>min</mi><mo></mo><mrow><munder><mo>∑</mo><mi>i</mi></munder><mo></mo><msup><mrow><mo>(</mo><mrow><mrow><msup><mrow><mo>(</mo><mrow><msub><mover><mi>x</mi><mo>></mo></mover><mi>i</mi></msub><mo></mo><mover><mi>b</mi><mo>></mo></mover></mrow><mo>)</mo></mrow><mi>T</mi></msup><mo></mo><mrow><mi>A</mi><mo></mo><mrow><mo>(</mo><mrow><msub><mover><mi>x</mi><mo>></mo></mover><mi>i</mi></msub><mo></mo><mover><mi>b</mi><mo>></mo></mover></mrow><mo>)</mo></mrow></mrow></mrow><mo></mo><mi>c</mi></mrow><mo>)</mo></mrow><mn>2</mn></msup></mrow></mrow><mo>,</mo></mrow></mtd><mtd><mrow><mo>(</mo><mn>11</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>
restricted to those A for which tr(A)=3. Now let A*,{right arrow over (b)}*, and c* be a solution to Eq. (11). Notice that we can make this an unrestricted problem by letting α<sub>11</sub>=3−α<sub>22</sub>−<sub>33</sub>. We can see that by scaling A* and c* by a constant k we get a solution for the optimization problem Eq. (11) restricted to the set where tr(A)=3/k. This is because<maths id="MATHUS00004" num="4"><math overflow="scroll"><mtable><mtr><mtd><mrow><mrow><mrow><mrow><msup><mrow><mo>(</mo><mrow><mover><mi>x</mi><mo>></mo></mover><mo></mo><mover><mi>b</mi><mo>></mo></mover></mrow><mo>)</mo></mrow><mi>T</mi></msup><mo></mo><mfrac><mi>A</mi><mi>k</mi></mfrac><mo></mo><mrow><mo>(</mo><mrow><mover><mi>x</mi><mo>></mo></mover><mo></mo><mover><mi>b</mi><mo>></mo></mover></mrow><mo>)</mo></mrow></mrow><mo></mo><mfrac><mi>c</mi><mi>k</mi></mfrac></mrow><mo>=</mo><mrow><mfrac><mn>1</mn><mi>k</mi></mfrac><mo></mo><mrow><mo>(</mo><mrow><mrow><msup><mrow><mo>(</mo><mrow><mover><mi>x</mi><mo>></mo></mover><mo></mo><mover><mi>b</mi><mo>></mo></mover></mrow><mo>)</mo></mrow><mi>T</mi></msup><mo></mo><mrow><mi>A</mi><mo></mo><mrow><mo>(</mo><mrow><mover><mi>x</mi><mo>></mo></mover><mo></mo><mover><mi>b</mi><mo>></mo></mover></mrow><mo>)</mo></mrow></mrow></mrow><mo></mo><mi>c</mi></mrow><mo>)</mo></mrow></mrow></mrow><mo>,</mo></mrow></mtd><mtd><mrow><mo>(</mo><mn>12</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>
so the partial derivatives (after making the problem unrestricted) are also just scaled, but at optimum the partial derivatives are zero, so the solutions are unaffected by this scaling.
In particular, we can select k=c* and get a solution (B*,b*,1) where<maths id="MATHUS00005" num="5"><math overflow="scroll"><mrow><mrow><msup><mi>B</mi><mo>*</mo></msup><mo>=</mo><mrow><mfrac><mn>1</mn><msup><mi>c</mi><mo>*</mo></msup></mfrac><mo></mo><msup><mi>A</mi><mo>*</mo></msup></mrow></mrow><mo>,</mo></mrow></math></maths>to<maths id="MATHUS00006" num="6"><math overflow="scroll"><mtable><mtr><mtd><mrow><mrow><mi>min</mi><mo></mo><mrow><munder><mo>∑</mo><mi>i</mi></munder><mo></mo><msup><mrow><mo>(</mo><mrow><mrow><msup><mrow><mo>(</mo><mrow><msub><mover><mi>x</mi><mo>></mo></mover><mi>i</mi></msub><mo></mo><mover><mi>b</mi><mo>></mo></mover></mrow><mo>)</mo></mrow><mi>T</mi></msup><mo></mo><mrow><msup><mi>B</mi><mo>*</mo></msup><mo></mo><mrow><mo>(</mo><mrow><msub><mover><mi>x</mi><mo>></mo></mover><mi>i</mi></msub><mo></mo><mover><mi>b</mi><mo>></mo></mover></mrow><mo>)</mo></mrow></mrow></mrow><mo></mo><mi>c</mi></mrow><mo>)</mo></mrow><mn>2</mn></msup></mrow></mrow><mo>,</mo></mrow></mtd><mtd><mrow><mo>(</mo><mn>13</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>
constrained to matrices B* such that<maths id="MATHUS00007" num="7"><math overflow="scroll"><mrow><mrow><mi>tr</mi><mo></mo><mrow><mo>(</mo><msup><mi>B</mi><mo>*</mo></msup><mo>)</mo></mrow></mrow><mo>=</mo><mrow><mfrac><mn>3</mn><msup><mi>c</mi><mo>*</mo></msup></mfrac><mo></mo><mrow><mi>A</mi><mo>.</mo></mrow></mrow></mrow></math></maths>In general, this may not be a solution for the problem in Eq. (11) constrained to c=1 instead of tr(A)=3. We call the problem constrained to c=1 the “transformation problem” herein. We would really want to find the solution of the transformation problem, but this problem generally cannot be linearized as above. We have, however, obtained a point that satisfies this constraint (e.g., of c=1) and is optimal in the sense that the solution cannot be improved without changing the trace of the matrix. In practice this will give good enough solutions to be used in calibration.
Note that equations (11) and (13) are identical, save that the matrix B in equation (13) replaces the matrix A in equation (11). The substantive difference lies in the scale. Solving equation (13) with the constraint<maths id="MATHUS00008" num="8"><math overflow="scroll"><mrow><mrow><mrow><mi>tr</mi><mo></mo><mrow><mo>(</mo><msup><mi>B</mi><mo>*</mo></msup><mo>)</mo></mrow></mrow><mo>=</mo><mfrac><mn>3</mn><msup><mi>c</mi><mo>*</mo></msup></mfrac></mrow><mo>,</mo></mrow></math></maths>where c* is from the solution to equation (11) with tr(A)=3, maps the samples close to the unit sphere. The solution for equation (11) maps the samples close to a sphere, just not necessarily the unit sphere. Solving equations (4) and (8) to (10) yields a transformation that draws the samples close to a unit sphere. Consider that the transformation obtained by taking A* by decomposing takes the samples close to some spherical surface, presumably having a radius in the range of √{square root over (c*)}. Scaling A* by c* yields a transformation, after decomposing, that takes the samples close the spherical surface with unit radius.
It is noted that the scaled solution of equation (11), as represented by equation (13), is also a direct solution to a similar problem, but using an alternate trace condition. Intuitively, this confirms that scaling as in equation (13) maintains whatever optimality is imposed by the more general, unscaled solution of equation (11). Since c* in the equation (13) condition of<maths id="MATHUS00009" num="9"><math overflow="scroll"><mrow><mrow><mi>tr</mi><mo></mo><mrow><mo>(</mo><mi>B</mi><mo>)</mo></mrow></mrow><mo>=</mo><mfrac><mn>3</mn><msup><mi>c</mi><mo>*</mo></msup></mfrac></mrow></math></maths>comes from the solution to equation (11), it is clear that equation (13) is but one of several possible scaling solutions, any of which scale for and solve equation (11) remaining within these teachings.
Now decompose B* so that B*=C<sup>T</sup>C and note that<FORM>({right arrow over (x)}−{right arrow over (c)})<sup>T</sup>B({right arrow over (x)}−{right arrow over (c)})=(x−c)<sup>T</sup>C<sup>T</sup>C({right arrow over (x)}−{right arrow over (c)})=∥C({right arrow over (x)}−{right arrow over (c)})∥<sup>2</sup>, (14) </FORM>
we see that in fact we have found (step 330 of FIG. 3) the transformation (i.e., the matrix C) that maps the data points close to a unit sphere. The transformation matrix C may be applied to samples (in step 340 of FIG. 3) to create transformed samples. However, as pointed out below, it may be beneficial to determine a different version of the transformation matrix C to be applied to samples.
It should be noted that it is possible to solve Eq. (4) directly, which gives the parameters of {right arrow over (s)}<sub>LS </sub>(e.g., step 320 of FIG. 3), and then to use the equations above to determine the transformation matrix C (step 330 of FIG. 3) and to apply the transformation matrix C to 3D field samples 130. However, a recursive method is disclosed in the section below and is more suitable to determining these elements based on realtime, incoming 3D field samples 130.
4. Recursive Updates
The solution to Eq. (4) is well known and is<FORM>{right arrow over (s)}<sub>LS</sub>=(Λ<sup>T</sup>Λ)<sup>−1</sup>Λ<sup>T</sup>{right arrow over (e)}. (15) </FORM>
This solution corresponds to the solution of<maths id="MATHUS00010" num="10"><math overflow="scroll"><mtable><mtr><mtd><mrow><msub><mover><mi>s</mi><mo>></mo></mover><mi>LS</mi></msub><mo>=</mo><mrow><munder><mi>min</mi><msub><mover><mi>s</mi><mo>></mo></mover><mi>LS</mi></msub></munder><mo></mo><mrow><munder><mo>∑</mo><mi>i</mi></munder><mo></mo><mrow><msup><mrow><mo>(</mo><mrow><mrow><msub><mi>Λ</mi><mi>i</mi></msub><mo></mo><msub><mover><mi>s</mi><mo>></mo></mover><mi>LS</mi></msub></mrow><mo></mo><msub><mi>e</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow><mn>2</mn></msup><mo>.</mo></mrow></mrow></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>16</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>
where Λ<sub>i </sub>is the ith row of Λ. In the spirit of the Recursive Least Squares (RLS) algorithm, a forgetting factor ε∈(0,1) is introduced, and the following is minimized:<maths id="MATHUS00011" num="11"><math overflow="scroll"><mtable><mtr><mtd><mtable><mtr><mtd><mrow><msub><mover><mi>s</mi><mo>></mo></mover><mi>LS</mi></msub><mo>=</mo><mrow><munder><mi>min</mi><msub><mover><mi>s</mi><mo>></mo></mover><mi>LS</mi></msub></munder><mo></mo><mrow><munder><mo>∑</mo><mi>i</mi></munder><mo></mo><msup><mrow><msup><mi>ρ</mi><mfrac><mrow><mi>m</mi><mo></mo><mi>i</mi></mrow><mn>2</mn></mfrac></msup><mo></mo><mrow><mo>(</mo><mrow><mrow><msub><mi>Λ</mi><mi>i</mi></msub><mo></mo><msub><mover><mi>s</mi><mo>></mo></mover><mi>LS</mi></msub></mrow><mo></mo><msub><mi>e</mi><mi>i</mi></msub></mrow><mo>)</mo></mrow></mrow><mn>2</mn></msup></mrow></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mo>=</mo><mrow><munder><mi>min</mi><msub><mover><mi>s</mi><mo>></mo></mover><mi>LS</mi></msub></munder><mo></mo><mrow><munder><mo>∑</mo><mi>i</mi></munder><mo></mo><msup><mrow><mo>(</mo><mrow><mrow><msup><mi>ρ</mi><mrow><mi>m</mi><mo></mo><mi>i</mi></mrow></msup><mo></mo><msub><mi>Λ</mi><mi>i</mi></msub><mo></mo><msub><mover><mi>s</mi><mo>></mo></mover><mi>LS</mi></msub></mrow><mo></mo><mrow><msub><mi>e</mi><mi>i</mi></msub><mo></mo><msup><mi>ρ</mi><mrow><mi>m</mi><mo></mo><mi>i</mi></mrow></msup></mrow></mrow><mo>)</mo></mrow><mn>2</mn></msup></mrow></mrow></mrow></mtd></mtr></mtable></mtd><mtd><mrow><mo>(</mo><mn>17</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>
It is then clear that we can simply weight the rows of Λ and the elements of {right arrow over (e)}, by ρ<sup>mi</sup>. So let<maths id="MATHUS00012" num="12"><math overflow="scroll"><mtable><mtr><mtd><mrow><mrow><msub><mi>Λ</mi><msup><mi>ρ</mi><mi>m</mi></msup></msub><mo>=</mo><mrow><mo>(</mo><mtable><mtr><mtd><mrow><msup><mi>ρ</mi><mfrac><mi>m</mi><mn>2</mn></mfrac></msup><mo></mo><msub><mi>Λ</mi><mn>1</mn></msub></mrow></mtd></mtr><mtr><mtd><mrow><msup><mi>ρ</mi><mfrac><mrow><mi>m</mi><mo></mo><mn>1</mn></mrow><mn>2</mn></mfrac></msup><mo></mo><msub><mi>Λ</mi><mn>2</mn></msub></mrow></mtd></mtr><mtr><mtd><mi>⋮</mi></mtd></mtr><mtr><mtd><msub><mi>Λ</mi><mi>m</mi></msub></mtd></mtr></mtable><mo>)</mo></mrow></mrow><mo>,</mo><mstyle><mtext></mtext></mstyle><mo></mo><mi>and</mi></mrow></mtd><mtd><mrow><mo>(</mo><mn>18</mn><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mrow><mrow><msubsup><mover><mi>e</mi><mo>></mo></mover><mi>ρ</mi><mi>m</mi></msubsup><mo>=</mo><msup><mrow><mo>(</mo><mrow><mrow><msub><mi>e</mi><mn>1</mn></msub><mo></mo><msup><mi>ρ</mi><mfrac><mi>m</mi><mn>2</mn></mfrac></msup></mrow><mo>,</mo><mrow><msub><mi>e</mi><mn>2</mn></msub><mo></mo><msup><mi>ρ</mi><mfrac><mrow><mi>m</mi><mo></mo><mn>1</mn></mrow><mn>2</mn></mfrac></msup></mrow><mo>,</mo><mi>…</mi><mo></mo><mstyle><mtext> </mtext></mstyle><mo>,</mo><msub><mi>e</mi><mi>m</mi></msub></mrow><mo>)</mo></mrow><mi>T</mi></msup></mrow><mo>,</mo></mrow></mtd><mtd><mrow><mo>(</mo><mn>19</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>
Then let<maths id="MATHUS00013" num="13"><math overflow="scroll"><mtable><mtr><mtd><mrow><mrow><msub><mi>Ξ</mi><mi>m</mi></msub><mo>=</mo><mstyle><mtext> </mtext></mstyle><mo></mo><mrow><msup><mrow><mo>(</mo><mrow><msubsup><mi>Λ</mi><msup><mi>ρ</mi><mi>m</mi></msup><mi>T</mi></msubsup><mo></mo><msub><mi>Λ</mi><msup><mi>ρ</mi><mi>m</mi></msup></msub></mrow><mo>)</mo></mrow><mrow><mo></mo><mn>1</mn></mrow></msup><mo>=</mo><mstyle><mtext> </mtext></mstyle><mo></mo><msup><mrow><mo>(</mo><mrow><munder><mo>∑</mo><mi>i</mi></munder><mo></mo><mrow><msup><mi>ρ</mi><mrow><mi>m</mi><mo></mo><mi>i</mi></mrow></msup><mo></mo><msubsup><mi>Λ</mi><mi>i</mi><mi>T</mi></msubsup><mo></mo><msub><mi>Λ</mi><mi>i</mi></msub></mrow></mrow><mo>)</mo></mrow><mrow><mo></mo><mn>1</mn></mrow></msup></mrow></mrow><mo>,</mo><mrow><mo>=</mo><mstyle><mtext> </mtext></mstyle><mo></mo><msup><mrow><mo>(</mo><mrow><msub><mi>ρΞ</mi><mrow><mi>m</mi><mo></mo><mn>1</mn></mrow></msub><mo>+</mo><mrow><msubsup><mi>Λ</mi><mi>i</mi><mi>T</mi></msubsup><mo></mo><msub><mi>Λ</mi><mi>i</mi></msub></mrow></mrow><mo>)</mo></mrow><mrow><mo></mo><mn>1</mn></mrow></msup></mrow><mo></mo><mstyle><mtext> </mtext></mstyle><mo>,</mo><mstyle><mtext></mtext></mstyle><mo></mo><mi>and</mi></mrow></mtd><mtd><mrow><mo>(</mo><mn>20</mn><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mrow><msub><mover><mi>ξ</mi><mo>></mo></mover><mi>m</mi></msub><mo>=</mo><mrow><msub><mi>Λ</mi><msup><mi>ρ</mi><mi>m</mi></msup></msub><mo></mo><mrow><msubsup><mover><mi>e</mi><mo>></mo></mover><mi>ρ</mi><mi>m</mi></msubsup><mo>.</mo></mrow></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>21</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>
Using the Woodbury matrix identity, we get a well known update for the least squares system, which in our notation is<maths id="MATHUS00014" num="14"><math overflow="scroll"><mtable><mtr><mtd><mrow><mrow><msub><mover><mi>ζ</mi><mo>></mo></mover><mi>m</mi></msub><mo>=</mo><mrow><mrow><mi>ρ</mi><mo></mo><msub><mover><mi>ξ</mi><mo>></mo></mover><mrow><mi>m</mi><mo></mo><mn>1</mn></mrow></msub></mrow><mo>+</mo><mrow><msub><mi>e</mi><mi>m</mi></msub><mo></mo><msubsup><mi>Λ</mi><mi>m</mi><mi>T</mi></msubsup></mrow></mrow></mrow><mo>,</mo></mrow></mtd><mtd><mrow><mo>(</mo><mn>22</mn><mo>)</mo></mrow></mtd></mtr><mtr><mtd><mrow><msub><mi>Ξ</mi><mi>m</mi></msub><mo>=</mo><mrow><mfrac><mn>1</mn><mi>ρ</mi></mfrac><mo></mo><mrow><mrow><msub><mi>Ξ</mi><mrow><mi>m</mi><mo></mo><mn>1</mn></mrow></msub><mo></mo><mrow><mo>(</mo><mrow><mi>I</mi><mo></mo><mrow><mfrac><mrow><msubsup><mi>Λ</mi><mi>m</mi><mi>T</mi></msubsup><mo></mo><msub><mi>Λ</mi><mi>m</mi></msub></mrow><mrow><mi>ρ</mi><mo>+</mo><mrow><msub><mi>Λ</mi><mi>m</mi></msub><mo></mo><msub><mi>Ξ</mi><mrow><mi>m</mi><mo></mo><mn>1</mn></mrow></msub><mo></mo><msubsup><mi>Λ</mi><mi>m</mi><mi>T</mi></msubsup></mrow></mrow></mfrac><mo></mo><msub><mi>Ξ</mi><mrow><mi>m</mi><mo></mo><mn>1</mn></mrow></msub></mrow></mrow><mo>)</mo></mrow></mrow><mo>.</mo></mrow></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>23</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>
A major advantage of this update being, of course, that no matrix inversions are required.
The parameters for the ellipsoid fit are then obtained from<FORM>{right arrow over (s)}<sub>LS</sub>=Ξ<sub>m</sub>{right arrow over (ξ)}<sub>m</sub>. (24) </FORM>
As shown above, the chosen parameters for the ellipsoid are selected such that the ellipsoid is the “best” fit to the 3D field samples 130. Earlier we saw that the linear least squares problem differs from the sphere transformation problem essentially in the way different sample points are weighted. The introduction of the forgetting factor, ρ, then further obscures the impact of the linearization to the optimality of the solution.
5. Matrix Decomposition
Typically, the parameters for the ellipsoid fit are determined using Eq. 24. Then, the matrix B and bias vector {right arrow over (b)} are obtained from equations (8)(10) and from<maths id="MATHUS00015" num="15"><math overflow="scroll"><mrow><mi>B</mi><mo>=</mo><mrow><mfrac><mn>1</mn><mi>c</mi></mfrac><mo></mo><mrow><mi>A</mi><mo>.</mo></mrow></mrow></mrow></math></maths>If the data is well distributed, this matrix, B, will be positive definite. Although the transformation matrix C from Eq. (14) is one possible solution to a transformation matrix, it turns out that there are innumerable matrices D such that D<sup>T</sup>D=B. This is easy to see since any rotation matrix R satisfies R<sup>T</sup>R=I and (RD)<sup>T</sup>(RD)=D<sup>T</sup>R<sup>T</sup>RD=D<sup>T</sup>D=B. Out of these possible matrices D, we would like to select one that is somehow “fixed” and that may be used for transformation.
One suitable technique for fixing D is to select a matrix D such that D{circumflex over (x)}=λ{circumflex over (x)}, where λ is an arbitrary constant and {circumflex over (x)}=(1,0,0)<sup>T</sup>. This way, the xaxis of the sensor is never rotated by the calibration. We know then that d<sub>21</sub>=d<sub>31</sub>=0. Also, since D<sup>T</sup>D=B we know that d<sub>11</sub>=√{square root over (b<sub>11</sub>)}. So we have that the first column of D is<maths id="MATHUS00016" num="16"><math overflow="scroll"><mtable><mtr><mtd><mrow><msub><mover><mi>d</mi><mo>></mo></mover><mn>1</mn></msub><mo>=</mo><mrow><mrow><mo>(</mo><mtable><mtr><mtd><msqrt><msub><mi>b</mi><mn>11</mn></msub></msqrt></mtd></mtr><mtr><mtd><mn>0</mn></mtd></mtr><mtr><mtd><mn>0</mn></mtd></mtr></mtable><mo>)</mo></mrow><mo>.</mo></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>25</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>
The second column of D, {right arrow over (d)}<sub>2</sub>, satisfies that {right arrow over (d)}<sub>1</sub><sup>T</sup>{right arrow over (d)}<sub>2</sub>=b<sub>12 </sub>and ∥{right arrow over (d)}<sub>2</sub>∥<sup>2</sup>=b<sub>22</sub>. These alone do not uniquely determine {right arrow over (d)}<sub>2</sub>, so we also require in this example that the last element of {right arrow over (d)}<sub>2 </sub>be zero. This way, any rotation to the yaxis of the sensor is constrained to the x,yplane. These conditions are satisfied by the vector<maths id="MATHUS00017" num="17"><math overflow="scroll"><mtable><mtr><mtd><mrow><msub><mover><mi>d</mi><mo>></mo></mover><mn>2</mn></msub><mo>=</mo><mrow><msqrt><msub><mi>b</mi><mn>22</mn></msub></msqrt><mo></mo><mrow><mrow><mo>(</mo><mtable><mtr><mtd><mfrac><msub><mi>b</mi><mn>12</mn></msub><msqrt><mrow><msub><mi>b</mi><mn>11</mn></msub><mo></mo><msub><mi>b</mi><mn>12</mn></msub></mrow></msqrt></mfrac></mtd></mtr><mtr><mtd><msqrt><mrow><mn>1</mn><mo></mo><mfrac><msubsup><mi>b</mi><mn>12</mn><mn>2</mn></msubsup><mrow><msub><mi>b</mi><mn>11</mn></msub><mo></mo><msub><mi>b</mi><mn>22</mn></msub></mrow></mfrac></mrow></msqrt></mtd></mtr><mtr><mtd><mn>0</mn></mtd></mtr></mtable><mo>)</mo></mrow><mo>.</mo></mrow></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>26</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>
Finally, the third column is now subject to the constraints {right arrow over (d)}<sub>1</sub><sup>T</sup>{right arrow over (d)}<sub>3</sub>=b<sub>13</sub>, {right arrow over (d)}<sub>2</sub><sup>T</sup>d<sub>3</sub>=b<sub>23</sub>, and ∥{right arrow over (d)}<sub>3</sub>∥<sup>2</sup>=b<sub>33</sub>.
Let<FORM>{right arrow over (d)}<sub>3</sub>=α{right arrow over (d)}<sub>1</sub>+β{right arrow over (d)}<sub>2</sub>+γ{circumflex over (z)}, (27) </FORM>
where {circumflex over (z)}=(0,0,1)<sup>T</sup>. Then<FORM>b<sub>13</sub>={right arrow over (d)}<sub>1</sub><sup>T</sup>{right arrow over (d)}<sub>3</sub>=αb<sub>11</sub>+βb<sub>12</sub>, (28) </FORM><FORM>b<sub>23</sub>={right arrow over (d)}<sub>2</sub><sup>T</sup>{right arrow over (d)}<sub>3</sub>=ab<sub>12</sub>+βb<sub>22</sub>. (29) </FORM>
This system of linear equations will give α and β when β is positive definite (because the determinant b<sub>11</sub>b<sub>22</sub>b<sub>12</sub><sup>2 </sup>has to be greater than zero in this example for positive definite matrices). Finally we get<maths id="MATHUS00018" num="18"><math overflow="scroll"><mtable><mtr><mtd><mrow><msub><mover><mi>d</mi><mo>></mo></mover><mn>3</mn></msub><mo>=</mo><mrow><mrow><mo>(</mo><mtable><mtr><mtd><mrow><mrow><mi>α</mi><mo></mo><msqrt><msub><mi>b</mi><mn>11</mn></msub></msqrt></mrow><mo>+</mo><mrow><mi>β</mi><mo></mo><mfrac><mrow><msqrt><msub><mi>b</mi><mn>22</mn></msub></msqrt><mo></mo><msub><mi>b</mi><mn>12</mn></msub></mrow><msqrt><mrow><msub><mi>b</mi><mn>11</mn></msub><mo></mo><msub><mi>b</mi><mn>22</mn></msub></mrow></msqrt></mfrac></mrow></mrow></mtd></mtr><mtr><mtd><mrow><mi>β</mi><mo></mo><msqrt><mrow><msub><mi>b</mi><mn>22</mn></msub><mo></mo><mfrac><mrow><msub><mi>b</mi><mn>22</mn></msub><mo></mo><msubsup><mi>b</mi><mn>12</mn><mn>2</mn></msubsup></mrow><mrow><msub><mi>b</mi><mn>11</mn></msub><mo></mo><msub><mi>b</mi><mn>22</mn></msub></mrow></mfrac></mrow></msqrt></mrow></mtd></mtr><mtr><mtd><msqrt><mtable><mtr><mtd><mrow><msub><mi>b</mi><mn>33</mn></msub><mo></mo><msup><mrow><mo>(</mo><mrow><mrow><mi>α</mi><mo></mo><msqrt><msub><mi>b</mi><mn>11</mn></msub></msqrt></mrow><mo>+</mo><mrow><mi>β</mi><mo></mo><mfrac><mrow><msqrt><msub><mi>b</mi><mn>22</mn></msub></msqrt><mo></mo><msub><mi>b</mi><mn>12</mn></msub></mrow><msqrt><mrow><msub><mi>b</mi><mn>11</mn></msub><mo></mo><msub><mi>b</mi><mn>22</mn></msub></mrow></msqrt></mfrac></mrow></mrow><mo>)</mo></mrow><mn>2</mn></msup><mo></mo></mrow></mtd></mtr><mtr><mtd><msup><mrow><mo>(</mo><mrow><mi>β</mi><mo></mo><msqrt><mrow><msub><mi>b</mi><mn>22</mn></msub><mo></mo><mfrac><mrow><msub><mi>b</mi><mn>22</mn></msub><mo></mo><msubsup><mi>b</mi><mn>12</mn><mn>2</mn></msubsup></mrow><mrow><msub><mi>b</mi><mn>11</mn></msub><mo></mo><msub><mi>b</mi><mn>22</mn></msub></mrow></mfrac></mrow></msqrt></mrow><mo>)</mo></mrow><mn>2</mn></msup></mtd></mtr></mtable></msqrt></mtd></mtr></mtable><mo>)</mo></mrow><mo>.</mo></mrow></mrow></mtd><mtd><mrow><mo>(</mo><mn>30</mn><mo>)</mo></mrow></mtd></mtr></mtable></math></maths>
The third element of {right arrow over (d)}<sub>3 </sub>is simply selected such that d<sub>13</sub><sup>2</sup>+d<sub>23</sub><sup>2</sup>+d<sub>33</sub><sup>2</sup>=b<sub>33</sub>, and is simply √{square root over (b<sub>33</sub>−d<sub>13</sub><sup>2</sup>−d<sub>23</sub><sup>2</sup>)}. The matrix D=└{right arrow over (d)}<sub>1</sub>,{right arrow over (d)}<sub>2</sub>,{right arrow over (d)}<sub>3</sub>┘ now satisfies that D<sup>T</sup>D=B. The matrix D can therefore be used as a transformation matrix and applied to samples in step 340 of FIG. 3.
It should be noted that in some conventional systems, the final calibration is separated into multiple matrices, one of which is used to remove linear dependence among the different axes. That matrix is called the decoupling matrix. In an exemplary embodiment of the disclosed invention, all the matrices are combined and thus the decoupling effect is included in the transformation matrix D. However, for instance in satellite systems, this decoupling portion may not be available exactly because the outside magnetic field changes during calibration. Thus, in other embodiments, the transformation matrix D can be modified to include multiple matrices, one of which would be a decoupling matrix.
An alternative technique for fixing D is to find the eigenvalue decomposition of B such that B=U*ΣU, where Σ is a diagonal matrix whose elements are the eigenvalues of B and the columns of U are the eigenvectors of B. Since B is symmetric and positive definite (i.e., unless the data are degenerate) the eigenvalues and eigenvectors are real, so U*=U<sup>T</sup>. Since Σ<sup>1/2 </sup>is easy to compute, we would select D=Σ<sup>1/2</sup>U. A (rotation) matrix R, such that R<sup>T</sup>R=I, can be found such that RD is close to the identity matrix. One way to do this is to let R=D<sup>−1 </sup>and then orthonormalize the rows of R, such as with the GramSchmidt method. The GramSchmidt method is guaranteed to not change the direction of the first processed row, and thus the corresponding row of RD is parallel to the corresponding coordinate axis. Said another way, if we start from the top, the first row of RD is parallel to the xcoordinate axis. Of course, some other orthonormalization method besides GramSchmidt may be used to guarantee the direction of a selected column of RD. In either of the above approaches to obtain R, the three column vectors of R are also orthonormalized.
It follows from the above derivation of R that not all rows will be orthogonal to all columns, because for a row of R and its corresponding column of D the vectors should be parallel (e.g., so that the diagonal elements of RD would be near the value of one and the offdiagonal elements near to the value zero within the constraint that R represents a rotoinversion. The transformation matrix applied (see step 340 of FIG. 3) to 3D field samples 130 would then be RD. The first more direct technique is superior to the eigenvalue decomposition technique in many ways. Illustratively, the first technique is numerically more stable than most eigenvalue decomposition methods and likely to be computationally cheaper.
6. Convergence Considerations
The update formulas in Section 4 are generally guaranteed to give the best fitting ellipsoid with respect to the residuals minimized. However, when the 3D field samples 130 are degenerate, for example in a plane, the matrix B may not be invertible and the techniques above may fail to produce a meaningful transformation matrix D.
There should be some method of determining which samples are fed to the autocalibration algorithm to reduce or eliminate the possibility of degenerate data. In an exemplary embodiment, orientation sensors such as gyroscopes and accelerometers may be relied upon to give an approximate idea of an orientation of an apparatus to help guarantee that 3D field samples 130 are well representative of an ellipsoid.
For instance, referring to FIG. 4, an exemplary apparatus 400 for calibrating a 3D field sensor 110 is shown. As in FIG. 2, the 3D field sensor 110 produces 3D field samples 130 and the calibration module 140 produces transformed 3D samples 150 using the 3D field samples 130. In the embodiment shown in FIG. 4, however, an orientation sensor 410 produces orientation samples 420. A filtering module 430 uses the orientation samples 420 to determine selected 3D field samples 440. Instead of using each 3D field sample 130 to perform calibration, the calibration module 140 uses only the selected 3D field samples 440 for calibration of the 3D field sensor 110. The filtering module 430 filters out some of the 3D field samples 130 so that the selected 3D field samples 440 that are passed to the calibration module 140 will decrease the likelihood of degenerate data used by the calibration module 140 during calibration.
The orientation sensor 410 can be one or several sensors that give orientation in three directions. For instance, one or more (typically three) gyroscopes or one or more (typically three) accelerometers may be used. Although the filtering module 430 is shown separately from the calibration module 140, portions or all of the filtering module 430 may be combined with the calibration module 140.
Turning now to FIG. 5, a flowchart is shown of an exemplary method 500 for filtering samples from a 3D field sensor 110. Method 500 is performed in this example by the filtering module 430. In step 510, an orientation sample 420 is accessed from the 3D field sensor 510. The orientation sample 420 is accessed by the filtering module 430 directly from the orientation sensor 410, but may also be accessed from memory (see FIG. 6, for instance), or through any other technique for accessing a sample. In step 520, an orientation is determined from one or more of the orientation samples 420. For a set of three gyroscopes, step 520 is performed by integrating the angular velocity (e.g., in an orientation samples 420) from the gyroscopes to get an orientation. In step 530, the orientation from the gyroscopes is quantized to one of the three possible directions. Similar techniques could be used for other orientation sensors 410. Each possible direction is associated with one of three “buckets”, which are used to hold selected 3D field samples 440 for calibration purposes. For calibration, all selected 3D field samples 440 in the buckets are used. For purposes of method 500, the buckets are used to reduce or prevent degenerate data.
In step 540, a current 3D field sample 130 (e.g., from a 3D field sensor 110 such as a magnetometer) is assigned to one of the buckets based on this quantized orientation. Each 3D field sample 130 is assigned to one of the three buckets, each of which represents one of the three orthogonal axes. In step 550, it is determined if the bucket into which the current 3D field sample 130 is placed has the most 3D field samples. A 3D field sample 130 that would fall into the bucket that already has the most samples (step 550=YES) is not used for calibration or added to the buckets (step 570). The 3D field sample 130 may be transformed, however, On the other hand, a 3D field sample 130 that does not fall into the bucket that already has the most samples (step 550=NO) is assigned to a bucket and used for calibration by the calibration module 140 (step 560). The 3D field samples 130 in the set of three buckets may be used for calibration and maintained in a memory, for instance.
It is beneficial at this point to review information about orientation sensors. Gyroscopes, for instance, can be used to track orientation of a device with some accuracy, but there is the problem of drift. A 3D magnetometer and 3D accelerometer together constitute a kind of orientation sensor. If the direction of gravitation and the earth's magnetic field are known relative to the device, then one can deduce the orientation of the device “relative to earth”. This is, of course, exactly what is being done herein with magnetometers and the reason why the magnetometers need to be calibrated somewhat accurately. The absolute orientation from the 3D magnetometer and 3D accelerometer can then be used to compensate for the drift in the gyroscopes. Gyroscopes on the other hand provide accurate tracking of fast movements, which is typically impossible with just the other two sensors. However, while a 3D accelerometer alone cannot give the orientation (e.g., because rotation around the axis of gravitation is invisible to the 3D accelerometer), a 3D accelerometer can give hints when the device surely has turned. So a 3D accelerometer would be another example of a sensor that could be used to select 3D field samples 440.
Although a gyroscope drifts, a gyroscope is still accurate enough to virtually guarantee that the data from the 3D field samples 130 will not be degenerate. If the range of the data (e.g., the 3D field samples 130 from the 3D field sensor 110) is well defined, it should be possible to devise a method to select the selected 3D field samples 440 for the calibration algorithm that would not need extra sensors in the orientation sensor 410.
7. Another Exemplary Embodiment
Turning now to FIG. 6, FIG. 6 is a block diagram of an apparatus 600 including a 3D field sensor 640 and modules (660, 670) for calibration of the 3D field sensor 110. In the example of FIG. 6, the apparatus 600 can be a mobile phone, a fitness monitoring device, a game controller, a wearable motion tracking device, or any device whose orientation in 3D space is of interest and which is rotated relatively frequently (e.g., enough to get an amount of 3D field samples 130 so that a ellipse can be fit to the 3D field samples 130).
Concerning the number 3D field samples 130 that can be used, in theory, one needs at least nine samples for the matrix Λ<sup>T</sup>Λ to be invertible. So nine samples are the minimum number of samples that can uniquely define an ellipsoid. In practice, one would probably need significantly more than that. In an exemplary implementation, a forgetting factor (e.g., ρ) of 0.95 was used, which results in a fairly fast adaptation. After 50 3D field samples 130 are used for calibration, the first sample has a corresponding weight of 0.077, so one would reason, then, that the first sample is practically completely forgotten after 50 samples. So the calibration is then basically based on the last 50 samples only, and 50 samples should certainly be enough. In an exemplary embodiment, the number of samples recommended would be somewhere between 10 and 50.
Apparatus 600 comprises a processor 610 that communicates with a memory 615, a 3D field sensor 620, and an orientation sensor 625. The memory 615 can be a static random access memory (RAM), dynamic RAM, registers, portions of a hard drive, or any memory that is readable and writable. Memory 615 comprises a calibration module 630, a filtering module 635, a current 3D field sample 640, orientation sample(s) 645, buckets of 3D field samples 650, a transformed 3D field sample 660, and a 3D application 665. The calibration module 630, filtering module 635, and 3D application 665 are implemented in this example as software (e.g., or firmware) comprising instructions and the instructions are loaded into memory 615 and executed by the processor 610. Processor 610 could be a general purpose processor, a digital signal processor, or any type of device suitable for executing instructions. The current 3D field sample 640 is transformed by the calibration module 630 into a transformed 3D field sample 660, which is used by 3D application 665. The 3D application 665 is any application that can use 3D information. The buckets of 3D field samples 650 contain three buckets (not shown), each of which contains selected 3D field samples, as described above in reference to FIGS. 4 and 5. The filtering module 635 is used as described above to filter 3D field samples for placement into the buckets of 3D field samples 650.
Portions or all of the filtering module 635 and the calibration module 630 could be combined, if desired. Additionally, the filtering module 635 and calibration module 630 could act on 3D field samples 130 and orientation samples 420 that have been previously determined and stored for later analysis.
Although being described in FIG. 6 as being software or firmware, calibration module 140 and filter module 430 may be implemented as software, firmware, hardware (e.g., discrete components that make up a circuit, an integrated circuit implemented on a semiconductor, and a programmable logic device), or some combination of these. It should be noted that the various blocks of the flow diagrams of FIGS. 3 and 5 might represent program steps, or interconnected logic circuits, blocks, and functions, or a combination of program steps and logic circuits, blocks and functions for performing the specified tasks. Furthermore, the disclosed invention may be implemented as a signal bearing medium tangibly embodying a program of machinereadable instructions executable by a processor to perform operations described herein.
The foregoing description has provided by way of exemplary and nonlimiting examples a full and informative description of the best method and apparatus presently contemplated by the inventors for carrying out the invention. However, various modifications and adaptations may become apparent to those skilled in the relevant arts in view of the foregoing description, when read in conjunction with the accompanying drawings and the appended claims. Nonetheless, all such and similar modifications of the teachings of this invention will still fall within the scope of this invention.
Furthermore, some of the features of the preferred embodiments of this invention could be used to advantage without the corresponding use of other features. As such, the foregoing description should be considered as merely illustrative of the principles of the present invention, and not in limitation thereof.