secs_dsp/
resample.rs

1//! Rate conversion for discrete impulse responses, preserving transfer gain.
2//!
3//! Independent windowed-sinc/polyphase implementation. Mathematical references:
4//! SciPy's resample_poly/upfirdn documentation and the symmetric four-term
5//! Blackman-Harris window (0.35875, 0.48829, 0.14128, 0.01168).
6//! [upfirdn](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.upfirdn.html)
7//! [Blackman-Harris](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.windows.blackmanharris.html)
8//!
9//! Unlike waveform resampling, the FIR gain is `down`, not `up`: an IR's
10//! discrete sum/transfer response must stay constant as sample density changes.
11//! Both filter tails are retained. The caller subtracts `leading / output_hz`
12//! from the returned samples' time origin.
13
14use std::f64::consts::PI;
15
16use crate::{DspError, invalid, pairwise_sum};
17
18pub const MIN_RESAMPLE_RATE: u32 = 8_000;
19pub const MAX_RESAMPLE_RATE: u32 = 384_000;
20const HALF_WIDTH: usize = 96;
21const CUTOFF: f64 = 0.96;
22
23/// Allocation-free shape and timing of a finite IR conversion.
24#[derive(Debug, Clone, Copy)]
25pub struct IrResamplePlan {
26    input_len: usize,
27    up: usize,
28    down: usize,
29    half: usize,
30    prepad: usize,
31    output_len: usize,
32    leading: usize,
33}
34
35impl IrResamplePlan {
36    pub fn new(input_len: usize, input_hz: u32, output_hz: u32) -> Result<Self, DspError> {
37        if input_len == 0 {
38            return Err(invalid("IR resampling needs nonempty input"));
39        }
40        if !(MIN_RESAMPLE_RATE..=MAX_RESAMPLE_RATE).contains(&input_hz)
41            || !(MIN_RESAMPLE_RATE..=MAX_RESAMPLE_RATE).contains(&output_hz)
42        {
43            return Err(invalid(
44                "IR resampling supports sample rates from 8000 to 384000 Hz",
45            ));
46        }
47        let mut a = input_hz;
48        let mut b = output_hz;
49        while b != 0 {
50            (a, b) = (b, a % b);
51        }
52        let up = (output_hz / a) as usize;
53        let down = (input_hz / a) as usize;
54        if up == down {
55            input_len
56                .checked_mul(size_of::<f64>())
57                .ok_or_else(|| invalid("IR resampling output byte size overflow"))?;
58            return Ok(Self {
59                input_len,
60                up,
61                down,
62                half: 0,
63                prepad: 0,
64                output_len: input_len,
65                leading: 0,
66            });
67        }
68        let half = HALF_WIDTH
69            .checked_mul(up.max(down))
70            .ok_or_else(|| invalid("IR resampling kernel length overflow"))?;
71        let prepad = (down - half % down) % down;
72        let last = (input_len - 1)
73            .checked_mul(up)
74            .and_then(|n| n.checked_add(2 * half + prepad))
75            .ok_or_else(|| invalid("IR resampling output length overflow"))?;
76        let output_len = (last / down)
77            .checked_add(1)
78            .ok_or_else(|| invalid("IR resampling output length overflow"))?;
79        output_len
80            .checked_mul(size_of::<f64>())
81            .ok_or_else(|| invalid("IR resampling output byte size overflow"))?;
82        (2 * half + 1)
83            .checked_mul(size_of::<f64>())
84            .ok_or_else(|| invalid("IR resampling kernel byte size overflow"))?;
85        Ok(Self {
86            input_len,
87            up,
88            down,
89            half,
90            prepad,
91            output_len,
92            leading: (half + prepad) / down,
93        })
94    }
95
96    pub fn output_len(&self) -> usize {
97        self.output_len
98    }
99
100    pub fn leading_output_samples(&self) -> usize {
101        self.leading
102    }
103
104    pub fn kernel_bytes(&self) -> usize {
105        if self.half == 0 {
106            0
107        } else {
108            (2 * self.half + 1) * size_of::<f64>()
109        }
110    }
111
112    /// Read samples without an intermediate f32-to-f64 allocation. The reader
113    /// supplies the same finite value for a given index throughout this call.
114    /// `check` allows the engine to propagate its own cancellation error.
115    pub fn process<E>(
116        &self,
117        sample: impl Fn(usize) -> f64,
118        mut check: impl FnMut() -> Result<(), E>,
119    ) -> Result<Vec<f64>, E>
120    where
121        E: From<DspError>,
122    {
123        check()?;
124        let mut output = Vec::new();
125        output
126            .try_reserve_exact(self.output_len)
127            .map_err(|_| E::from(invalid("IR resampling output allocation failed")))?;
128        if self.half == 0 {
129            for n in 0..self.input_len {
130                if n % 256 == 0 {
131                    check()?;
132                }
133                output.push(sample(n));
134            }
135            check()?;
136            return Ok(output);
137        }
138        let cutoff = CUTOFF / self.up.max(self.down) as f64;
139        let mut kernel = Vec::new();
140        kernel
141            .try_reserve_exact(2 * self.half + 1)
142            .map_err(|_| E::from(invalid("IR resampling kernel allocation failed")))?;
143        for i in 0..=2 * self.half {
144            if i % 256 == 0 {
145                check()?;
146            }
147            let x = cutoff * (i as f64 - self.half as f64);
148            let sinc = if x == 0.0 {
149                1.0
150            } else {
151                (PI * x).sin() / (PI * x)
152            };
153            let phase = PI * i as f64 / self.half as f64;
154            let window = 0.35875 - 0.48829 * phase.cos() + 0.14128 * (2.0 * phase).cos()
155                - 0.01168 * (3.0 * phase).cos();
156            kernel.push(cutoff * sinc * window);
157        }
158        check()?;
159        let scale = self.down as f64 / pairwise_sum(&kernel);
160        check()?;
161        for (i, value) in kernel.iter_mut().enumerate() {
162            if i % 256 == 0 {
163                check()?;
164            }
165            *value *= scale;
166        }
167        for k in 0..self.output_len {
168            if k % 256 == 0 {
169                check()?;
170            }
171            let time = k * self.down;
172            let mut value = 0.0;
173            if time >= self.prepad {
174                let time = time - self.prepad;
175                let first = time.saturating_sub(2 * self.half).div_ceil(self.up);
176                let last = (time / self.up).min(self.input_len - 1);
177                for n in first..=last {
178                    value += sample(n) * kernel[time - n * self.up];
179                }
180            }
181            if !value.is_finite() {
182                return Err(E::from(invalid(
183                    "IR resampling produced a nonfinite sample",
184                )));
185            }
186            output.push(value);
187        }
188        check()?;
189        Ok(output)
190    }
191}
secs_dsp/resample.rs

secs_dsp/
resample.rs
