secs_dsp/
numeric.rs

1//! Complex division adapted from Impulcifer-py313
2//! `crates/impulcifer-dsp/src/stages/room_tuning2/peq.rs` (MIT,
3//! Copyright (c) 2024- 115dkk); see docs/THIRD_PARTY.md.
4//! NumPy arithmetic order used by the primitives.
5use crate::Complex64;
6
7/// NumPy's contiguous float64 pairwise reduction (128-element blocks).
8pub fn sum(x: &[f64]) -> f64 {
9    if x.len() < 8 {
10        x.iter().fold(-0.0, |a, b| a + b)
11    } else if x.len() <= 128 {
12        let mut r: [f64; 8] = x[..8].try_into().unwrap();
13        let stop = x.len() / 8 * 8;
14        for i in (8..stop).step_by(8) {
15            for j in 0..8 {
16                r[j] += x[i + j];
17            }
18        }
19        let mut s = ((r[0] + r[1]) + (r[2] + r[3])) + ((r[4] + r[5]) + (r[6] + r[7]));
20        for v in &x[stop..] {
21            s += v;
22        }
23        s
24    } else {
25        let n = (x.len() / 2) / 8 * 8;
26        sum(&x[..n]) + sum(&x[n..])
27    }
28}
29
30// Adapted from Impulcifer-py313 crates/impulcifer-dsp/src/stages/room_tuning2/peq.rs
31// (MIT, Copyright (c) 2024- 115dkk); see docs/THIRD_PARTY.md.
32/// NumPy's complex division (Smith's scaled formula, including a division by
33/// zero), so quotients match NumPy to the last bit.
34pub fn divide(a: Complex64, b: Complex64) -> Complex64 {
35    if b.re.abs() >= b.im.abs() {
36        if b.re == 0.0 && b.im == 0.0 {
37            return Complex64::new(a.re / b.re.abs(), a.im / b.re.abs());
38        }
39        let r = b.im / b.re;
40        let d = 1. / (b.re + b.im * r);
41        Complex64::new((a.re + a.im * r) * d, (a.im - a.re * r) * d)
42    } else {
43        let r = b.re / b.im;
44        let d = 1. / (b.re * r + b.im);
45        Complex64::new((a.re * r + a.im) * d, (a.im * r - a.re) * d)
46    }
47}
48
49/// NumPy linear interpolation with endpoint clamping (nonempty sorted grid).
50pub fn interp(x: f64, g: &[f64], v: &[f64]) -> f64 {
51    if x.is_nan() {
52        return f64::NAN;
53    }
54    if x < g[0] {
55        return v[0];
56    }
57    if x > g[g.len() - 1] {
58        return v[v.len() - 1];
59    }
60    let k = g.partition_point(|a| *a <= x).saturating_sub(1);
61    if k == g.len() - 1 || x == g[k] {
62        return v[k];
63    }
64    let slope = (v[k + 1] - v[k]) / (g[k + 1] - g[k]);
65    let y = slope * (x - g[k]) + v[k];
66    if !y.is_nan() {
67        return y;
68    }
69    let y = slope * (x - g[k + 1]) + v[k + 1];
70    if y.is_nan() && v[k] == v[k + 1] {
71        v[k]
72    } else {
73        y
74    }
75}
secs_dsp/numeric.rs

secs_dsp/
numeric.rs
