harmos_signal/waveform/
resample.rs

1//! Reading a waveform between the samples it actually has.
2//!
3//! One interpolation, used four ways: at a point, onto an axis a caller
4//! supplies, onto a uniform grid, and onto a target sample count. It is linear
5//! and it *clamps* — a query outside the sampled span answers the nearest
6//! sample rather than extrapolating a trend that was never measured.
7
8/// The value at `at`, linearly interpolated between the samples bracketing it.
9///
10/// `x` must be sorted ascending. Outside the sampled span the nearest end
11/// sample is answered; an empty run answers zero.
12#[must_use]
13pub fn value_at(x: &[f64], y: &[f64], at: f64) -> f64 {
14    let length = x.len().min(y.len());
15    match length {
16        0 => return 0.0,
17        1 => return y[0],
18        _ => {}
19    }
20
21    if at <= x[0] {
22        return y[0];
23    }
24    if at >= x[length - 1] {
25        return y[length - 1];
26    }
27
28    let above = x[..length]
29        .partition_point(|sample| *sample < at)
30        .clamp(1, length - 1);
31    let below = above - 1;
32    let span = x[above] - x[below];
33
34    if span.abs() < f64::EPSILON {
35        return (y[below] + y[above]) / 2.0;
36    }
37
38    y[below] + (at - x[below]) / span * (y[above] - y[below])
39}
40
41/// The waveform read at each coordinate of `at`.
42#[must_use]
43pub fn resample(x: &[f64], y: &[f64], at: &[f64]) -> Vec<f64> {
44    at.iter().map(|at| value_at(x, y, *at)).collect()
45}
46
47/// The waveform read onto a uniform grid of `step`, spanning what it covers.
48///
49/// The grid starts at `x[0]` and takes as many steps as fit inside the sampled
50/// span, so the answer is a uniformly sampled copy of the same waveform. A
51/// non-positive step, or an empty run, resamples to nothing.
52#[must_use]
53pub fn resample_uniform(x: &[f64], y: &[f64], step: f64) -> Vec<f64> {
54    let length = x.len().min(y.len());
55    if length == 0 || step <= 0.0 {
56        return Vec::new();
57    }
58
59    let span = x[length - 1] - x[0];
60    let count = (span / step).floor() as usize + 1;
61
62    (0..count)
63        .map(|index| value_at(x, y, x[0] + index as f64 * step))
64        .collect()
65}
66
67/// The run stretched or squeezed to exactly `count` samples.
68///
69/// Index space, with no axis at all: the first and last samples are preserved
70/// and everything between is interpolated evenly. This is the reading a
71/// characterization curve wants, where the samples are ordered but the
72/// coordinate they sit on is nobody's business.
73#[must_use]
74pub fn resample_count(values: &[f64], count: usize) -> Vec<f64> {
75    if values.is_empty() || count == 0 {
76        return Vec::new();
77    }
78    if values.len() == 1 || count == 1 {
79        return vec![values[0]; count];
80    }
81
82    let step = (values.len() - 1) as f64 / (count - 1) as f64;
83    let axis: Vec<f64> = (0..values.len()).map(|index| index as f64).collect();
84
85    (0..count)
86        .map(|index| value_at(&axis, values, index as f64 * step))
87        .collect()
88}
89
90/// The index of the sample nearest `at`, when there is one.
91///
92/// `x` must be sorted ascending. A query outside the sampled span answers the
93/// end nearest it.
94#[must_use]
95pub fn nearest(x: &[f64], at: f64) -> Option<usize> {
96    if x.is_empty() {
97        return None;
98    }
99
100    let above = x.partition_point(|sample| *sample < at);
101    if above == 0 {
102        return Some(0);
103    }
104    if above >= x.len() {
105        return Some(x.len() - 1);
106    }
107
108    let below = above - 1;
109    Some(if at - x[below] <= x[above] - at {
110        below
111    } else {
112        above
113    })
114}
115
116#[cfg(test)]
117mod tests {
118    use super::*;
119
120    #[test]
121    fn a_midpoint_reads_the_average_of_its_neighbours() {
122        let x = vec![0.0, 1.0, 2.0];
123        let y = vec![0.0, 10.0, 20.0];
124
125        assert!((value_at(&x, &y, 0.5) - 5.0).abs() < 1e-12);
126        assert!((value_at(&x, &y, 1.5) - 15.0).abs() < 1e-12);
127        assert!((value_at(&x, &y, 1.0) - 10.0).abs() < 1e-12);
128    }
129
130    #[test]
131    fn a_quarter_of_the_way_reads_a_quarter_of_the_rise() {
132        assert!((value_at(&[0.0, 4.0], &[0.0, 100.0], 1.0) - 25.0).abs() < 1e-12);
133        assert!((value_at(&[0.0, 4.0], &[0.0, 100.0], 3.0) - 75.0).abs() < 1e-12);
134    }
135
136    #[test]
137    fn outside_the_span_the_nearest_end_is_the_answer() {
138        let x = vec![0.0, 1.0, 2.0];
139        let y = vec![0.0, 10.0, 20.0];
140
141        assert_eq!(value_at(&x, &y, -5.0), 0.0);
142        assert_eq!(value_at(&x, &y, 10.0), 20.0);
143    }
144
145    #[test]
146    fn a_run_too_short_to_interpolate_answers_what_it_has() {
147        assert_eq!(value_at(&[], &[], 5.0), 0.0);
148        assert_eq!(value_at(&[1.0], &[10.0], 5.0), 10.0);
149    }
150
151    #[test]
152    fn resampling_onto_an_axis_reads_it_point_by_point() {
153        let x = vec![0.0, 1.0, 2.0];
154        let y = vec![0.0, 10.0, 20.0];
155
156        assert_eq!(resample(&x, &y, &[0.0, 0.5, 2.0]), vec![0.0, 5.0, 20.0]);
157        assert!(resample(&x, &y, &[]).is_empty());
158    }
159
160    #[test]
161    fn a_uniform_grid_covers_the_span_it_was_given() {
162        // Ten samples a second apart, read back twice as often.
163        let x: Vec<f64> = (0..10).map(f64::from).collect();
164        let y: Vec<f64> = x.iter().map(|at| at * 2.0).collect();
165
166        let resampled = resample_uniform(&x, &y, 0.5);
167
168        assert_eq!(resampled.len(), 19);
169        // The waveform is a straight line, so every reading lands on it.
170        for (index, value) in resampled.iter().enumerate() {
171            assert!(
172                (value - index as f64).abs() < 1e-12,
173                "sample {index} read {value}"
174            );
175        }
176    }
177
178    #[test]
179    fn a_grid_with_no_step_resamples_to_nothing() {
180        assert!(resample_uniform(&[0.0, 1.0], &[0.0, 1.0], 0.0).is_empty());
181        assert!(resample_uniform(&[], &[], 0.5).is_empty());
182    }
183
184    #[test]
185    fn a_count_resample_keeps_both_ends() {
186        assert_eq!(
187            resample_count(&[0.0, 100.0], 5),
188            vec![0.0, 25.0, 50.0, 75.0, 100.0]
189        );
190        assert_eq!(
191            resample_count(&[0.0, 25.0, 50.0, 75.0, 100.0], 3),
192            vec![0.0, 50.0, 100.0]
193        );
194        assert_eq!(resample_count(&[0.0, 10.0, 20.0], 3), vec![0.0, 10.0, 20.0]);
195    }
196
197    #[test]
198    fn a_count_resample_of_nothing_is_nothing() {
199        assert!(resample_count(&[], 10).is_empty());
200        assert!(resample_count(&[1.0, 2.0], 0).is_empty());
201        assert_eq!(resample_count(&[5.0], 3), vec![5.0, 5.0, 5.0]);
202        assert_eq!(resample_count(&[1.0, 2.0, 3.0], 1), vec![1.0]);
203    }
204
205    #[test]
206    fn the_nearest_sample_is_the_one_the_query_is_closest_to() {
207        let x = vec![0.0, 1.0, 2.0, 3.0, 4.0];
208
209        assert_eq!(nearest(&x, 0.0), Some(0));
210        assert_eq!(nearest(&x, 1.3), Some(1));
211        assert_eq!(nearest(&x, 1.7), Some(2));
212        assert_eq!(nearest(&x, 2.5), Some(2));
213        assert_eq!(nearest(&x, -1.0), Some(0));
214        assert_eq!(nearest(&x, 10.0), Some(4));
215        assert_eq!(nearest(&[], 1.0), None);
216    }
217
218    #[test]
219    fn the_nearest_sample_of_an_uneven_axis_is_still_the_closest() {
220        let x = vec![0.0, 0.3, 1.0];
221
222        assert_eq!(nearest(&x, 0.1), Some(0));
223        assert_eq!(nearest(&x, 0.2), Some(1));
224        assert_eq!(nearest(&x, 0.5), Some(1));
225        assert_eq!(nearest(&x, 0.8), Some(2));
226    }
227}
harmos_signal/waveform/resample.rs

harmos_signal/waveform/
resample.rs
