1#[derive(Clone, Copy, PartialEq, Debug)]
2pub enum Peak {
3 Low,
4 High,
5}
6
7pub struct PeaksDetector {
8 threshold: f64,
9 influence: f64,
10 window: Vec<f64>,
11}
12
13impl PeaksDetector {
14 pub fn new(lag: usize, threshold: f64, influence: f64) -> PeaksDetector {
15 PeaksDetector {
16 threshold,
17 influence,
18 window: Vec::with_capacity(lag),
19 }
20 }
21
22 pub fn signal(&mut self, value: f64) -> Option<Peak> {
23 if self.window.len() < self.window.capacity() {
24 self.window.push(value);
25 None
26 } else if let (Some((mean, stddev)), Some(&window_last)) = (self.stats(), self.window.last()) {
27 self.window.remove(0);
28 if (value - mean).abs() > (self.threshold * stddev) {
29 let next_value =
30 (value * self.influence) + ((1. - self.influence) * window_last);
31 self.window.push(next_value);
32 Some(if value > mean { Peak::High } else { Peak::Low })
33 } else {
34 self.window.push(value);
35 None
36 }
37 } else {
38 None
39 }
40 }
41
42 pub fn stats(&self) -> Option<(f64, f64)> {
43 if self.window.is_empty() {
44 None
45 } else {
46 let window_len = self.window.len() as f64;
47 let mean = self.window.iter().fold(0., |a, v| a + v) / window_len;
48 let sq_sum = self.window.iter().fold(0., |a, v| a + ((v - mean) * (v - mean)));
49 let stddev = (sq_sum / window_len).sqrt();
50 Some((mean, stddev))
51 }
52 }
53}
54
55pub struct PeaksIter<I, F> {
56 source: I,
57 signal: F,
58 detector: PeaksDetector,
59}
60
61pub trait PeaksFilter<I> where I: Iterator {
62 fn peaks<F>(self, detector: PeaksDetector, signal: F) -> PeaksIter<I, F>
63 where F: FnMut(&I::Item) -> f64;
64}
65
66impl<I> PeaksFilter<I> for I where I: Iterator {
67 fn peaks<F>(self, detector: PeaksDetector, signal: F) -> PeaksIter<I, F>
68 where F: FnMut(&I::Item) -> f64
69 {
70 PeaksIter { source: self, signal, detector, }
71 }
72}
73
74impl<I, F> Iterator for PeaksIter<I, F> where I: Iterator, F: FnMut(&I::Item) -> f64 {
75 type Item = (I::Item, Peak);
76
77 fn next(&mut self) -> Option<Self::Item> {
78 while let Some(item) = self.source.next() {
79 let value = (self.signal)(&item);
80 if let Some(peak) = self.detector.signal(value) {
81 return Some((item, peak));
82 }
83 }
84 None
85 }
86}
87
88#[cfg(test)]
89mod tests {
90 use super::{Peak, PeaksDetector, PeaksFilter};
91
92 #[test]
93 fn sample_data() {
94 let input = vec![
95 1.0, 1.0, 1.1, 1.0, 0.9, 1.0, 1.0, 1.1, 1.0, 0.9, 1.0, 1.1, 1.0, 1.0, 0.9, 1.0, 1.0, 1.1, 1.0,
96 1.0, 1.0, 1.0, 1.1, 0.9, 1.0, 1.1, 1.0, 1.0, 0.9, 1.0, 1.1, 1.0, 1.0, 1.1, 1.0, 0.8, 0.9, 1.0,
97 1.2, 0.9, 1.0, 1.0, 1.1, 1.2, 1.0, 1.5, 1.0, 3.0, 2.0, 5.0, 3.0, 2.0, 1.0, 1.0, 1.0, 0.9, 1.0,
98 1.0, 3.0, 2.6, 4.0, 3.0, 3.2, 2.0, 1.0, 1.0, 0.8, 4.0, 4.0, 2.0, 2.5, 1.0, 1.0, 1.0
99 ];
100 let output: Vec<_> = input
101 .into_iter()
102 .enumerate()
103 .peaks(PeaksDetector::new(30, 5.0, 0.0), |e| e.1)
104 .map(|((i, _), p)| (i, p))
105 .collect();
106 assert_eq!(output, vec![
107 (45, Peak::High), (47, Peak::High), (48, Peak::High), (49, Peak::High),
108 (50, Peak::High), (51, Peak::High), (58, Peak::High), (59, Peak::High),
109 (60, Peak::High), (61, Peak::High), (62, Peak::High), (63, Peak::High),
110 (67, Peak::High), (68, Peak::High), (69, Peak::High), (70, Peak::High),
111 ]);
112 }
113}