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

nyx_space/od/noise/
white.rs
