sva_samples/measure/score/
voices.rs

1// Concern: what sounds in one span: voices by a harmonic sieve, periodicity, coherence, a hit's kind | Non-concern: where a span starts (strikes.rs) | IO: (Clip, peaks or span) -> voices, period, kind
2
3use crate::fft::{dft, fft, idft, ifft};
4
5use super::spectra::{Clip, cents, db_of_power};
6
7const HARMONIC_SHARE: f64 = 0.035;
8const HARMONIC_CENTS: f64 = 40.0;
9const START_DB: f64 = -30.0;
10const MOST_VOICES: usize = 6;
11const OVER_DB: f64 = 6.0;
12const LEAST_SHARE: f64 = 0.01;
13const LONE_SHARE: f64 = 0.15;
14const PERIOD_HZ: (f64, f64) = (40.0, 1500.0);
15const PERIOD_CLOSE: f64 = 0.9;
16const COHERENCE_LAG_SECS: f64 = 0.01;
17const COHERENCE_LEAST_HZ: f64 = 20.0;
18const COHERENCE_SHARE: f64 = 0.03;
19const HIT_SAMPLES_48K: f64 = 1024.0;
20const KICK_HZ: f64 = 150.0;
21const HAT_HZ: f64 = 5000.0;
22
23#[derive(Clone, Debug, PartialEq)]
24pub struct Partial {
25    pub hz: f64,
26    pub db: f64,
27    pub harmonic: usize,
28}
29
30#[derive(Clone, Debug, PartialEq)]
31pub struct Voice {
32    pub hz: f64,
33    pub db: f64,
34    pub share: f64,
35    pub partials: Vec<Partial>,
36}
37
38pub fn harmonic(hz: f64, f0: f64) -> usize {
39    let k = (hz / f0).round_ties_even();
40    let ok = k >= 1.0
41        && (hz / f0 - k).abs() <= HARMONIC_SHARE
42        && cents(hz, k * f0).abs() <= HARMONIC_CENTS;
43    if ok { k as usize } else { 0 }
44}
45
46struct Open {
47    f0: f64,
48    ks: Vec<usize>,
49    partials: Vec<Partial>,
50    loudest: f64,
51}
52
53/// The lowest peak no voice explains starts one, a heard fundamental; a higher peak joins
54/// the first voice whose series it continues, k f0 with k - 1 or k - 2 in it, unless it
55/// stands 6 dB over that voice's loudest, a loud partial over a weak sub being its own voice.
56pub fn voices(peaks: &[(f64, f64)]) -> Vec<Voice> {
57    if peaks.is_empty() {
58        return Vec::new();
59    }
60    let top = peaks.iter().map(|p| p.1).fold(f64::NEG_INFINITY, f64::max);
61    let mut by_hz = peaks.to_vec();
62    by_hz.sort_by(|a, b| a.0.total_cmp(&b.0));
63    let mut open: Vec<Open> = Vec::new();
64    for (hz, db) in by_hz {
65        let home = open.iter().position(|v| {
66            let k = harmonic(hz, v.f0);
67            k >= 2
68                && (v.ks.contains(&(k - 1)) || v.ks.contains(&(k - 2)) || k <= 2)
69                && db <= v.loudest + OVER_DB
70        });
71        match home {
72            Some(i) => {
73                let v = &mut open[i];
74                let k = harmonic(hz, v.f0);
75                if !v.ks.contains(&k) {
76                    v.ks.push(k);
77                }
78                v.partials.push(Partial {
79                    hz,
80                    db,
81                    harmonic: k,
82                });
83                v.loudest = v.loudest.max(db);
84            }
85            None if db > top + START_DB && open.len() < MOST_VOICES => open.push(Open {
86                f0: hz,
87                ks: vec![1],
88                partials: vec![Partial {
89                    hz,
90                    db,
91                    harmonic: 1,
92                }],
93                loudest: db,
94            }),
95            None => {}
96        }
97    }
98    let total: f64 = peaks.iter().map(|p| libm::pow(10.0, p.1 / 10.0)).sum();
99    let mut out: Vec<Voice> = open
100        .into_iter()
101        .filter_map(|v| {
102            let w = |p: &Partial| libm::pow(10.0, p.db / 10.0);
103            let held: f64 = v.partials.iter().map(w).sum();
104            let k = |p: &Partial| p.harmonic as f64;
105            let num: f64 = v.partials.iter().map(|p| w(p) * p.hz * k(p)).sum();
106            let den: f64 = v.partials.iter().map(|p| w(p) * k(p) * k(p)).sum();
107            let share = held / total;
108            let lone = v.partials.len() < 2 && share < LONE_SHARE;
109            (share >= LEAST_SHARE && !lone).then(|| Voice {
110                hz: num / den,
111                db: db_of_power(held),
112                share,
113                partials: v.partials,
114            })
115        })
116        .collect();
117    out.sort_by(|a, b| b.share.total_cmp(&a.share));
118    out
119}
120
121/// McLeod's normalised square difference over lags for 40 Hz to 1.5 kHz: the first peak
122/// within 0.9 of the highest, the shortest period that repeats; its frequency and height.
123pub fn period(seg: &[f64], sr: f64) -> Option<(f64, f64)> {
124    let n = seg.len();
125    let most = ((sr / PERIOD_HZ.0) as usize).min(n / 2);
126    let least = (sr / PERIOD_HZ.1) as usize;
127    if most < least + 3 {
128        return None;
129    }
130    let size = (2 * n).next_power_of_two();
131    let (mut re, mut im) = (vec![0.0; size], vec![0.0; size]);
132    re[..n].copy_from_slice(seg);
133    fft(&mut re, &mut im);
134    for k in 0..size {
135        re[k] = re[k] * re[k] + im[k] * im[k];
136        im[k] = 0.0;
137    }
138    ifft(&mut re, &mut im);
139    let mut c = Vec::with_capacity(n);
140    let mut run = 0.0;
141    for v in seg {
142        run += v * v;
143        c.push(run);
144    }
145    let d: Vec<f64> = (0..most)
146        .map(|lag| {
147            let m = c[n - 1] - if lag > 0 { c[lag - 1] } else { 0.0 } + c[n - 1 - lag];
148            2.0 * re[lag] / (m + 1e-20)
149        })
150        .collect();
151    let peaks: Vec<usize> = (least.max(1)..most - 1)
152        .filter(|&i| d[i] >= d[i - 1] && d[i] > d[i + 1] && d[i] > 0.0)
153        .collect();
154    let top = peaks
155        .iter()
156        .map(|&i| d[i])
157        .fold(f64::NEG_INFINITY, f64::max);
158    let i = *peaks.iter().find(|&&i| d[i] >= PERIOD_CLOSE * top)?;
159    Some((sr / i as f64, d[i]))
160}
161
162/// A tone's band stays near 1 against itself 10 ms on; noise of bandwidth B falls as
163/// e^(-pi B tau).
164pub fn coherence(clip: &Clip, (a, b): (f64, f64), partials: &[&Partial]) -> Option<f64> {
165    let (a, b) = (clip.at(a), clip.at(b));
166    let seg = &clip.x[a..b.max(a)];
167    let lag = (COHERENCE_LAG_SECS * clip.sr) as usize;
168    if seg.len() < 3 * lag || lag == 0 {
169        return None;
170    }
171    let n = seg.len();
172    let (re, im) = dft(seg, &vec![0.0; n]);
173    let (mut sum, mut weights) = (0.0, 0.0);
174    for p in partials {
175        let half = COHERENCE_LEAST_HZ.max(COHERENCE_SHARE * p.hz);
176        let (mut zr, mut zi) = (vec![0.0; n], vec![0.0; n]);
177        for k in 0..n.div_ceil(2) {
178            let f = k as f64 * clip.sr / n as f64;
179            if f > p.hz - half && f < p.hz + half {
180                zr[k] = 2.0 * re[k];
181                zi[k] = 2.0 * im[k];
182            }
183        }
184        let (zr, zi) = idft(&zr, &zi);
185        let (mut cr, mut ci, mut e0, mut e1) = (0.0, 0.0, 0.0, 0.0);
186        for m in 0..seg.len() - lag {
187            let (ar, ai, br, bi) = (zr[m], zi[m], zr[m + lag], zi[m + lag]);
188            cr += ar * br + ai * bi;
189            ci += ar * bi - ai * br;
190            e0 += ar * ar + ai * ai;
191            e1 += br * br + bi * bi;
192        }
193        let w = libm::pow(10.0, p.db / 10.0);
194        sum += w * (cr * cr + ci * ci).sqrt() / ((e0 * e1).sqrt() + 1e-30);
195        weights += w;
196    }
197    (weights > 0.0).then(|| sum / weights)
198}
199
200#[derive(Clone, Copy, Debug, PartialEq, Eq)]
201pub enum Hit {
202    Kick,
203    Snare,
204    Hat,
205}
206
207impl Hit {
208    pub fn key(self) -> &'static str {
209        match self {
210            Hit::Kick => "kick",
211            Hit::Snare => "snare",
212            Hit::Hat => "hat",
213        }
214    }
215}
216
217pub fn hit(clip: &Clip, (a, b): (usize, usize)) -> Option<Hit> {
218    let n = super::spectra::sized(HIT_SAMPLES_48K, clip.sr);
219    let p = clip.spectrum_of(a, b, n)?;
220    let total = p.iter().sum::<f64>() + 1e-30;
221    let hz = |k: usize| k as f64 * clip.sr / n as f64;
222    let low: f64 = (0..p.len())
223        .filter(|&k| hz(k) < KICK_HZ)
224        .map(|k| p[k])
225        .sum();
226    let high: f64 = (0..p.len())
227        .filter(|&k| hz(k) >= HAT_HZ)
228        .map(|k| p[k])
229        .sum();
230    Some(match (low / total > 0.5, high / total > 0.4) {
231        (true, _) => Hit::Kick,
232        (false, true) => Hit::Hat,
233        _ => Hit::Snare,
234    })
235}
sva_samples/measure/score/voices.rs

sva_samples/measure/score/
voices.rs
