sva_samples/
resample.rs

1// Concern: a recording at any instant: its own sample at its own rate, else a windowed sinc of proven bound | Non-concern: decoding, which node reads it | IO: (planes, rates, t) -> a sample, a bound
2
3use std::f64::consts::{FRAC_PI_2, PI};
4use std::sync::Arc;
5
6use crate::grid::{Grid, Round};
7
8/// Relative to a tone's amplitude: -120 dB.
9pub const BOUND: f64 = 1e-6;
10/// Fractions of the lower rate.
11pub const PASSBAND: f64 = 0.45;
12pub const STOPBAND: f64 = 0.5;
13/// `cos^(2P)` windows the sinc.
14const P: i32 = 8;
15
16#[derive(Clone, Copy, Debug, PartialEq)]
17pub struct Resampled {
18    pub from_hz: u32,
19    pub to_hz: u32,
20    pub passband_hz: f64,
21    pub stopband_hz: f64,
22    pub ripple_db: f64,
23    pub stopband_db: f64,
24}
25
26impl Resampled {
27    pub fn of(from_hz: u32, to_hz: u32) -> Resampled {
28        let lower = f64::from(from_hz.min(to_hz));
29        let (_, pass, stop) = design(lower / f64::from(from_hz));
30        Resampled {
31            from_hz,
32            to_hz,
33            passband_hz: PASSBAND * lower,
34            stopband_hz: STOPBAND * lower,
35            ripple_db: 20.0 * libm::log10(1.0 + pass),
36            stopband_db: -20.0 * libm::log10(stop),
37        }
38    }
39}
40
41/// Poisson summation reads a tone at `f` back as `sum_m Hh(f - m)`, the passband box
42/// convolved with the window's transform `W`; `2P` integrations by parts give
43/// `|W(s)| <= TV(w^(2P)) / (2 pi s)^(2P+1)`, so `W`'s mass past `d` is at most `K / d^(2P)`.
44fn design(lower: f64) -> (f64, f64, f64) {
45    let (fp, fs) = (PASSBAND * lower, STOPBAND * lower);
46    let c = (fp + fs) / 2.0;
47    let n = f64::from(2 * P);
48    let bounds = |h: f64| {
49        let tv: f64 = (0..=2 * P)
50            .map(|j| {
51                let a = f64::from((P - j).abs()) * PI / h;
52                binomial(2 * P, j) * (2.0 * libm::pow(a, n) + 2.0 * h * libm::pow(a, n + 1.0))
53            })
54            .sum::<f64>()
55            / 4f64.powi(P);
56        let k = tv / (libm::pow(2.0 * PI, n + 1.0) * n);
57        let tail = |d: f64| k / libm::pow(d, n);
58        let images =
59            |from: f64| tail(1.0 - from) + k / ((n - 1.0) * libm::pow(1.0 - from, n - 1.0));
60        (
61            2.0 * tail(c - fp) + 4.0 * images(fp + c),
62            2.0 * tail(fs - c) + 4.0 * images(0.5 + c),
63        )
64    };
65    let mut h = 1.0;
66    loop {
67        let (pass, stop) = bounds(h);
68        if pass <= BOUND && stop <= BOUND {
69            return (h, pass, stop);
70        }
71        h += 1.0;
72    }
73}
74
75fn binomial(n: i32, k: i32) -> f64 {
76    (0..k).fold(1.0, |acc, i| acc * f64::from(n - i) / f64::from(i + 1))
77}
78
79#[derive(Debug)]
80pub struct Recorded {
81    planes: Arc<Vec<Vec<f32>>>,
82    rate: u32,
83    to: u32,
84    frames: usize,
85    cutoff: f64,
86    half: f64,
87}
88
89impl Recorded {
90    pub fn new(planes: Arc<Vec<Vec<f32>>>, rate: u32, to: u32) -> Recorded {
91        let frames = planes.iter().map(Vec::len).min().unwrap_or(0);
92        let lower = f64::from(rate.min(to)) / f64::from(rate);
93        let (half, ..) = design(lower);
94        Recorded {
95            planes,
96            rate,
97            to,
98            frames,
99            cutoff: (PASSBAND + STOPBAND) / 2.0 * lower,
100            half,
101        }
102    }
103
104    pub fn width(&self) -> usize {
105        self.planes.len()
106    }
107
108    pub fn resampled(&self) -> Option<Resampled> {
109        (self.rate != self.to).then(|| Resampled::of(self.rate, self.to))
110    }
111
112    pub fn secs(&self) -> f64 {
113        self.frames as f64 / f64::from(self.rate)
114    }
115
116    pub fn at(&self, c: usize, t: f64) -> f64 {
117        if !(0.0..self.secs()).contains(&t) {
118            return 0.0;
119        }
120        let plane = &self.planes[c];
121        let own = Grid::of(self.rate);
122        if self.rate == self.to
123            && let Some(n) = own.step_at(t, Round::Even)
124            && own.instant(n) == t
125        {
126            return usize::try_from(n)
127                .ok()
128                .and_then(|n| plane.get(n))
129                .map_or(0.0, |v| f64::from(*v));
130        }
131        let u = t * f64::from(self.rate);
132        let first = (u - self.half).ceil().max(0.0) as usize;
133        let last = ((u + self.half).floor() as usize).min(self.frames - 1);
134        (first..=last)
135            .map(|k| f64::from(plane[k]) * self.kernel(u - k as f64))
136            .sum()
137    }
138
139    fn kernel(&self, x: f64) -> f64 {
140        if x.abs() >= self.half {
141            return 0.0;
142        }
143        let y = 2.0 * self.cutoff * x;
144        let sinc = match y == 0.0 {
145            true => 1.0,
146            false => libm::sin(PI * y) / (PI * y),
147        };
148        2.0 * self.cutoff * sinc * libm::cos(FRAC_PI_2 * x / self.half).powi(2 * P)
149    }
150}
sva_samples/resample.rs

sva_samples/
resample.rs
