sva_samples/measure/harmony/
measured.rs

1// Concern: the partials of samples, a Hann frame centred on each instant, leakage peaks left out | Non-concern: picking peaks (peaks.rs) | IO: (planes, rate, start, hop) -> moments, framing
2
3use super::{Framing, Moment, Partial, Window};
4use crate::measure::peaks::{kernel, leakage, peaks};
5use crate::measure::spectrum::{magnitudes_of, pinned_frame};
6
7pub const WINDOW_SECS: f64 = 0.05;
8
9pub const METHOD: &str = "measured: a Hann frame of `window.size` samples of each channel \
10                          centred on each instant, zero past its ends, their power per bin \
11                          averaged over the channels; its local \
12                          maxima, each refined by the two-bin Hann \
13                          interpolation of Grandke, T. (1983), IEEE Trans. Instrum. Meas. 32(2) \
14                          350-355, exact for a lone sinusoid; a peak no louder than the Hann \
15                          kernel's bound 1/(pi m (m^2 - 1)) at m bins (Harris, F. J. (1978), \
16                          Proc. IEEE 66(1) 51-83) lets the louder peaks kept leak there is \
17                          dropped as leakage. Each partial's bounds are Grandke's estimate over \
18                          its two bins each widened by what that bound lets the other kept \
19                          partials, at their estimated levels, leak into it: they hold where the \
20                          frame is stationary sinusoids at the kept partials, and are a point \
21                          for a lone one";
22
23fn offset(a: f64) -> f64 {
24    match a.is_finite() {
25        true => ((2.0 * a - 1.0) / (a + 1.0)).clamp(-1.0, 1.0),
26        false => 1.0,
27    }
28}
29
30/// A peak's own bin and the side its louder neighbour stands on.
31struct Peaked {
32    k: usize,
33    sign: f64,
34    hz: f64,
35    amplitude: f64,
36}
37
38/// Grandke's offset d = (2a - 1)/(a + 1) toward the louder neighbour, a = |X(k+-1)|/|X(k)|,
39/// and amplitude |X(k)| over the kernel there.
40fn refined(mags: &[f64], bin_hz: f64, hz: f64) -> Peaked {
41    let near = (hz / bin_hz).round() as usize;
42    let k = (near.saturating_sub(1)..=near + 1)
43        .filter(|k| (1..mags.len() - 1).contains(k))
44        .max_by(|a, b| mags[*a].total_cmp(&mags[*b]))
45        .expect("a peak stands off both ends");
46    let sign = match mags[k + 1] >= mags[k - 1] {
47        true => 1.0,
48        false => -1.0,
49    };
50    let d = offset(mags[(k as f64 + sign) as usize] / mags[k]);
51    Peaked {
52        k,
53        sign,
54        hz: (k as f64 + sign * d) * bin_hz,
55        amplitude: mags[k] / kernel(d),
56    }
57}
58
59/// The same estimate over each of its two bins widened by `leaked`, what the other kept
60/// partials can put there: the bounds METHOD states.
61fn bounded(p: &Peaked, mags: &[f64], bin_hz: f64, leaked: impl Fn(f64) -> f64) -> Partial {
62    let toward = p.k as f64 + p.sign;
63    let widened = |x: f64, by: f64| ((x - by).max(0.0), x + by);
64    let own = widened(mags[p.k], leaked(p.k as f64));
65    let side = widened(mags[toward as usize], leaked(toward));
66    let a_hi = match own.0 > 0.0 {
67        true => side.1 / own.0,
68        false => f64::INFINITY,
69    };
70    let d = (offset(side.0 / own.1), offset(a_hi));
71    let ends = [
72        (p.k as f64 + p.sign * d.0) * bin_hz,
73        (p.k as f64 + p.sign * d.1) * bin_hz,
74    ];
75    let nearest = match d.0 <= 0.0 && d.1 >= 0.0 {
76        true => 0.0,
77        false => d.0.abs().min(d.1.abs()),
78    };
79    Partial {
80        hz: p.hz,
81        amplitude: p.amplitude,
82        hz_bounds: (ends[0].min(ends[1]), ends[0].max(ends[1])),
83        amplitude_bounds: (
84            own.0 / kernel(nearest),
85            own.1 / kernel(d.0.abs().max(d.1.abs())),
86        ),
87    }
88}
89
90pub fn moments(
91    planes: &[&[f64]],
92    sample_rate: f64,
93    start_secs: f64,
94    hop_secs: f64,
95) -> (Vec<Moment>, Framing) {
96    let len = crate::measure::held(planes);
97    let hop = ((hop_secs * sample_rate).round() as usize).max(1);
98    let size = pinned_frame(WINDOW_SECS, sample_rate);
99    let half = (size / 2) as i64;
100    let mut held = Vec::new();
101    for centre in (0..len).step_by(hop) {
102        let frames: Vec<Vec<f64>> = planes
103            .iter()
104            .map(|samples| {
105                (0..size as i64)
106                    .map(|k| {
107                        let at = centre as i64 - half + k;
108                        usize::try_from(at)
109                            .ok()
110                            .and_then(|at| samples.get(at))
111                            .copied()
112                            .unwrap_or(0.0)
113                    })
114                    .collect()
115            })
116            .collect();
117        let frames: Vec<&[f64]> = frames.iter().map(Vec::as_slice).collect();
118        let framed = magnitudes_of(&frames, sample_rate, Some(WINDOW_SECS));
119        let (mags, bin_hz) = (framed.mags, framed.bin_hz);
120        let mut kept: Vec<Peaked> = Vec::new();
121        for peak in peaks(&mags, bin_hz, usize::MAX) {
122            let found = refined(&mags, bin_hz, peak.hz);
123            let leaked: f64 = kept
124                .iter()
125                .map(|q| q.amplitude * leakage((found.hz - q.hz).abs() / bin_hz))
126                .sum();
127            if found.amplitude > leaked {
128                kept.push(found);
129            }
130        }
131        let partials = kept
132            .iter()
133            .enumerate()
134            .map(|(i, p)| {
135                let leaked = |bin: f64| -> f64 {
136                    kept.iter()
137                        .enumerate()
138                        .filter(|(j, _)| *j != i)
139                        .map(|(_, q)| q.amplitude * leakage((bin - q.hz / bin_hz).abs()))
140                        .sum()
141                };
142                bounded(p, &mags, bin_hz, leaked)
143            })
144            .collect();
145        held.push(Moment {
146            t_secs: start_secs + centre as f64 / sample_rate,
147            partials,
148        });
149    }
150    let framing = Framing {
151        hop_secs: hop as f64 / sample_rate,
152        window: Some(Window {
153            size,
154            secs: size as f64 / sample_rate,
155            resolution_hz: sample_rate / size as f64,
156        }),
157    };
158    (held, framing)
159}
sva_samples/measure/harmony/measured.rs

sva_samples/measure/harmony/
measured.rs
