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

sva_samples/measure/
pitch.rs
