sva_samples/measure/
beating.rs

1// Concern: how far each band's level swings about its own trend per sliding window | Non-concern: specific loudness (loudness_pattern.rs) | IO: (planes, rate, start, hop, window, spl) -> Beating
2
3use super::loudness_pattern::{self, LoudnessPattern};
4use super::sliding;
5
6pub const MODEL: &str = "per window of `window_secs` centred on each frame of \
7                         `loudness_pattern` at the same `frame` and `spl`, clipped to the \
8                         target at its ends, `from`/`to` its span: each band's level y = 10 \
9                         log2(N' + T) over the window's frames, T the pattern's \
10                         `threshold_sone`; the frames within 20 phon of the band's loudest \
11                         there, at least 20 of them, fitted by least squares with a quadratic \
12                         in time, its trend; the band's swing the rms of y about that trend; \
13                         `phon` the bands' swings averaged, each weighted by its mean specific \
14                         loudness over the window, a band no frame of the window holds over T \
15                         left out, null where no band is left. A quadratic follows a note's \
16                         rise and its fall in phon, so what remains moves faster than the \
17                         window's length, beating and tremolo, up to the pattern's window \
18                         rate; a swing between two partials inside one band reads as the \
19                         beating it is heard as. A window holding under 20 frames refuses.";
20
21/// The fewest frames a band's trend is fitted over: a correctness limit, not a cost one, so
22/// it has no override; a quadratic through fewer follows the swing it should leave.
23pub const LEAST_FRAMES: usize = 20;
24/// How far under a band's loudest a frame may fall and still count: below, a note's silence
25/// before or after it, not its swing.
26const RANGE_PHON: f64 = 20.0;
27
28#[derive(Clone, Copy, Debug, PartialEq)]
29pub struct BeatingWindow {
30    pub t_secs: f64,
31    pub from_secs: f64,
32    pub to_secs: f64,
33    pub phon: Option<f64>,
34}
35
36#[derive(Clone, Debug, PartialEq)]
37pub struct Beating {
38    pub frame_secs: f64,
39    pub window_secs: f64,
40    pub spl_db: f64,
41    pub windows: Vec<BeatingWindow>,
42}
43
44/// The window holds fewer than `LEAST_FRAMES` frames; `least_secs` is the least that holds
45/// them at this frame.
46#[derive(Clone, Copy, Debug, PartialEq)]
47pub struct ShortWindow {
48    pub window_secs: f64,
49    pub frame_secs: f64,
50    pub least_secs: f64,
51}
52
53pub fn measure(
54    planes: &[&[f64]],
55    sr: f64,
56    start: f64,
57    (hop_secs, window_secs): (f64, f64),
58    spl_db: f64,
59) -> Result<Beating, ShortWindow> {
60    let pattern = loudness_pattern::measure(planes, sr, start, hop_secs, spl_db);
61    let rate = 1.0 / pattern.frame_secs;
62    if ((window_secs * rate).round() as usize) < LEAST_FRAMES {
63        return Err(ShortWindow {
64            window_secs,
65            frame_secs: pattern.frame_secs,
66            least_secs: LEAST_FRAMES as f64 * pattern.frame_secs,
67        });
68    }
69    let grid = sliding::spans(
70        super::held(planes),
71        sr,
72        start,
73        pattern.frame_secs,
74        rate,
75        window_secs,
76        pattern.frames.len(),
77    );
78    let windows = grid
79        .spans
80        .iter()
81        .map(|s| BeatingWindow {
82            t_secs: s.t_secs,
83            from_secs: s.from_secs,
84            to_secs: s.to_secs,
85            phon: swing(&pattern, s.from..s.to),
86        })
87        .collect();
88    Ok(Beating {
89        frame_secs: pattern.frame_secs,
90        window_secs: grid.window_secs,
91        spl_db,
92        windows,
93    })
94}
95
96fn swing(pattern: &LoudnessPattern, frames: std::ops::Range<usize>) -> Option<f64> {
97    let held = &pattern.frames[frames];
98    let floor = pattern.threshold_sone;
99    let (mut sum, mut weights) = (0.0, 0.0);
100    for z in 0..pattern.centres_hz.len() {
101        if held.iter().all(|f| f.specific[z] <= floor) {
102            continue;
103        }
104        let y: Vec<f64> = held
105            .iter()
106            .map(|f| 10.0 * libm::log2(f.specific[z] + floor))
107            .collect();
108        let top = y.iter().copied().fold(f64::NEG_INFINITY, f64::max);
109        let kept: Vec<(f64, f64)> = y
110            .iter()
111            .enumerate()
112            .filter(|(_, v)| **v > top - RANGE_PHON)
113            .map(|(i, v)| (2.0 * i as f64 / held.len() as f64 - 1.0, *v))
114            .collect();
115        if kept.len() < LEAST_FRAMES {
116            continue;
117        }
118        let w = held.iter().map(|f| f.specific[z]).sum::<f64>() / held.len() as f64;
119        sum += about_quadratic(&kept) * w;
120        weights += w;
121    }
122    (weights > 0.0).then(|| sum / weights)
123}
124
125/// The rms of `points` about their least-squares quadratic, `t` on [-1, 1) so the normal
126/// equations stay as well conditioned at a minute's frames as at a second's.
127fn about_quadratic(points: &[(f64, f64)]) -> f64 {
128    let mut s = [0.0; 5];
129    let mut r = [0.0; 3];
130    for &(t, y) in points {
131        let mut p = 1.0;
132        for (k, slot) in s.iter_mut().enumerate() {
133            *slot += p;
134            if k < 3 {
135                r[k] += p * y;
136            }
137            p *= t;
138        }
139    }
140    let m = [[s[0], s[1], s[2]], [s[1], s[2], s[3]], [s[2], s[3], s[4]]];
141    let c = solve3(m, r);
142    let square: f64 = points
143        .iter()
144        .map(|&(t, y)| {
145            let d = y - (c[0] + c[1] * t + c[2] * t * t);
146            d * d
147        })
148        .sum();
149    (square / points.len() as f64).sqrt()
150}
151
152/// Cramer's rule.
153fn solve3(m: [[f64; 3]; 3], r: [f64; 3]) -> [f64; 3] {
154    let det = |a: [[f64; 3]; 3]| {
155        a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1])
156            - a[0][1] * (a[1][0] * a[2][2] - a[1][2] * a[2][0])
157            + a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0])
158    };
159    let d = det(m);
160    let mut out = [0.0; 3];
161    for (k, slot) in out.iter_mut().enumerate() {
162        let mut a = m;
163        for row in 0..3 {
164            a[row][k] = r[row];
165        }
166        *slot = det(a) / d;
167    }
168    out
169}
sva_samples/measure/beating.rs

sva_samples/measure/
beating.rs
