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 crate::io::ConfigError;
22use anise::constants::SPEED_OF_LIGHT_KM_S;
23use der::Sequence;
24use hifitime::{Duration, Epoch};
25use rand::{Rng, RngExt};
26use rand_distr::Normal;
27use serde_derive::{Deserialize, Serialize};
28
29#[cfg(feature = "python")]
30use pyo3::prelude::*;
31
32use super::Stochastics;
33
34/// White noise is an uncorrelated random variable.
35#[derive(Copy, Clone, Debug, Default, PartialEq, Serialize, Deserialize, Sequence)]
36#[cfg_attr(feature = "python", pyclass(get_all, set_all))]
37pub struct WhiteNoise {
38    /// Mean value of this white noise
39    pub mean: f64,
40    /// Process noise as a one-sigma of the Normal distribution.
41    pub sigma: f64,
42}
43
44impl WhiteNoise {
45    /// Initializes a new random walk stochastic noise model from the process noise and the integration time.
46    /// This will compute the process noise per second automatically.
47    pub fn new(process_noise: f64, integration_time: Duration) -> Result<Self, ConfigError> {
48        if process_noise.is_sign_negative() {
49            return Err(ConfigError::InvalidConfig {
50                msg: format!("process noise must be positive: {process_noise}"),
51            });
52        }
53        if integration_time.to_seconds() <= 0.0 {
54            return Err(ConfigError::InvalidConfig {
55                msg: format!("integration time must be positive: {integration_time}"),
56            });
57        }
58        Ok(Self {
59            sigma: process_noise / integration_time.to_seconds(),
60            ..Default::default()
61        })
62    }
63
64    /// Initializes a new random walk stochastic noise model from the provided process noise, assuming that the noise level
65    /// is fixed regardless of the integration time.
66    pub fn constant_white_noise(process_noise: f64) -> Self {
67        Self {
68            sigma: process_noise,
69            ..Default::default()
70        }
71    }
72
73    // Initializes a new time-uncorrelated white noise process, using only the Pr/N0 value and the bandwidth.
74    // This returns a white noise sigma in kilometers.
75    //
76    // # Equation
77    // σ = c / (2 * B * √(Pr/N0))
78    //
79    // Where c is the speed of light, B is the bandwidth in Hz, and the Pr/N0 is the signal-to-noise ratio.
80    pub fn from_pr_n0(pr_n0: f64, bandwidth_hz: f64) -> Self {
81        Self {
82            sigma: SPEED_OF_LIGHT_KM_S / (2.0 * bandwidth_hz * (pr_n0).sqrt()),
83            mean: 0.0,
84        }
85    }
86}
87
88impl Stochastics for WhiteNoise {
89    fn covariance(&self, _epoch: Epoch) -> f64 {
90        self.sigma.powi(2)
91    }
92
93    fn sample<R: Rng>(&mut self, _epoch: Epoch, rng: &mut R) -> f64 {
94        rng.sample(Normal::new(self.mean, self.sigma).unwrap())
95    }
96}
97
98impl Mul<f64> for WhiteNoise {
99    type Output = Self;
100
101    /// Scale the white noise sigmas by a constant.
102    fn mul(mut self, rhs: f64) -> Self::Output {
103        self.sigma *= rhs;
104        self
105    }
106}
107
108impl MulAssign<f64> for WhiteNoise {
109    fn mul_assign(&mut self, rhs: f64) {
110        *self = *self * rhs;
111    }
112}
113
114#[cfg(test)]
115mod ut_wn {
116    use hifitime::{Epoch, TimeUnits};
117    use rand_pcg::Pcg64Mcg;
118
119    use super::{Stochastics, WhiteNoise};
120
121    #[test]
122    fn white_noise_test() {
123        let sigma = 10.0_f64;
124        let mut wn = WhiteNoise { mean: 0.0, sigma };
125
126        let mut larger_wn = WhiteNoise {
127            mean: 0.0,
128            sigma: sigma * 10.0,
129        };
130
131        let epoch = Epoch::now().unwrap();
132
133        let mut rng = Pcg64Mcg::new(1000);
134        let mut cnt_above_3sigma = 0;
135        let mut cnt_below_3sigma = 0;
136        let mut larger_cnt_above_3sigma = 0;
137        let mut larger_cnt_below_3sigma = 0;
138        for seconds in 0..1000_i64 {
139            let bias = wn.sample(epoch + seconds.seconds(), &mut rng);
140
141            if bias > 3.0 * sigma {
142                cnt_above_3sigma += 1;
143            } else if bias < -3.0 * sigma {
144                cnt_below_3sigma += 1;
145            }
146
147            let larger_bias = larger_wn.sample(epoch + seconds.seconds(), &mut rng);
148            if larger_bias > 30.0 * sigma {
149                larger_cnt_above_3sigma += 1;
150            } else if larger_bias < -30.0 * sigma {
151                larger_cnt_below_3sigma += 1;
152            }
153        }
154
155        assert!(dbg!(cnt_above_3sigma) <= 3);
156        assert!(dbg!(cnt_below_3sigma) <= 3);
157
158        assert!(dbg!(larger_cnt_above_3sigma) <= 3);
159        assert!(dbg!(larger_cnt_below_3sigma) <= 3);
160    }
161}
nyx_space/od/noise/white.rs

nyx_space/od/noise/
white.rs
