use crate::config::{ALIGNED_HEADING, Config};
use crate::frames::{Body, Ned};
use crate::gsf::YawEstimator;
use crate::math::wrap_pi;
use crate::observation::heading::{has_heading, heading_of};
use crate::propagate::ImuSample;
use crate::units::{Position, Radians, Timestamp, Velocity, VelocityNoise};
use nalgebra::Vector3;
use super::Eskf;
const YAW_SIGMA_MAX: Radians = Radians::from_degrees(15.0);
const YAW_FAILURE: Radians = Radians::from_degrees(25.0);
#[derive(Clone, Copy, Debug, Default, PartialEq, Eq)]
pub(super) struct YawReplaced {
position: bool,
velocity: bool,
magnetometer: bool,
}
impl YawReplaced {
const OWING: Self = Self {
position: true,
velocity: true,
magnetometer: true,
};
pub(super) fn owes_position(&self) -> bool {
self.position
}
pub(super) fn owes_velocity(&self) -> bool {
self.velocity
}
pub(super) fn overrules_magnetometer(&self) -> bool {
self.magnetometer
}
pub(super) fn settle_position(&mut self) {
self.position = false;
}
pub(super) fn settle_velocity(&mut self) {
self.velocity = false;
}
pub(super) fn settle_magnetometer(&mut self) {
self.magnetometer = false;
}
}
const _: () = assert!(YAW_SIGMA_MAX.as_radians() < ALIGNED_HEADING.as_radians());
impl Eskf {
#[inline(never)]
pub(super) fn step_yaw_estimator(&mut self, imu: ImuSample) {
let Config {
yaw_estimator: true,
gravity,
..
} = self.config
else {
return;
};
let gyro_bias = self.estimate.state().gyro_bias.vector();
self.yaw_estimator.predict(imu, gravity, gyro_bias);
}
pub(super) fn restart_yaw_estimator(&mut self) {
if self.config.yaw_estimator {
self.yaw_estimator.restart();
}
}
#[inline(never)]
pub(super) fn weigh_yaw(
&mut self,
time: Timestamp,
velocity: Velocity<Ned>,
noise: VelocityNoise<Ned>,
antenna: Position<Body>,
) {
if !(self.config.yaw_estimator && self.yaw_estimator.is_leveled()) {
return;
}
let [north, east, _] = noise.variance();
let Some(variance) = YawEstimator::variance_of(north, east) else {
return;
};
let Some(carried) = self.carried_velocity(Velocity::zero(), antenna, time) else {
return;
};
let heading = heading_of(self.estimate.state().attitude.quaternion());
self.diagnostics.yaw_estimator.note_arrival(time);
self.yaw_estimator.fuse_velocity(
time,
velocity.vector().xy(),
variance,
carried.vector().xy(),
heading,
self.config.gravity,
);
}
pub(super) fn adopt_first_yaw(&mut self) {
if !(self.unestablished.heading && self.yaw_estimator.is_settled()) {
return;
}
if self.turn_to_estimated_yaw(|_| true) {
self.diagnostics.yaw_estimator.record_adopted();
self.note_alignment();
self.magnetic_north = false;
}
}
#[inline(never)]
pub(super) fn replace_failed_yaw(&mut self) {
if !self.config.yaw_estimator || !self.yaw_estimator.is_settled() {
return;
}
let after = self.unless_accepted(
self.config.recovery.yaw_estimator,
&[&self.diagnostics.gnss_heading],
);
let since_initialized = self.diagnostics.since_initialized;
let velocity = &self.diagnostics.gnss_velocity;
if !velocity.locked_out(after, since_initialized) || velocity.consecutive_rejections == 0 {
return;
}
let disagrees = |y: f32| y.abs() > YAW_FAILURE.as_radians();
if self.turn_to_estimated_yaw(disagrees) {
self.diagnostics.yaw_estimator.record_recovered();
self.note_alignment();
self.magnetic_north = false;
self.yaw_replaced = YawReplaced::OWING;
}
}
fn turn_to_estimated_yaw(&mut self, wanted: impl FnOnce(f32) -> bool) -> bool {
let Some((yaw, variance)) = self.yaw_estimator.yaw() else {
return false;
};
let bar = YAW_SIGMA_MAX.as_radians();
let converged = variance < bar * bar;
if !converged || !has_heading(self.estimate.state()) {
return false;
}
let attitude = self.estimate.state().attitude.quaternion();
let y = wrap_pi(yaw - heading_of(attitude));
if !(y.is_finite() && wanted(y)) {
return false;
}
self.reset_heading_by(y, variance, attitude.inverse() * Vector3::z());
true
}
}
#[cfg(test)]
mod tests {
use nalgebra::{UnitQuaternion, Vector3};
use super::super::fixtures::*;
use crate::config::{Config, GRAVITY, Recovery};
use crate::eskf::Eskf;
use crate::gsf::YawEstimator;
use crate::gsf::tests::{circling, legs};
use crate::health::{Fusion, Propagation, Status};
use crate::math::wrap_pi;
use crate::observation::heading::heading_of;
use crate::propagate::ImuSample;
use crate::units::{
Acceleration, AngularRate, HeadingNoise, MagField, Position, PositionNoise, Seconds,
Velocity, VelocityNoise,
};
struct Flight {
yaw: f32,
velocity: Vector3<f32>,
sigma: f32,
age: usize,
last: Option<Fusion>,
}
impl Flight {
fn heading(yaw: f32) -> Self {
Self {
yaw,
velocity: Vector3::zeros(),
sigma: 0.1,
age: 0,
last: None,
}
}
fn fly(
&mut self,
filter: &mut Eskf,
seconds: f32,
acceleration: impl Fn(f32) -> Vector3<f32>,
mut also: impl FnMut(&mut Eskf),
) {
let attitude = UnitQuaternion::from_euler_angles(0.0, 0.0, self.yaw);
let dt = DT.as_secs();
let mut past = std::collections::VecDeque::from([self.velocity]);
for k in 0..(seconds / dt).round() as u32 {
let a = acceleration(k as f32 * dt);
let force = attitude.inverse() * (a - Vector3::z() * GRAVITY);
let imu = ImuSample::reading(AngularRate::zero(), Acceleration::from_vector(force));
assert_eq!(filter.step(imu, DT), Propagation::Propagated);
self.velocity += a * dt;
past.push_back(self.velocity);
if past.len() > self.age + 1 {
past.pop_front();
}
if k % 20 == 19 {
let taken = filter
.now()
.before(Seconds::from_secs(self.age as f32 * dt));
self.last = Some(filter.fuse_gnss_velocity(
taken,
Velocity::from_vector(past[0]),
VelocityNoise::from_speed_accuracy(self.sigma),
Position::zero(),
));
also(filter);
}
}
}
fn error(&self, filter: &Eskf) -> f32 {
wrap_pi(heading_of(filter.state().attitude.quaternion()) - self.yaw).abs()
}
}
fn still(_: f32) -> Vector3<f32> {
Vector3::zeros()
}
fn north(filter: &mut Eskf) {
let _ = filter.fuse_mag_heading(
filter.now(),
MagField::body(0.22, 0.0, 0.44),
HeadingNoise::from_sigma(0.1),
);
}
fn pointing_north(config: Config) -> Eskf {
let mut filter = Eskf::new(config).unwrap();
let _ = filter
.initialize_over(&window_with_mag(), Seconds::from_secs(0.25))
.unwrap();
assert!(filter.state().validity.heading);
filter
}
#[test]
fn a_multirotor_with_no_heading_sensor_takes_its_heading_from_its_legs() {
let mut filter = initialized();
let mut flight = Flight::heading(2.2);
flight.fly(&mut filter, 2.0, still, |_| {});
assert_eq!(filter.state().status, Status::Aligning);
assert_eq!(filter.diagnostics().yaw_estimator.adopted, 0);
flight.fly(&mut filter, 12.0, legs, |_| {});
let yaw = filter.diagnostics().yaw_estimator;
assert_eq!((yaw.adopted, yaw.recovered), (1, 0));
assert!(filter.state().validity.heading);
assert_eq!(filter.state().status, Status::Healthy);
assert!(flight.error(&filter) < 0.05, "{}", flight.error(&filter));
assert!(flight.last.is_some_and(|fusion| fusion.is_accepted()));
}
#[test]
fn a_cruise_with_no_acceleration_establishes_nothing() {
let mut filter = initialized();
let mut flight = Flight::heading(2.2);
flight.fly(&mut filter, 2.0, still, |_| {});
flight.velocity = Vector3::new(2.0, 0.0, 0.0);
let noise = VelocityNoise::from_speed_accuracy(0.1);
assert!(filter.reset_velocity_to(Velocity::from_vector(flight.velocity), noise));
flight.fly(&mut filter, 20.0, still, |_| {});
assert!(filter.yaw_estimator.yaw().is_some());
assert_eq!(filter.diagnostics().yaw_estimator.adopted, 0);
assert_eq!(filter.state().status, Status::Aligning);
}
#[test]
fn with_the_estimator_off_nothing_is_stepped_weighed_or_adopted() {
let mut filter = Eskf::new(Config {
yaw_estimator: false,
..Config::default()
})
.unwrap();
let _ = filter
.initialize_over(&[crate::init::tests::still(); 8], Seconds::from_secs(0.25))
.unwrap();
let mut flight = Flight::heading(2.2);
flight.fly(&mut filter, 2.0, still, |_| {});
flight.fly(&mut filter, 12.0, legs, |_| {});
assert_eq!(filter.yaw_estimator, YawEstimator::default());
assert_eq!(filter.diagnostics().yaw_estimator, Default::default());
assert!(!filter.state().validity.heading);
}
#[test]
fn a_velocity_at_half_a_meter_a_second_of_sigma_is_not_weighed() {
let mut filter = initialized();
let mut flight = Flight::heading(2.2);
flight.sigma = 0.5;
flight.fly(&mut filter, 2.0, still, |_| {});
flight.fly(&mut filter, 12.0, legs, |_| {});
assert_eq!(filter.diagnostics().yaw_estimator.period(), None);
assert_eq!(filter.yaw_estimator.yaw(), None);
}
#[test]
fn a_heading_gnss_contradicts_is_replaced_and_the_velocity_adopted_with_it() {
let mut filter = pointing_north(Config::default());
let mut flight = Flight::heading(1.4);
flight.fly(&mut filter, 2.0, still, north);
assert!(flight.error(&filter) > 1.3);
let mut owed = None;
flight.fly(&mut filter, 12.0, legs, |filter| {
north(filter);
if filter.diagnostics().yaw_estimator.recovered == 1 && owed.is_none() {
owed = Some(filter.yaw_replaced.owes_velocity());
}
});
assert_eq!(
owed,
Some(false),
"the adoption settles what the replacement owed"
);
let d = filter.diagnostics();
assert_eq!((d.yaw_estimator.adopted, d.yaw_estimator.recovered), (1, 1));
assert!(
d.gnss_velocity.rejected >= 4,
"{}",
d.gnss_velocity.rejected
);
assert_eq!(d.gnss_velocity.recovered, 1);
assert!(flight.error(&filter) < 0.05, "{}", flight.error(&filter));
assert!(d.mag_heading.consecutive_rejections > 10);
assert!(flight.last.is_some_and(|fusion| fusion.is_accepted()));
}
#[test]
fn a_receiver_fault_under_a_good_heading_leaves_the_yaw_alone() {
let mut filter = pointing_north(Config::default());
let mut flight = Flight::heading(0.0);
flight.fly(&mut filter, 2.0, still, north);
flight.fly(&mut filter, 12.0, legs, north);
let noise = VelocityNoise::from_speed_accuracy(0.6);
let wild = Velocity::from_vector(flight.velocity + Vector3::new(30.0, 0.0, 0.0));
hold(&mut filter, 3.0, 20, |filter| {
let fusion = filter.fuse_gnss_velocity(filter.now(), wild, noise, Position::zero());
assert!(matches!(fusion, Fusion::Rejected { .. }), "{fusion:?}");
});
assert_eq!(filter.diagnostics().yaw_estimator.adopted, 0);
assert!(filter.yaw_estimator.yaw().is_some());
}
#[test]
fn after_the_yaw_is_replaced_the_first_fix_the_gate_turns_down_is_adopted() {
let mut filter = pointing_north(Config {
recovery: Recovery {
gnss_position: Some(Seconds::from_secs(1000.0)),
..Recovery::default()
},
..Config::default()
});
let mut flight = Flight::heading(1.4);
flight.fly(&mut filter, 2.0, still, north);
flight.fly(&mut filter, 12.0, legs, north);
assert_eq!(filter.diagnostics().yaw_estimator.recovered, 1);
let far =
Position::from_vector(filter.state().position.vector() + Vector3::new(500.0, 0.0, 0.0));
let fix = |filter: &mut Eskf| {
filter.fuse_gnss_position(filter.now(), far, one_metre(), Position::zero())
};
assert_eq!(fix(&mut filter).horizontal, Fusion::Reset);
let _ = filter.reset_position_to(Position::zero(), one_metre());
assert!(matches!(
fix(&mut filter).horizontal,
Fusion::Rejected { .. }
));
let mut filter = pointing_north(Config {
recovery: Recovery {
gnss_position: Some(Seconds::from_secs(1000.0)),
..Recovery::default()
},
..Config::default()
});
let mut flight = Flight::heading(1.4);
flight.fly(&mut filter, 2.0, still, north);
flight.fly(&mut filter, 12.0, legs, north);
let here = filter.state().position;
let wide = PositionNoise::from_sigma(50.0, 50.0, 50.0);
let passed = filter.fuse_gnss_position(filter.now(), here, wide, Position::zero());
assert!(passed.horizontal.is_accepted(), "{passed:?}");
let far = Position::from_vector(here.vector() + Vector3::new(500.0, 0.0, 0.0));
let turned_down =
filter.fuse_gnss_position(filter.now(), far, one_metre(), Position::zero());
assert!(matches!(turned_down.horizontal, Fusion::Rejected { .. }));
let mut filter = pointing_north(Config::default());
let mut flight = Flight::heading(1.4);
flight.fly(&mut filter, 2.0, still, north);
flight.fly(&mut filter, 12.0, legs, north);
assert!(filter.yaw_replaced.owes_position());
let _ = filter.initialize_over(&window_with_mag(), Seconds::from_secs(0.25));
assert_eq!(filter.yaw_replaced, super::YawReplaced::default());
}
#[test]
fn an_estimator_begun_again_in_flight_does_not_overrule_until_it_settles() {
let mut filter = pointing_north(Config::default());
let mut flight = Flight::heading(1.4);
flight.fly(&mut filter, 2.0, still, north);
let gap = filter.step(crate::init::tests::still().imu, Seconds::from_secs(0.5));
assert!(matches!(gap, Propagation::Coasted { .. }));
flight.fly(&mut filter, 8.0, legs, north);
assert!(filter.yaw_estimator.yaw().is_some_and(|(_, v)| v < 0.01));
assert!(filter.diagnostics().gnss_velocity.rejected > 5);
assert_eq!(filter.diagnostics().yaw_estimator.recovered, 0);
}
#[test]
fn an_old_velocity_reaches_the_estimator_carried_to_now() {
let mut filter = initialized();
let mut flight = Flight::heading(0.9);
flight.age = 20;
flight.fly(&mut filter, 2.0, still, |_| {});
flight.fly(&mut filter, 15.0, circling, |_| {});
let (yaw, _) = filter.yaw_estimator.yaw().unwrap();
let error = wrap_pi(yaw - 0.9);
assert!(error.abs() < 0.03, "{error}");
}
#[test]
fn a_second_antenna_that_is_being_accepted_is_not_overruled() {
let mut filter = pointing_north(Config::default());
let mut flight = Flight::heading(1.4);
let both = |filter: &mut Eskf| {
north(filter);
let zero = crate::units::Radians::from_radians(0.0);
let _ = filter.fuse_gnss_heading(filter.now(), zero, HeadingNoise::from_sigma(0.02));
};
flight.fly(&mut filter, 2.0, still, both);
flight.fly(&mut filter, 12.0, legs, both);
assert!(filter.diagnostics().gnss_velocity.rejected > 5);
assert_eq!(filter.diagnostics().yaw_estimator.adopted, 0);
}
#[test]
fn with_its_recovery_off_a_wrong_heading_is_reported_and_kept() {
let mut filter = pointing_north(Config {
recovery: Recovery {
yaw_estimator: None,
..Recovery::default()
},
..Config::default()
});
let mut flight = Flight::heading(1.4);
flight.fly(&mut filter, 2.0, still, north);
flight.fly(&mut filter, 6.0, legs, north);
assert_eq!(filter.diagnostics().yaw_estimator.adopted, 0);
assert!(filter.diagnostics().gnss_velocity.rejected > 0);
assert!(flight.error(&filter) > 1.0);
}
#[test]
fn a_coasted_gap_starts_the_estimator_over_and_it_must_settle_again() {
let mut filter = pointing_north(Config::default());
let mut flight = Flight::heading(0.0);
flight.fly(&mut filter, 2.0, still, north);
flight.fly(&mut filter, 6.0, legs, north);
assert!(filter.yaw_estimator.yaw().is_some());
let gap = filter.step(crate::init::tests::still().imu, Seconds::from_secs(0.5));
assert!(matches!(gap, Propagation::Coasted { .. }));
assert_eq!(filter.yaw_estimator.yaw(), None);
assert!(!filter.yaw_estimator.is_settled());
}
#[test]
fn a_step_refused_after_the_clock_moved_starts_the_estimator_over() {
let mut filter = pointing_north(Config {
coast: None,
..Config::default()
});
let mut flight = Flight::heading(0.0);
flight.fly(&mut filter, 2.0, still, north);
flight.fly(&mut filter, 6.0, legs, north);
assert!(filter.yaw_estimator.yaw().is_some());
let gap = filter.step(crate::init::tests::still().imu, Seconds::from_secs(5.0));
assert!(matches!(gap, Propagation::StepTooLong { .. }));
assert_eq!(filter.yaw_estimator.yaw(), None);
assert!(!filter.yaw_estimator.is_settled());
flight.fly(&mut filter, 4.0, legs, north);
let held = filter.yaw_estimator.clone();
let stale = crate::init::tests::still().imu.timed(filter.now(), DT);
assert!(matches!(
filter.predict(stale),
Propagation::InvalidStep { .. }
));
assert_eq!(filter.yaw_estimator, held);
}
#[test]
fn a_first_heading_waits_for_an_estimator_begun_in_flight_to_settle() {
let mut filter = initialized();
let mut flight = Flight::heading(2.2);
flight.fly(&mut filter, 2.0, still, |_| {});
let gap = filter.step(crate::init::tests::still().imu, Seconds::from_secs(0.5));
assert!(matches!(gap, Propagation::Coasted { .. }));
flight.fly(&mut filter, 8.0, legs, |_| {});
assert!(filter.yaw_estimator.yaw().is_some_and(|(_, v)| v < 0.01));
assert_eq!(filter.diagnostics().yaw_estimator.adopted, 0);
flight.fly(&mut filter, 8.0, legs, |_| {});
assert_eq!(filter.diagnostics().yaw_estimator.adopted, 1);
}
#[test]
fn one_velocity_turned_down_is_not_a_yaw_failure() {
let mut filter = pointing_north(Config::default());
let mut flight = Flight::heading(0.0);
flight.fly(&mut filter, 2.0, still, north);
flight.fly(&mut filter, 12.0, legs, north);
let down = filter.state().attitude.quaternion().inverse() * Vector3::z();
filter.reset_heading_by(1.4, 0.01, down);
hold(&mut filter, 1.5, 1000, |_| {});
let noise = VelocityNoise::from_speed_accuracy(0.6);
let wild = Velocity::from_vector(flight.velocity + Vector3::new(30.0, 0.0, 0.0));
let first = filter.fuse_gnss_velocity(filter.now(), wild, noise, Position::zero());
assert!(matches!(first, Fusion::Rejected { .. }), "{first:?}");
assert_eq!(filter.diagnostics().yaw_estimator.recovered, 0);
let second = filter.fuse_gnss_velocity(filter.now(), wild, noise, Position::zero());
assert_eq!(second, Fusion::Reset);
assert_eq!(filter.diagnostics().yaw_estimator.recovered, 1);
}
#[test]
fn an_overruled_magnetometer_is_not_adopted_back_when_gnss_goes() {
let mut filter = pointing_north(Config::default());
let mut flight = Flight::heading(1.4);
flight.fly(&mut filter, 2.0, still, north);
flight.fly(&mut filter, 12.0, legs, north);
assert_eq!(filter.diagnostics().yaw_estimator.recovered, 1);
hold(&mut filter, 12.0, 20, north);
assert_eq!(filter.diagnostics().mag_heading.recovered, 0);
assert!(flight.error(&filter) < 0.3, "{}", flight.error(&filter));
}
#[test]
fn a_magnetometer_that_agrees_again_is_trusted_again() {
let mut filter = pointing_north(Config::default());
let mut flight = Flight::heading(1.4);
flight.fly(&mut filter, 2.0, still, north);
flight.fly(&mut filter, 12.0, legs, north);
assert!(filter.yaw_replaced.overrules_magnetometer());
let (sin, cos) = 1.4_f32.sin_cos();
let field = MagField::body(0.22 * cos, -0.22 * sin, 0.44);
let fusion = filter.fuse_mag_heading(filter.now(), field, HeadingNoise::from_sigma(0.1));
assert!(fusion.is_accepted(), "{fusion:?}");
assert!(!filter.yaw_replaced.overrules_magnetometer());
}
#[test]
fn a_velocity_before_the_estimator_has_leveled_is_not_an_arrival() {
let mut filter = initialized();
let mut flight = Flight::heading(2.2);
flight.fly(&mut filter, 0.5, still, |_| {});
assert_eq!(filter.diagnostics().yaw_estimator, Default::default());
flight.fly(&mut filter, 5.0, still, |_| {});
assert!(filter.diagnostics().yaw_estimator.period().is_some());
}
#[test]
fn a_new_start_forgets_the_estimator() {
let mut filter = initialized();
let mut flight = Flight::heading(2.2);
flight.fly(&mut filter, 2.0, still, |_| {});
flight.fly(&mut filter, 12.0, legs, |_| {});
let _ = filter
.initialize_over(&[crate::init::tests::still(); 8], Seconds::from_secs(0.25))
.unwrap();
assert_eq!(filter.yaw_estimator, YawEstimator::default());
assert_eq!(filter.yaw_replaced, super::YawReplaced::default());
}
}