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#[derive(Copy, Clone, Debug, Default, PartialEq, Serialize, Deserialize, Sequence)]
36#[cfg_attr(feature = "python", pyclass(from_py_object, 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 white noise 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
88#[cfg(feature = "python")]
89#[cfg_attr(feature = "python", pymethods)]
90impl WhiteNoise {
91    #[new]
92    fn py_new(mean: f64, sigma: f64) -> Self {
93        Self { mean, sigma }
94    }
95
96    fn __str__(&self) -> String {
97        format!("{self:?}")
98    }
99
100    fn __repr__(&self) -> String {
101        format!("{self:?} @ {self:p}")
102    }
103}
104
105impl Stochastics for WhiteNoise {
106    fn covariance(&self, _epoch: Epoch) -> f64 {
107        self.sigma.powi(2)
108    }
109
110    fn sample<R: Rng>(&mut self, _epoch: Epoch, rng: &mut R) -> f64 {
111        rng.sample(Normal::new(self.mean, self.sigma).unwrap())
112    }
113}
114
115impl Mul<f64> for WhiteNoise {
116    type Output = Self;
117
118    /// Scale the white noise sigmas by a constant.
119    fn mul(mut self, rhs: f64) -> Self::Output {
120        self.sigma *= rhs;
121        self
122    }
123}
124
125impl MulAssign<f64> for WhiteNoise {
126    fn mul_assign(&mut self, rhs: f64) {
127        *self = *self * rhs;
128    }
129}
130
131#[cfg(test)]
132mod ut_wn {
133    use hifitime::{Epoch, TimeUnits};
134    use rand_pcg::Pcg64Mcg;
135
136    use super::{Stochastics, WhiteNoise};
137
138    #[test]
139    fn white_noise_test() {
140        let sigma = 10.0_f64;
141        let mut wn = WhiteNoise { mean: 0.0, sigma };
142
143        let mut larger_wn = WhiteNoise {
144            mean: 0.0,
145            sigma: sigma * 10.0,
146        };
147
148        let epoch = Epoch::now().unwrap();
149
150        let mut rng = Pcg64Mcg::new(1000);
151        let mut cnt_above_3sigma = 0;
152        let mut cnt_below_3sigma = 0;
153        let mut larger_cnt_above_3sigma = 0;
154        let mut larger_cnt_below_3sigma = 0;
155        for seconds in 0..1000_i64 {
156            let bias = wn.sample(epoch + seconds.seconds(), &mut rng);
157
158            if bias > 3.0 * sigma {
159                cnt_above_3sigma += 1;
160            } else if bias < -3.0 * sigma {
161                cnt_below_3sigma += 1;
162            }
163
164            let larger_bias = larger_wn.sample(epoch + seconds.seconds(), &mut rng);
165            if larger_bias > 30.0 * sigma {
166                larger_cnt_above_3sigma += 1;
167            } else if larger_bias < -30.0 * sigma {
168                larger_cnt_below_3sigma += 1;
169            }
170        }
171
172        assert!(dbg!(cnt_above_3sigma) <= 3);
173        assert!(dbg!(cnt_below_3sigma) <= 3);
174
175        assert!(dbg!(larger_cnt_above_3sigma) <= 3);
176        assert!(dbg!(larger_cnt_below_3sigma) <= 3);
177    }
178}
nyx_space/od/noise/white.rs

nyx_space/od/noise/
white.rs
