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

sva_samples/measure/
pitch.rs
