stabilizer_stream/
var.rs

1use core::f32::consts::PI;
2use derive_builder::Builder;
3
4#[derive(Debug, Builder, Clone, Copy, PartialEq)]
5pub struct Var {
6    /// exponent of `pi*f*tau` in the variance frequency response (-2 for AVAR, -4 for MVAR)
7    #[builder(default = "-2")]
8    x_exp: i32,
9    /// Exponent of `sin(pi*f*tau)` in the variance frequency response (4 for AVAR, 6 for MVAR)
10    #[builder(default = "4")]
11    sinx_exp: i32,
12    /// Response clip (infinite for AVAR and MVAR, 1 for the main lobe: FVAR)
13    #[builder(default = "f32::MAX")]
14    clip: f32,
15    /// skip the first `dc_cut` bins to suppress DC window leakage
16    #[builder(default = "2")]
17    dc_cut: usize,
18}
19
20impl Var {
21    /// Compute statistical variance estimator (AVAR, MVAR, FVAR...) from Phase PSD
22    ///
23    /// # Args
24    /// * `phase_psd`: Phase noise PSD vector from Nyquist down
25    /// * `frequencies`: PSD bin frequencies, Nyquist first
26    pub fn eval(&self, phase_psd: &[f32], frequencies: &[f32], tau: f32) -> f32 {
27        phase_psd
28            .iter()
29            .rev()
30            .zip(
31                frequencies
32                    .iter()
33                    .rev()
34                    .take_while(|&f| f * tau <= self.clip),
35            )
36            .skip(self.dc_cut)
37            // force DC bin to 0
38            .fold((0.0, (0.0, 0.0)), |(accu, (a0, f0)), (&sp, &f)| {
39                // frequency PSD
40                let sy = sp * f * f;
41                let pft = PI * (f * tau);
42                // Allan variance transfer function (rectangular window: sinc**2 and differencing: 2*sin**2)
43                // Cancel the 2 here with the 0.5 in the trapezoidal rule
44                let hahd = pft.sin().powi(self.sinx_exp) * pft.powi(self.x_exp);
45                let a = sy * hahd;
46                // trapezoidal integration
47                (accu + (a + a0) * (f - f0), (a, f))
48            })
49            .0
50    }
51}
52
53#[cfg(test)]
54mod test {
55    use super::*;
56
57    #[test]
58    fn basic() {
59        let mut p = [1000.0, 100.0, 1.2, 3.4, 5.6];
60        let mut f = [0.0, 1.0, 3.0, 6.0, 9.0];
61        p.reverse();
62        f.reverse();
63        let var = VarBuilder::default().build().unwrap();
64        let v = var.eval(&p, &f, 2.7);
65        println!("{}", v);
66        assert!((0.13478442 - v).abs() < 1e-6);
67    }
68}
stabilizer_stream/var.rs

stabilizer_stream/
var.rs
