use serde::Serialize;
pub fn combine_power_dbm(powers_dbm: &[f64]) -> f64 {
if powers_dbm.is_empty() {
return f64::NEG_INFINITY;
}
let linear: f64 = powers_dbm.iter().map(|p| 10.0_f64.powf(p / 10.0)).sum();
10.0 * linear.log10()
}
#[derive(Clone, Copy, Debug)]
pub struct AgcMonitor {
pub expected_dbm: f64,
pub alert_margin_db: f64,
}
impl AgcMonitor {
pub fn new(expected_dbm: f64) -> Self {
Self {
expected_dbm,
alert_margin_db: 3.0,
}
}
pub fn excess_db(&self, measured_dbm: f64) -> f64 {
measured_dbm - self.expected_dbm
}
pub fn alert(&self, measured_dbm: f64) -> bool {
self.excess_db(measured_dbm) > self.alert_margin_db
}
}
pub fn bpsk_autocorr(tau_chips: f64) -> f64 {
(1.0 - tau_chips.abs()).max(0.0)
}
pub fn early_late_ideal(spacing_chips: f64) -> (f64, f64) {
let r = bpsk_autocorr(spacing_chips);
(r, r)
}
#[derive(Clone, Copy, Debug)]
pub struct SqmMonitor {
pub el_tolerance: f64,
}
impl SqmMonitor {
pub fn new() -> Self {
Self { el_tolerance: 0.10 }
}
pub fn el_metric(&self, early: f64, late: f64) -> f64 {
let sum = early + late;
if sum.abs() < 1e-300 {
0.0
} else {
(early - late) / sum
}
}
pub fn alert(&self, early: f64, late: f64) -> bool {
self.el_metric(early, late).abs() > self.el_tolerance
}
}
impl Default for SqmMonitor {
fn default() -> Self {
Self::new()
}
}
use crate::detection::chi2_inv_cdf;
use crate::fusion::ukf::inverse;
#[derive(Clone, Copy, Debug, Serialize)]
pub struct RaimConsistency {
pub statistic: f64,
pub threshold: f64,
pub dof: usize,
pub alert: bool,
}
pub fn parity_raim_test(
geometry: &[[f64; 4]],
residuals: &[f64],
sigma: f64,
p_fa: f64,
) -> Option<RaimConsistency> {
let m = geometry.len();
if m < 5 || residuals.len() != m || sigma <= 0.0 {
return None;
}
let n = 4;
let mut nmat = vec![vec![0.0; n]; n];
let mut b = vec![0.0; n];
for (g, &z) in geometry.iter().zip(residuals) {
for a in 0..n {
b[a] += g[a] * z;
for (c, ncell) in nmat[a].iter_mut().enumerate() {
*ncell += g[a] * g[c];
}
}
}
let ninv = inverse(&nmat)?;
let x: Vec<f64> = (0..n)
.map(|a| (0..n).map(|c| ninv[a][c] * b[c]).sum())
.collect();
let t: f64 = geometry
.iter()
.zip(residuals)
.map(|(g, &z)| {
let pred: f64 = (0..n).map(|a| g[a] * x[a]).sum();
let r = z - pred;
r * r
})
.sum::<f64>()
/ (sigma * sigma);
let dof = m - n;
let threshold = chi2_inv_cdf(1.0 - p_fa, dof as f64);
Some(RaimConsistency {
statistic: t,
threshold,
dof,
alert: t > threshold,
})
}
#[derive(Clone, Copy, Debug, Default, PartialEq, Eq, Serialize)]
pub struct SpoofLayers {
pub raim: bool,
pub agc: bool,
pub sqm: bool,
}
impl SpoofLayers {
pub fn count(self) -> usize {
self.raim as usize + self.agc as usize + self.sqm as usize
}
}
#[derive(Clone, Copy, Debug, Serialize)]
pub struct FusedSpoofDecision {
pub layers: SpoofLayers,
pub score: f64,
pub alert: bool,
}
pub fn fuse_spoof_layers(
raim: bool,
agc: bool,
sqm: bool,
weights: [f64; 3],
threshold: f64,
) -> FusedSpoofDecision {
let layers = SpoofLayers { raim, agc, sqm };
let score = raim as u8 as f64 * weights[0]
+ agc as u8 as f64 * weights[1]
+ sqm as u8 as f64 * weights[2];
FusedSpoofDecision {
layers,
score,
alert: score >= threshold,
}
}
#[derive(Clone, Debug)]
pub struct SpoofEpoch {
pub geometry: Vec<[f64; 4]>,
pub residuals: Vec<f64>,
pub sigma_m: f64,
pub measured_dbm: f64,
pub early: f64,
pub late: f64,
}
#[derive(Clone, Copy, Debug)]
pub struct CombinedSpoofDetector {
pub agc: AgcMonitor,
pub sqm: SqmMonitor,
pub raim_p_fa: f64,
pub weights: [f64; 3],
pub fusion_threshold: f64,
}
impl CombinedSpoofDetector {
pub fn new(expected_dbm: f64) -> Self {
Self {
agc: AgcMonitor::new(expected_dbm),
sqm: SqmMonitor::new(),
raim_p_fa: 1.0e-3,
weights: [0.5, 0.3, 0.2],
fusion_threshold: 0.5,
}
}
pub fn evaluate(&self, epoch: &SpoofEpoch) -> CombinedSpoofDecision {
let raim = parity_raim_test(
&epoch.geometry,
&epoch.residuals,
epoch.sigma_m,
self.raim_p_fa,
);
let raim_alert = raim.map(|r| r.alert).unwrap_or(false);
let agc_excess_db = self.agc.excess_db(epoch.measured_dbm);
let agc_alert = self.agc.alert(epoch.measured_dbm);
let sqm_el_metric = self.sqm.el_metric(epoch.early, epoch.late);
let sqm_alert = self.sqm.alert(epoch.early, epoch.late);
let fused = fuse_spoof_layers(
raim_alert,
agc_alert,
sqm_alert,
self.weights,
self.fusion_threshold,
);
CombinedSpoofDecision {
raim,
agc_excess_db,
sqm_el_metric,
fused,
}
}
}
#[derive(Clone, Copy, Debug, Serialize)]
pub struct CombinedSpoofDecision {
pub raim: Option<RaimConsistency>,
pub agc_excess_db: f64,
pub sqm_el_metric: f64,
pub fused: FusedSpoofDecision,
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn power_combines_incoherently() {
let two = combine_power_dbm(&[-130.0, -130.0]);
assert!((two - (-126.9897)).abs() < 1e-3, "two = {two}");
let eight = combine_power_dbm(&[-130.0; 8]);
assert!(
(eight - (-130.0 + 10.0 * 8.0_f64.log10())).abs() < 1e-9,
"eight = {eight}"
);
assert_eq!(combine_power_dbm(&[]), f64::NEG_INFINITY);
}
#[test]
fn agc_flags_only_a_real_power_excess() {
let nominal: Vec<f64> = vec![-130.0; 8];
let floor = combine_power_dbm(&nominal);
let mon = AgcMonitor::new(floor);
assert!((mon.excess_db(floor)).abs() < 1e-9);
assert!(!mon.alert(floor));
assert!(!mon.alert(floor + 2.0));
assert!(mon.alert(floor + 4.0));
let mut spoofed = nominal.clone();
spoofed.push(-115.0);
let measured = combine_power_dbm(&spoofed);
assert!(mon.alert(measured), "excess {} dB", mon.excess_db(measured));
}
#[test]
fn bpsk_autocorr_is_the_triangle() {
assert!((bpsk_autocorr(0.0) - 1.0).abs() < 1e-12);
assert!((bpsk_autocorr(0.1) - 0.9).abs() < 1e-12);
assert!((bpsk_autocorr(-0.5) - 0.5).abs() < 1e-12);
assert_eq!(bpsk_autocorr(1.0), 0.0);
assert_eq!(bpsk_autocorr(2.0), 0.0);
}
#[test]
fn sqm_flags_only_real_correlation_asymmetry() {
let mon = SqmMonitor::new();
let (e, l) = early_late_ideal(0.1);
assert!((mon.el_metric(e, l)).abs() < 1e-12);
assert!(!mon.alert(e, l));
assert!(!mon.alert(0.95, 0.85));
assert!(mon.alert(1.0, 0.8));
assert_eq!(mon.el_metric(0.0, 0.0), 0.0);
}
fn geometry8() -> Vec<[f64; 4]> {
let azels = [
(0.0, 80.0),
(45.0, 30.0),
(100.0, 55.0),
(150.0, 20.0),
(200.0, 60.0),
(255.0, 25.0),
(300.0, 45.0),
(340.0, 15.0),
];
azels
.iter()
.map(|&(az, el): &(f64, f64)| {
let (a, e) = (az.to_radians(), el.to_radians());
[e.cos() * a.sin(), e.cos() * a.cos(), e.sin(), 1.0]
})
.collect()
}
#[test]
fn parity_raim_flags_an_inconsistent_satellite() {
let g = geometry8();
let x_true = [12.0, -8.0, 5.0, 30.0]; let consistent: Vec<f64> = g
.iter()
.map(|row| (0..4).map(|a| row[a] * x_true[a]).sum())
.collect();
let clean = parity_raim_test(&g, &consistent, 5.0, 0.001).expect("redundant");
assert!(clean.statistic < 1e-9 && !clean.alert, "{clean:?}");
let mut spoofed = consistent.clone();
spoofed[3] += 60.0;
let attacked = parity_raim_test(&g, &spoofed, 5.0, 0.001).expect("redundant");
assert!(attacked.alert, "single-SV bias not flagged: {attacked:?}");
}
#[test]
fn common_mode_bias_is_raim_invisible() {
let g = geometry8();
let base: Vec<f64> = (0..8).map(|i| (i as f64 - 3.5) * 1.3).collect();
let s0 = parity_raim_test(&g, &base, 5.0, 0.001).expect("redundant");
let ramped: Vec<f64> = base.iter().map(|&z| z + 250.0).collect();
let s1 = parity_raim_test(&g, &ramped, 5.0, 0.001).expect("redundant");
assert!(
(s0.statistic - s1.statistic).abs() < 1e-6,
"common-mode bias changed the statistic: {} vs {}",
s0.statistic,
s1.statistic
);
}
#[test]
fn fusion_combines_independent_layers() {
let w = [0.5, 0.3, 0.2];
let lone_sqm = fuse_spoof_layers(false, false, true, w, 0.5);
assert!(!lone_sqm.alert && lone_sqm.layers.count() == 1);
let raim_agc = fuse_spoof_layers(true, true, false, w, 0.5);
assert!(raim_agc.alert && raim_agc.layers.count() == 2);
assert!(raim_agc.layers.raim && raim_agc.layers.agc && !raim_agc.layers.sqm);
}
#[test]
fn parity_raim_false_alert_and_missed_detection_are_characterized() {
use rand::SeedableRng;
use rand_chacha::ChaCha8Rng;
use rand_distr::{Distribution, Normal};
let g = geometry8();
let sigma = 5.0;
let p_fa = 0.05;
let noise = Normal::new(0.0, sigma).unwrap();
let trials = 400;
let mut rng = ChaCha8Rng::seed_from_u64(0x5F00_F123);
let run = |bias: f64, rng: &mut ChaCha8Rng| {
let mut alerts = 0;
for _ in 0..trials {
let mut z: Vec<f64> = (0..8).map(|_| noise.sample(rng)).collect();
z[3] += bias;
if parity_raim_test(&g, &z, sigma, p_fa).unwrap().alert {
alerts += 1;
}
}
alerts as f64 / trials as f64
};
let pfa = run(0.0, &mut rng);
assert!((0.02..0.10).contains(&pfa), "empirical P_fa = {pfa}");
let pmd_small = 1.0 - run(2.0 * sigma, &mut rng);
let pmd_large = 1.0 - run(8.0 * sigma, &mut rng);
assert!(
pmd_large < pmd_small,
"P_md did not fall with bias: {pmd_large} vs {pmd_small}"
);
assert!(pmd_large < 0.2, "P_md at 8σ still high: {pmd_large}");
}
fn consistent_residuals(g: &[[f64; 4]]) -> Vec<f64> {
let x_true = [9.0, -4.0, 6.0, 25.0];
g.iter()
.map(|row| (0..4).map(|a| row[a] * x_true[a]).sum())
.collect()
}
#[test]
fn combined_detector_passes_a_clean_epoch() {
let g = geometry8();
let residuals = consistent_residuals(&g);
let floor = combine_power_dbm(&[-130.0; 8]);
let det = CombinedSpoofDetector::new(floor);
let (early, late) = early_late_ideal(0.1);
let epoch = SpoofEpoch {
geometry: g,
residuals,
sigma_m: 5.0,
measured_dbm: floor,
early,
late,
};
let d = det.evaluate(&epoch);
assert!(!d.fused.alert, "clean epoch raised a spoof alert: {d:?}");
assert_eq!(d.fused.layers.count(), 0, "no layer should fire: {d:?}");
}
#[test]
fn combined_detector_catches_a_single_sv_bias_via_raim() {
let g = geometry8();
let mut residuals = consistent_residuals(&g);
residuals[3] += 60.0; let floor = combine_power_dbm(&[-130.0; 8]);
let det = CombinedSpoofDetector::new(floor);
let (early, late) = early_late_ideal(0.1);
let epoch = SpoofEpoch {
geometry: g,
residuals,
sigma_m: 5.0,
measured_dbm: floor,
early,
late,
};
let d = det.evaluate(&epoch);
assert!(
d.raim.map(|r| r.alert).unwrap_or(false),
"RAIM did not fire on a single-SV bias: {d:?}"
);
assert!(
d.fused.layers.raim && !d.fused.layers.agc && !d.fused.layers.sqm,
"only RAIM should fire: {d:?}"
);
assert!(
d.fused.alert,
"single-SV bias not flagged by the fused detector"
);
}
#[test]
fn combined_detector_catches_a_raim_invisible_meaconer_via_agc_and_sqm() {
let g = geometry8();
let base: Vec<f64> = (0..8).map(|i| (i as f64 - 3.5) * 1.1).collect();
let residuals: Vec<f64> = base.iter().map(|&z| z + 300.0).collect(); let floor = combine_power_dbm(&[-130.0; 8]);
let det = CombinedSpoofDetector::new(floor);
let epoch = SpoofEpoch {
geometry: g,
residuals,
sigma_m: 5.0,
measured_dbm: floor + 6.0, early: 1.0,
late: 0.72, };
let d = det.evaluate(&epoch);
assert!(
!d.raim.map(|r| r.alert).unwrap_or(false),
"RAIM must be blind to a common-mode bias: {d:?}"
);
assert!(
d.fused.layers.agc && d.fused.layers.sqm,
"AGC + SQM should both fire on the over-powered meaconer: {d:?}"
);
assert!(
d.fused.alert,
"RAIM-invisible meaconer escaped the fused detector: {d:?}"
);
}
#[test]
fn combined_detector_reports_per_layer_diagnostics() {
let g = geometry8();
let residuals = consistent_residuals(&g);
let floor = combine_power_dbm(&[-130.0; 8]);
let det = CombinedSpoofDetector::new(floor);
let epoch = SpoofEpoch {
geometry: g,
residuals,
sigma_m: 5.0,
measured_dbm: floor + 1.5,
early: 1.0,
late: 0.95,
};
let d = det.evaluate(&epoch);
assert!(
d.raim.is_some(),
"redundant geometry should yield a RAIM statistic"
);
assert!(
(d.agc_excess_db - 1.5).abs() < 1e-9,
"AGC excess should be +1.5 dB, got {}",
d.agc_excess_db
);
assert!(
d.sqm_el_metric > 0.0,
"Early > Late should give a positive imbalance, got {}",
d.sqm_el_metric
);
}
}