pamoja_kit/
median.rs

1//! A rolling median filter for rejecting spikes.
2
3/// A median filter over the most recent `N` readings.
4///
5/// A single bad sample - a spike from electrical noise or a flaky contact - drags a mean or
6/// an exponential average off course, because it is blended into the result. The median
7/// ignores it: one outlier cannot move the middle value of a sorted window. That makes a
8/// [`Median`] the right filter when the noise is occasional spikes rather than steady
9/// jitter; for steady jitter reach for [`Smoother`](crate::Smoother) instead. Keep the
10/// window small and odd (3, 5, 7) so there is a single middle reading; with an even `N` the
11/// median is the average of the two middle readings.
12///
13/// # Examples
14///
15/// ```
16/// use pamoja_kit::Median;
17///
18/// let mut filtered = Median::<5>::new();
19/// // A lone spike among steady readings is rejected.
20/// for reading in [10.0, 10.0, 99.0, 10.0, 10.0] {
21///     filtered.update(reading);
22/// }
23/// assert_eq!(filtered.median(), Some(10.0));
24/// ```
25#[derive(Clone, Copy, Debug)]
26pub struct Median<const N: usize> {
27    samples: [f32; N],
28    len: usize,
29    next: usize,
30}
31
32impl<const N: usize> Median<N> {
33    /// Creates an empty median filter over a window of `N` readings.
34    ///
35    /// # Returns
36    ///
37    /// A filter holding no readings yet.
38    pub fn new() -> Self {
39        Self {
40            samples: [0.0; N],
41            len: 0,
42            next: 0,
43        }
44    }
45
46    /// Adds a reading and returns the median of the current window.
47    ///
48    /// # Arguments
49    ///
50    /// * `reading` - the latest raw reading.
51    ///
52    /// # Returns
53    ///
54    /// The median of the readings now in the window. With a zero-length window (`N` is `0`)
55    /// the reading passes through unchanged.
56    pub fn update(&mut self, reading: f32) -> f32 {
57        self.push(reading);
58        self.median().unwrap_or(reading)
59    }
60
61    /// Adds a reading to the window, evicting the oldest once it is full.
62    ///
63    /// # Arguments
64    ///
65    /// * `reading` - the latest raw reading.
66    pub fn push(&mut self, reading: f32) {
67        if N == 0 {
68            return;
69        }
70        self.samples[self.next] = reading;
71        self.next = (self.next + 1) % N;
72        if self.len < N {
73            self.len += 1;
74        }
75    }
76
77    /// Returns the median of the readings in the window, or [`None`] if it is empty.
78    ///
79    /// # Returns
80    ///
81    /// The middle reading of the sorted window for an odd count, the average of the two
82    /// middle readings for an even count, or [`None`] before any reading.
83    pub fn median(&self) -> Option<f32> {
84        if self.len == 0 {
85            return None;
86        }
87        let mut sorted = [0.0f32; N];
88        sorted[..self.len].copy_from_slice(&self.samples[..self.len]);
89        let window = &mut sorted[..self.len];
90        window.sort_unstable_by(f32::total_cmp);
91        let mid = self.len / 2;
92        if self.len % 2 == 1 {
93            Some(window[mid])
94        } else {
95            Some((window[mid - 1] + window[mid]) / 2.0)
96        }
97    }
98
99    /// Returns the number of readings currently held, at most `N`.
100    pub fn len(&self) -> usize {
101        self.len
102    }
103
104    /// Returns `true` if the filter holds no readings.
105    pub fn is_empty(&self) -> bool {
106        self.len == 0
107    }
108}
109
110impl<const N: usize> Default for Median<N> {
111    fn default() -> Self {
112        Self::new()
113    }
114}
115
116#[cfg(test)]
117mod tests {
118    use super::*;
119
120    #[test]
121    fn odd_window_takes_the_middle_value() {
122        let mut median = Median::<3>::new();
123        median.update(3.0);
124        median.update(1.0);
125        assert_eq!(median.update(2.0), 2.0); // sorted [1, 2, 3] -> 2
126    }
127
128    #[test]
129    fn even_window_averages_the_two_middle_values() {
130        // The standard convention: the median of {1, 2, 3, 4} is (2 + 3) / 2 = 2.5.
131        let mut median = Median::<4>::new();
132        for reading in [1.0, 2.0, 3.0, 4.0] {
133            median.push(reading);
134        }
135        assert_eq!(median.median(), Some(2.5));
136    }
137
138    #[test]
139    fn a_single_spike_is_rejected() {
140        let mut median = Median::<5>::new();
141        for reading in [10.0, 10.0, 99.0, 10.0, 10.0] {
142            median.push(reading);
143        }
144        assert_eq!(median.median(), Some(10.0));
145    }
146
147    #[test]
148    fn an_empty_filter_has_no_median() {
149        let median = Median::<3>::new();
150        assert!(median.is_empty());
151        assert_eq!(median.median(), None);
152    }
153
154    #[test]
155    fn the_window_evicts_oldest_when_full() {
156        let mut median = Median::<3>::new();
157        for reading in [1.0, 2.0, 3.0, 100.0, 100.0] {
158            median.push(reading);
159        }
160        // The window holds the last three readings, 3, 100, 100: median 100.
161        assert_eq!(median.median(), Some(100.0));
162        assert_eq!(median.len(), 3);
163    }
164}
pamoja_kit/median.rs

pamoja_kit/
median.rs
