harmos_signal/waveform/
calculus.rs

1//! The rate a waveform changes at, and the quantity it accumulates to.
2
3/// The derivative of `y` with respect to `x`, sample for sample.
4///
5/// Central differences inside, one-sided differences at the two ends, so the
6/// answer is as long as the run it came from. A run of fewer than two samples
7/// has no rate of change.
8#[must_use]
9pub fn derivative(x: &[f64], y: &[f64]) -> Vec<f64> {
10    let length = x.len().min(y.len());
11    if length < 2 {
12        return Vec::new();
13    }
14
15    let mut slope = Vec::with_capacity(length);
16    slope.push((y[1] - y[0]) / (x[1] - x[0]));
17    slope.extend(
18        (1..length - 1).map(|index| (y[index + 1] - y[index - 1]) / (x[index + 1] - x[index - 1])),
19    );
20    slope.push((y[length - 1] - y[length - 2]) / (x[length - 1] - x[length - 2]));
21
22    slope
23}
24
25/// The cumulative integral of `y` over `x` by the trapezoidal rule.
26///
27/// As long as the run it came from and starting at zero, so sample `i` is the
28/// area under the waveform up to `x[i]` — a charge from a current, a distance
29/// from a speed, an energy from a power.
30#[must_use]
31pub fn integral(x: &[f64], y: &[f64]) -> Vec<f64> {
32    let length = x.len().min(y.len());
33    if length == 0 {
34        return Vec::new();
35    }
36
37    let mut accumulated = 0.0;
38    let mut area = Vec::with_capacity(length);
39    area.push(0.0);
40
41    for index in 1..length {
42        accumulated += (y[index] + y[index - 1]) / 2.0 * (x[index] - x[index - 1]);
43        area.push(accumulated);
44    }
45
46    area
47}
48
49#[cfg(test)]
50mod tests {
51    use super::*;
52
53    /// A unit-spaced axis of `count` coordinates.
54    fn axis(count: usize) -> Vec<f64> {
55        (0..count).map(|index| index as f64).collect()
56    }
57
58    #[test]
59    fn a_straight_line_differentiates_to_its_slope() {
60        let x = axis(100);
61        let y: Vec<f64> = x.iter().map(|at| 2.0 * at).collect();
62
63        let slope = derivative(&x, &y);
64
65        assert_eq!(slope.len(), 100);
66        assert!(slope.iter().all(|slope| (slope - 2.0).abs() < 1e-12));
67    }
68
69    #[test]
70    fn a_constant_integrates_to_a_ramp() {
71        let x = axis(11);
72        let area = integral(&x, &[2.0; 11]);
73
74        assert!((area[0] - 0.0).abs() < 1e-12);
75        assert!((area[5] - 10.0).abs() < 1e-12);
76        assert!((area[10] - 20.0).abs() < 1e-12);
77    }
78
79    #[test]
80    fn a_triangle_integrates_to_its_own_area() {
81        // A ramp from zero to ten over ten seconds encloses fifty.
82        let x = axis(11);
83        let y: Vec<f64> = x.clone();
84
85        let area = integral(&x, &y);
86
87        assert!((area[10] - 50.0).abs() < 1e-12);
88    }
89
90    #[test]
91    fn a_derivative_integrated_back_recovers_the_shape() {
92        let x: Vec<f64> = (0..100).map(|index| f64::from(index) * 0.1).collect();
93        let y: Vec<f64> = x.iter().map(|at| at.sin()).collect();
94
95        let recovered = integral(&x, &derivative(&x, &y));
96
97        // The constant of integration is gone, so the shape is what survives.
98        for index in 10..90 {
99            let original = y[index + 1] - y[index];
100            let round_trip = recovered[index + 1] - recovered[index];
101            assert!(
102                (original - round_trip).abs() < 1e-3,
103                "sample {index} drifted"
104            );
105        }
106    }
107
108    #[test]
109    fn a_run_too_short_to_differentiate_has_no_answer() {
110        assert!(derivative(&[], &[]).is_empty());
111        assert!(derivative(&[1.0], &[1.0]).is_empty());
112        assert!(integral(&[], &[]).is_empty());
113        assert_eq!(integral(&[1.0], &[5.0]), vec![0.0]);
114    }
115}
harmos_signal/waveform/calculus.rs

harmos_signal/waveform/
calculus.rs
