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 * (1.0 + pass).log10(),
36            stopband_db: -20.0 * stop.log10(),
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 * a.powf(n) + 2.0 * h * a.powf(n + 1.0))
53            })
54            .sum::<f64>()
55            / 4f64.powi(P);
56        let k = tv / ((2.0 * PI).powf(n + 1.0) * n);
57        let tail = |d: f64| k / d.powf(n);
58        let images = |from: f64| tail(1.0 - from) + k / ((n - 1.0) * (1.0 - from).powf(n - 1.0));
59        (
60            2.0 * tail(c - fp) + 4.0 * images(fp + c),
61            2.0 * tail(fs - c) + 4.0 * images(0.5 + c),
62        )
63    };
64    let mut h = 1.0;
65    loop {
66        let (pass, stop) = bounds(h);
67        if pass <= BOUND && stop <= BOUND {
68            return (h, pass, stop);
69        }
70        h += 1.0;
71    }
72}
73
74fn binomial(n: i32, k: i32) -> f64 {
75    (0..k).fold(1.0, |acc, i| acc * f64::from(n - i) / f64::from(i + 1))
76}
77
78#[derive(Debug)]
79pub struct Recorded {
80    planes: Arc<Vec<Vec<f32>>>,
81    rate: u32,
82    to: u32,
83    frames: usize,
84    cutoff: f64,
85    half: f64,
86}
87
88impl Recorded {
89    pub fn new(planes: Arc<Vec<Vec<f32>>>, rate: u32, to: u32) -> Recorded {
90        let frames = planes.iter().map(Vec::len).min().unwrap_or(0);
91        let lower = f64::from(rate.min(to)) / f64::from(rate);
92        let (half, ..) = design(lower);
93        Recorded {
94            planes,
95            rate,
96            to,
97            frames,
98            cutoff: (PASSBAND + STOPBAND) / 2.0 * lower,
99            half,
100        }
101    }
102
103    pub fn width(&self) -> usize {
104        self.planes.len()
105    }
106
107    pub fn resampled(&self) -> Option<Resampled> {
108        (self.rate != self.to).then(|| Resampled::of(self.rate, self.to))
109    }
110
111    pub fn secs(&self) -> f64 {
112        self.frames as f64 / f64::from(self.rate)
113    }
114
115    pub fn at(&self, c: usize, t: f64) -> f64 {
116        if !(0.0..self.secs()).contains(&t) {
117            return 0.0;
118        }
119        let plane = &self.planes[c];
120        let own = Grid::of(self.rate);
121        if self.rate == self.to
122            && let Some(n) = own.step_at(t, Round::Even)
123            && own.instant(n) == t
124        {
125            return usize::try_from(n)
126                .ok()
127                .and_then(|n| plane.get(n))
128                .map_or(0.0, |v| f64::from(*v));
129        }
130        let u = t * f64::from(self.rate);
131        let first = (u - self.half).ceil().max(0.0) as usize;
132        let last = ((u + self.half).floor() as usize).min(self.frames - 1);
133        (first..=last)
134            .map(|k| f64::from(plane[k]) * self.kernel(u - k as f64))
135            .sum()
136    }
137
138    fn kernel(&self, x: f64) -> f64 {
139        if x.abs() >= self.half {
140            return 0.0;
141        }
142        let y = 2.0 * self.cutoff * x;
143        let sinc = match y == 0.0 {
144            true => 1.0,
145            false => (PI * y).sin() / (PI * y),
146        };
147        2.0 * self.cutoff * sinc * (FRAC_PI_2 * x / self.half).cos().powi(2 * P)
148    }
149}
sva_samples/resample.rs

sva_samples/
resample.rs
