sva_samples/measure/
pitch.rs

1// Concern: names each frame's spectral peaks as notes and flags likely harmonics | Non-concern: computing the spectrum (spectrum.rs), chord or key identity | IO: (planes, sample rate) -> frames
2
3use crate::measure::peaks::{Peak, resolved};
4use crate::measure::spectrum::magnitudes_of;
5use crate::measure::stiff::{fit, stretch};
6
7const NAMES: [&str; 12] = [
8    "C", "Cs", "D", "Ds", "E", "F", "Fs", "G", "Gs", "A", "As", "B",
9];
10pub const LOWEST_HZ: f64 = 25.0;
11pub const HIGHEST_HZ: f64 = 5000.0;
12const HARMONIC_TOLERANCE: f64 = 0.04;
13const HIGHEST_HARMONIC: f64 = 8.0;
14const COMPRESSION: f64 = 0.84;
15
16#[derive(Clone, Debug, PartialEq)]
17pub struct Note {
18    pub hz: f64,
19    pub name: String,
20    pub cents: f64,
21    pub db: f64,
22    pub harmonic_of: Option<f64>,
23}
24
25#[derive(Clone, Debug, PartialEq)]
26pub struct PitchFrame {
27    pub t_secs: f64,
28    pub notes: Vec<Note>,
29}
30
31pub fn track(
32    planes: &[&[f64]],
33    sample_rate: f64,
34    start_secs: f64,
35    frame_secs: f64,
36    max_notes: usize,
37) -> Vec<PitchFrame> {
38    let stride = ((frame_secs * sample_rate).round() as usize).max(1);
39    let len = crate::measure::held(planes);
40    (0..len)
41        .step_by(stride)
42        .enumerate()
43        .map(|(n, from)| {
44            let chunk: Vec<&[f64]> = planes
45                .iter()
46                .map(|p| &p[from..len.min(from + stride)])
47                .collect();
48            let framed = magnitudes_of(&chunk, sample_rate, None);
49            let found = resolved(&framed.mags, framed.bin_hz, framed.spanned, max_notes * 4);
50            let resolution_hz = sample_rate / chunk[0].len().max(1) as f64;
51            PitchFrame {
52                t_secs: start_secs + (n * stride) as f64 / sample_rate,
53                notes: name_all(&found, max_notes, resolution_hz / 2.0),
54            }
55        })
56        .collect()
57}
58
59/// The same naming over peaks a reading already holds exactly, loudest first.
60pub fn name_peaks(found: &[Peak], max_notes: usize) -> Vec<Note> {
61    name_all(found, max_notes, 0.0)
62}
63
64/// `spread_hz`: how far a peak may stand from where it was measured, 0 where exact.
65fn name_all(found: &[Peak], max_notes: usize, spread_hz: f64) -> Vec<Note> {
66    let listed: Vec<&Peak> = found
67        .iter()
68        .filter(|p| (LOWEST_HZ..=HIGHEST_HZ).contains(&p.hz))
69        .take(max_notes)
70        .collect();
71    let roots = series_roots(&listed, spread_hz);
72
73    listed
74        .iter()
75        .zip(roots)
76        .map(|(p, harmonic_of)| {
77            let midi = midi_of(p.hz);
78            let nearest = midi.round();
79            Note {
80                hz: p.hz,
81                name: note_name(nearest as i64),
82                cents: (midi - nearest) * 100.0,
83                db: p.db,
84                harmonic_of,
85            }
86        })
87        .collect()
88}
89
90pub fn midi_of(hz: f64) -> f64 {
91    69.0 + 12.0 * (hz / 440.0).log2()
92}
93
94pub fn note_name(midi: i64) -> String {
95    format!("{}{}", class_name(midi), midi.div_euclid(12) - 1)
96}
97
98pub fn class_name(midi: i64) -> &'static str {
99    NAMES[midi.rem_euclid(12) as usize]
100}
101
102/// Klapuri (ISMIR 2006): the greatest subharmonic sum (Hermes, JASA 83(1), 1988) roots a
103/// series; a note 2 to 8 times a louder member joins it, or the root if their pair's sum there
104/// reaches its amplitude. Exact peaks may stand as sharp as the series' stiff string fits.
105fn series_roots(notes: &[&Peak], spread_hz: f64) -> Vec<Option<f64>> {
106    let mut joined = vec![None; notes.len()];
107    let mut left: Vec<usize> = (0..notes.len()).collect();
108    while let Some((root, f0)) = most_salient(notes, &left, spread_hz) {
109        left.retain(|&i| i != root);
110        let mut series = vec![Member {
111            at: root,
112            n: 1,
113            hz: f0,
114        }];
115        loop {
116            let points: Vec<(u32, f64)> = series.iter().map(|m| (m.n, m.hz)).collect();
117            let b = match spread_hz > 0.0 {
118                true => 0.0,
119                // No fit, or a negative slope, is no stiffness: the plain whole-multiple window.
120                false => fit(&points).unwrap_or(0.0).max(0.0),
121            };
122            let joins = |m: &Member, i: usize| {
123                over(notes[i].hz, m, b).filter(|&q| {
124                    let (under, over) = (amp(notes[m.at].db), amp(notes[i].db));
125                    under > over
126                        || m.at == root && under + COMPRESSION.powi(q as i32 - 1) * over >= over
127                })
128            };
129            let Some((at, m, q)) = left
130                .iter()
131                .enumerate()
132                .find_map(|(at, &i)| series.iter().find_map(|m| joins(m, i).map(|q| (at, *m, q))))
133            else {
134                break;
135            };
136            let i = left.remove(at);
137            joined[i] = Some(notes[m.at].hz);
138            series.push(Member {
139                at: i,
140                n: q * m.n,
141                hz: notes[i].hz,
142            });
143        }
144    }
145    joined
146}
147
148#[derive(Clone, Copy)]
149struct Member {
150    at: usize,
151    n: u32,
152    hz: f64,
153}
154
155/// q from 2 to 8, `hz` lying from q times `m` to where stiffness `b` puts partial q n over n.
156fn over(hz: f64, m: &Member, b: f64) -> Option<u32> {
157    let ratio = hz / m.hz;
158    let q = ratio.round();
159    let (n, k) = (f64::from(m.n), q * f64::from(m.n));
160    let sharpest = q * stretch(k, b) / stretch(n, b);
161    ((2.0..=HIGHEST_HARMONIC).contains(&q)
162        && ratio > q - HARMONIC_TOLERANCE
163        && ratio < sharpest.max(q) + HARMONIC_TOLERANCE)
164        .then_some(q as u32)
165}
166
167fn amp(db: f64) -> f64 {
168    10f64.powf(db / 20.0)
169}
170
171fn salience(notes: &[&Peak], among: &[usize], hz: f64) -> f64 {
172    let mut strongest = [0.0f64; HIGHEST_HARMONIC as usize];
173    for &i in among {
174        if let Some(k) = multiple(notes[i].hz, hz) {
175            let at = &mut strongest[k as usize - 1];
176            *at = at.max(amp(notes[i].db));
177        }
178    }
179    strongest
180        .iter()
181        .enumerate()
182        .map(|(k, a)| COMPRESSION.powi(k as i32) * a)
183        .sum::<f64>()
184}
185
186fn most_salient(notes: &[&Peak], left: &[usize], spread_hz: f64) -> Option<(usize, f64)> {
187    left.iter()
188        .map(|&c| {
189            let f0 = fundamental(notes, left, c, spread_hz);
190            (c, f0, salience(notes, left, f0))
191        })
192        .max_by(|(a, _, sa), (b, _, sb)| {
193            sa.total_cmp(sb).then(notes[*b].hz.total_cmp(&notes[*a].hz))
194        })
195        .map(|(c, f0, _)| (c, f0))
196}
197
198/// Leakage pulls a measured peak further than stiffness: the fundamental in the root's spread
199/// explaining every note its peak does, and more, with the greatest sum.
200fn fundamental(notes: &[&Peak], left: &[usize], root: usize, spread_hz: f64) -> f64 {
201    let at = notes[root].hz;
202    if spread_hz <= 0.0 {
203        return at;
204    }
205    let explained = |f: f64| -> Vec<bool> {
206        left.iter()
207            .map(|&i| multiple(notes[i].hz, f).is_some())
208            .collect()
209    };
210    let own = explained(at);
211    let keeps = |f: f64| explained(f).iter().zip(&own).all(|(e, o)| *e || !o);
212    let near = |f: f64| (f - at).abs() <= spread_hz && (at / f - 1.0).abs() < HARMONIC_TOLERANCE;
213    left.iter()
214        .flat_map(|&i| (1..=HIGHEST_HARMONIC as u32).map(move |k| notes[i].hz / f64::from(k)))
215        .filter(|f| near(*f) && keeps(*f))
216        .map(|f| (f, salience(notes, left, f)))
217        .fold((at, salience(notes, left, at)), |held, next| {
218            let nearer = (held.0 - at).abs().total_cmp(&(next.0 - at).abs());
219            match next.1.total_cmp(&held.1).then(nearer) {
220                std::cmp::Ordering::Greater => next,
221                _ => held,
222            }
223        })
224        .0
225}
226
227fn multiple(hz: f64, root: f64) -> Option<u32> {
228    let n = hz / root;
229    let k = n.round();
230    ((1.0..=HIGHEST_HARMONIC).contains(&k) && (n - k).abs() < HARMONIC_TOLERANCE)
231        .then_some(k as u32)
232}
sva_samples/measure/pitch.rs

sva_samples/measure/
pitch.rs
