nyx_space/od/noise/
white.rs

1/*
2    Nyx, blazing fast astrodynamics
3    Copyright (C) 2018-onwards Christopher Rabotin <christopher.rabotin@gmail.com>
4
5    This program is free software: you can redistribute it and/or modify
6    it under the terms of the GNU Affero General Public License as published
7    by the Free Software Foundation, either version 3 of the License, or
8    (at your option) any later version.
9
10    This program is distributed in the hope that it will be useful,
11    but WITHOUT ANY WARRANTY; without even the implied warranty of
12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13    GNU Affero General Public License for more details.
14
15    You should have received a copy of the GNU Affero General Public License
16    along with this program.  If not, see <https://www.gnu.org/licenses/>.
17*/
18
19use std::ops::{Mul, MulAssign};
20
21use anise::constants::SPEED_OF_LIGHT_KM_S;
22use hifitime::{Duration, Epoch};
23use rand::Rng;
24use rand_distr::Normal;
25use serde_derive::{Deserialize, Serialize};
26
27use super::Stochastics;
28
29/// White noise is an uncorrelated random variable.
30#[derive(Copy, Clone, Debug, Default, PartialEq, Serialize, Deserialize)]
31pub struct WhiteNoise {
32    /// Mean value of this white noise
33    pub mean: f64,
34    /// Process noise as a one-sigma of the Normal distribution.
35    pub sigma: f64,
36}
37
38impl WhiteNoise {
39    /// Initializes a new random walk stochastic noise model from the process noise and the integration time.
40    /// This will compute the process noise per second automatically.
41    pub fn new(process_noise: f64, integration_time: Duration) -> Self {
42        Self {
43            sigma: process_noise / integration_time.to_seconds(),
44            ..Default::default()
45        }
46    }
47
48    /// Initializes a new random walk stochastic noise model from the provided process noise, assuming that the noise level
49    /// is fixed regardless of the integration time.
50    pub fn constant_white_noise(process_noise: f64) -> Self {
51        Self {
52            sigma: process_noise,
53            ..Default::default()
54        }
55    }
56
57    // Initializes a new time-uncorrelated white noise process, using only the Pr/N0 value and the bandwidth.
58    // This returns a white noise sigma in kilometers.
59    //
60    // # Equation
61    // σ = c / (2 * B * √(Pr/N0))
62    //
63    // Where c is the speed of light, B is the bandwidth in Hz, and the Pr/N0 is the signal-to-noise ratio.
64    pub fn from_pr_n0(pr_n0: f64, bandwidth_hz: f64) -> Self {
65        Self {
66            sigma: SPEED_OF_LIGHT_KM_S / (2.0 * bandwidth_hz * (pr_n0).sqrt()),
67            mean: 0.0,
68        }
69    }
70}
71
72impl Stochastics for WhiteNoise {
73    fn covariance(&self, _epoch: Epoch) -> f64 {
74        self.sigma.powi(2)
75    }
76
77    fn sample<R: Rng>(&mut self, _epoch: Epoch, rng: &mut R) -> f64 {
78        rng.sample(Normal::new(self.mean, self.sigma).unwrap())
79    }
80}
81
82impl Mul<f64> for WhiteNoise {
83    type Output = Self;
84
85    /// Scale the white noise sigmas by a constant.
86    fn mul(mut self, rhs: f64) -> Self::Output {
87        self.sigma *= rhs;
88        self
89    }
90}
91
92impl MulAssign<f64> for WhiteNoise {
93    fn mul_assign(&mut self, rhs: f64) {
94        *self = *self * rhs;
95    }
96}
97
98#[cfg(test)]
99mod ut_wn {
100    use hifitime::{Epoch, TimeUnits};
101    use rand_pcg::Pcg64Mcg;
102
103    use super::{Stochastics, WhiteNoise};
104
105    #[test]
106    fn white_noise_test() {
107        let sigma = 10.0_f64;
108        let mut wn = WhiteNoise { mean: 0.0, sigma };
109
110        let mut larger_wn = WhiteNoise {
111            mean: 0.0,
112            sigma: sigma * 10.0,
113        };
114
115        let epoch = Epoch::now().unwrap();
116
117        let mut rng = Pcg64Mcg::new(1000);
118        let mut cnt_above_3sigma = 0;
119        let mut cnt_below_3sigma = 0;
120        let mut larger_cnt_above_3sigma = 0;
121        let mut larger_cnt_below_3sigma = 0;
122        for seconds in 0..1000_i64 {
123            let bias = wn.sample(epoch + seconds.seconds(), &mut rng);
124
125            if bias > 3.0 * sigma {
126                cnt_above_3sigma += 1;
127            } else if bias < -3.0 * sigma {
128                cnt_below_3sigma += 1;
129            }
130
131            let larger_bias = larger_wn.sample(epoch + seconds.seconds(), &mut rng);
132            if larger_bias > 30.0 * sigma {
133                larger_cnt_above_3sigma += 1;
134            } else if larger_bias < -30.0 * sigma {
135                larger_cnt_below_3sigma += 1;
136            }
137        }
138
139        assert!(dbg!(cnt_above_3sigma) <= 3);
140        assert!(dbg!(cnt_below_3sigma) <= 3);
141
142        assert!(dbg!(larger_cnt_above_3sigma) <= 3);
143        assert!(dbg!(larger_cnt_below_3sigma) <= 3);
144    }
145}
nyx_space/od/noise/white.rs

nyx_space/od/noise/
white.rs
