1use std::ops::Range;
10
11use crate::audio::SampleRate;
12use crate::cast;
13use crate::dsp::{ActiveFrames, StftParams, Window, frame_rms};
14use crate::units::{Decibels, Hz, MidiPitch, Seconds};
15
16use super::NoteList;
17
18pub const QUIET_FRAME: Seconds = Seconds::from_const(0.025);
20
21pub const QUIET_GATE: Decibels = crate::dsp::ACTIVE_GATE_DB;
24
25#[derive(Clone, Debug)]
29pub struct StemLevels<'a> {
30 mono: &'a [f64],
31 rate: SampleRate,
32 params: StftParams,
33 rms: Vec<f64>,
34 window: Vec<f64>,
35}
36
37impl<'a> StemLevels<'a> {
38 #[must_use]
40 pub fn new(mono: &'a [f64], rate: SampleRate) -> Self {
41 let frame = cast::round_to_usize(QUIET_FRAME.get() * f64::from(rate.get()))
42 .unwrap_or(4)
43 .max(4);
44 let params = StftParams::fixed(frame);
45 Self {
46 mono,
47 rate,
48 params,
49 rms: frame_rms(mono, params),
50 window: Window::Hann.samples(params.frame_len()),
51 }
52 }
53
54 #[must_use]
56 pub fn rms(&self) -> &[f64] {
57 &self.rms
58 }
59
60 pub(crate) fn ac_rms(&self, frame: usize) -> f64 {
63 self.ac_level(frame * self.params.hop(), self.params.frame_len())
64 }
65
66 pub(crate) fn frame_time(&self, frame: usize) -> Seconds {
68 Seconds::clamped(cast::usize_to_f64(frame * self.params.hop()) / f64::from(self.rate.get()))
69 }
70
71 fn hops(&self, time: Seconds) -> f64 {
72 time.get() * f64::from(self.rate.get()) / cast::usize_to_f64(self.params.hop())
73 }
74
75 fn sample(&self, centre: usize, offset: usize, len: usize) -> f64 {
78 (centre + offset)
79 .checked_sub(len / 2)
80 .and_then(|index| self.mono.get(index))
81 .copied()
82 .unwrap_or(0.0)
83 }
84
85 fn ac_level(&self, centre: usize, len: usize) -> f64 {
88 let (sum, energy) = (0..len)
89 .map(|offset| self.sample(centre, offset, len))
90 .fold((0.0, 0.0), |(sum, energy), x| (sum + x, energy + x * x));
91 let count = cast::usize_to_f64(len);
92 let mean = sum / count;
93 (energy / count - mean * mean).max(0.0).sqrt()
94 }
95
96 #[must_use]
100 pub fn envelope(&self, time: Seconds, window: Seconds) -> f64 {
101 let rate = f64::from(self.rate.get());
102 let len = cast::round_to_usize(window.get() * rate)
103 .unwrap_or(1)
104 .max(1);
105 self.ac_level(cast::round_to_usize(time.get() * rate).unwrap_or(0), len)
106 }
107
108 #[must_use]
110 pub fn frames_between(&self, from: Seconds, to: Seconds) -> Range<usize> {
111 let first = cast::ceil_to_usize(self.hops(from)).unwrap_or(0);
112 let end = cast::floor_to_usize(self.hops(to))
113 .map_or(0, |last| last + 1)
114 .min(self.rms.len());
115 first.min(end)..end
116 }
117
118 #[must_use]
127 pub fn harmonic_level(&self, frame: usize, pitch: MidiPitch) -> f64 {
128 let len = self.params.frame_len();
129 let rate = f64::from(self.rate.get());
130 let centre = frame * self.params.hop();
131 let raw: Vec<f64> = (0..len)
132 .map(|offset| self.sample(centre, offset, len))
133 .collect();
134 let mean = raw.iter().sum::<f64>() / cast::usize_to_f64(len);
135 let windowed: Vec<f64> = raw
136 .iter()
137 .zip(&self.window)
138 .map(|(x, weight)| (x - mean) * weight)
139 .collect();
140 let f0 = Hz::from(pitch).get();
141 (1..=HARMONICS)
142 .map(|h| f0 * f64::from(h))
143 .take_while(|&frequency| frequency < rate / 2.0)
144 .map(|frequency| bin_power(&windowed, frequency, rate))
145 .sum::<f64>()
146 .sqrt()
147 }
148
149 #[must_use]
160 pub fn fundamental_level(&self, time: Seconds, pitch: MidiPitch) -> f64 {
161 let rate = f64::from(self.rate.get());
162 let f0 = Hz::from(pitch).get();
163 let len = cast::round_to_usize(rate / f0).unwrap_or(1).max(1);
164 let centre = cast::round_to_usize(time.get() * rate).unwrap_or(0);
165 let raw: Vec<f64> = (0..len)
166 .map(|offset| self.sample(centre, offset, len))
167 .collect();
168 let mean = raw.iter().sum::<f64>() / cast::usize_to_f64(len);
169 let centred: Vec<f64> = raw.iter().map(|x| x - mean).collect();
170 bin_power(¢red, f0, rate).sqrt()
171 }
172}
173
174fn bin_power(samples: &[f64], frequency: f64, rate: f64) -> f64 {
177 let coefficient = 2.0 * (std::f64::consts::TAU * frequency / rate).cos();
178 let (mut s1, mut s2) = (0.0, 0.0);
179 for x in samples {
180 let s0 = x + coefficient * s1 - s2;
181 s2 = s1;
182 s1 = s0;
183 }
184 (s1 * s1 + s2 * s2 - coefficient * s1 * s2).max(0.0)
185}
186
187pub const HARMONICS: u8 = 8;
189
190#[must_use]
198pub fn drop_quiet(
199 notes: &NoteList,
200 mono: &[f64],
201 rate: SampleRate,
202 gate: Decibels,
203) -> (NoteList, usize) {
204 drop_quiet_with(notes, &StemLevels::new(mono, rate), gate)
205}
206
207#[must_use]
209pub fn drop_quiet_with(
210 notes: &NoteList,
211 levels: &StemLevels<'_>,
212 gate: Decibels,
213) -> (NoteList, usize) {
214 let active = ActiveFrames::from_rms(levels.rms(), gate);
215 let before = notes.notes().len();
216 let kept: Vec<_> = notes
217 .notes()
218 .iter()
219 .filter(|note| {
220 let first = cast::round_to_usize(levels.hops(note.onset())).unwrap_or(0);
221 let last = cast::ceil_to_usize(levels.hops(note.offset()))
222 .unwrap_or(first)
223 .max(first + 1);
224 (first..last).any(|m| active.is_active(m))
225 })
226 .cloned()
227 .collect();
228 let dropped = before - kept.len();
229 (NoteList::new(kept), dropped)
230}
231
232#[cfg(test)]
233mod tests {
234 use super::*;
235 use crate::notes::Note;
236 use crate::units::MidiPitch;
237
238 fn note(onset: f64, offset: f64) -> Note {
239 Note::new(
240 Seconds::new(onset).unwrap(),
241 Seconds::new(offset).unwrap(),
242 MidiPitch::new(60.0).unwrap(),
243 )
244 .unwrap()
245 }
246
247 #[test]
248 fn a_note_in_a_passage_50_db_down_is_dropped_and_one_at_30_db_kept() {
249 let rate = SampleRate::new(8000).unwrap();
250 let mono: Vec<f64> = (0..24_000)
252 .map(|n| {
253 let level = match n / 8000 {
254 0 => 1.0,
255 1 => 10_f64.powf(-50.0 / 20.0),
256 _ => 10_f64.powf(-30.0 / 20.0),
257 };
258 let t = f64::from(n) / 8000.0;
259 level * (std::f64::consts::TAU * 220.0 * t).sin()
260 })
261 .collect();
262 let notes = NoteList::new(vec![
263 note(0.2, 0.8),
264 note(1.2, 1.8),
265 note(2.2, 2.8),
266 note(0.9, 1.5),
267 ]);
268 let (kept, dropped) = drop_quiet(¬es, &mono, rate, QUIET_GATE);
269 assert_eq!(dropped, 1);
270 let onsets: Vec<f64> = kept.notes().iter().map(|n| n.onset().get()).collect();
271 assert_eq!(onsets, vec![0.2, 0.9, 2.2]);
272 }
273
274 fn silence_then_tone() -> Vec<f64> {
277 (0..16_000)
278 .map(|n| {
279 let t = f64::from(n) / 8000.0;
280 if n < 8000 {
281 0.0
282 } else {
283 (std::f64::consts::TAU * 220.0 * t).sin()
284 }
285 })
286 .collect()
287 }
288
289 #[test]
290 fn a_note_that_reaches_into_a_loud_passage_is_kept() {
291 let rate = SampleRate::new(8000).unwrap();
292 let notes = NoteList::new(vec![note(0.5, 1.03), note(0.5, 0.9)]);
293 let (kept, dropped) = drop_quiet(¬es, &silence_then_tone(), rate, QUIET_GATE);
294 assert_eq!(dropped, 1);
295 assert!((kept.notes()[0].offset().get() - 1.03).abs() < 1e-9);
296 }
297
298 #[test]
299 fn a_note_in_a_10_ms_gap_is_kept_by_the_25_ms_frames_over_it() {
300 let rate = SampleRate::new(8000).unwrap();
301 let tone = silence_then_tone().split_off(8000);
303 let mut mono = tone.clone();
304 mono.extend(std::iter::repeat_n(0.0, 80));
305 mono.extend(&tone);
306 let notes = NoteList::new(vec![note(1.001, 1.008)]);
307 let (kept, dropped) = drop_quiet(¬es, &mono, rate, QUIET_GATE);
308 assert_eq!((kept.notes().len(), dropped), (1, 0));
309 }
310
311 #[test]
312 fn a_note_shorter_than_a_hop_takes_the_frame_at_its_onset() {
313 let rate = SampleRate::new(8000).unwrap();
314 let hop = 0.006_25;
317 let notes = NoteList::new(vec![note(240.6 * hop, 240.9 * hop)]);
318 let (kept, dropped) = drop_quiet(¬es, &silence_then_tone(), rate, QUIET_GATE);
319 assert_eq!((kept.notes().len(), dropped), (1, 0));
320 }
321
322 fn sine_from(rate: u32, frequency: f64, start: usize) -> Vec<f64> {
325 (0..rate)
326 .map(|n| {
327 let t = f64::from(n) / f64::from(rate);
328 if usize::try_from(n).unwrap() < start {
329 0.0
330 } else {
331 (std::f64::consts::TAU * frequency * t).sin()
332 }
333 })
334 .collect()
335 }
336
337 #[test]
338 fn frames_between_holds_both_ends() {
339 let rate = SampleRate::new(8000).unwrap();
340 let mono = vec![0.0; 8000];
341 let levels = StemLevels::new(&mono, rate);
342 let s = |value| Seconds::new(value).unwrap();
343 assert_eq!(levels.frames_between(s(0.0), s(0.0125)), 0..3);
344 assert_eq!(levels.frames_between(s(0.001), s(0.0125)), 1..3);
345 assert_eq!(
346 levels.frames_between(s(0.99), s(2.0)).end,
347 levels.rms().len()
348 );
349 }
350
351 #[test]
352 fn the_harmonic_level_reads_each_harmonic_in_its_own_frame() {
353 let rate = SampleRate::new(8000).unwrap();
356 let a4 = MidiPitch::new(69.0).unwrap();
357 for frequency in [440.0, 1320.0] {
358 let mono = sine_from(8000, frequency, 4000);
359 let levels = StemLevels::new(&mono, rate);
360 assert!(levels.harmonic_level(78, a4).abs() < 1e-12);
362 let level = levels.harmonic_level(82, a4);
363 assert!((level - 50.0).abs() < 0.5, "{frequency} Hz reads {level}");
364 }
365 }
366
367 #[test]
368 fn the_envelope_reads_the_rms_of_its_window_about_its_mean() {
369 let rate = SampleRate::new(8000).unwrap();
372 let s = |value| Seconds::new(value).unwrap();
373 let mono = sine_from(8000, 440.0, 4000);
374 let offset: Vec<f64> = mono.iter().map(|x| x + 0.5).collect();
375 for stem in [&mono, &offset] {
376 let levels = StemLevels::new(stem, rate);
377 let full = levels.envelope(s(0.75), s(0.01));
378 assert!(
379 (full - std::f64::consts::FRAC_1_SQRT_2).abs() < 0.02,
380 "{full}"
381 );
382 }
383 let half = StemLevels::new(&mono, rate).envelope(s(0.5), s(0.01));
385 assert!((half - 0.5).abs() < 0.03, "{half}");
386 }
387
388 #[test]
389 fn a_dc_offset_leaves_the_harmonic_level_alone() {
390 let rate = SampleRate::new(8000).unwrap();
391 let a1 = MidiPitch::new(33.0).unwrap();
392 let mono = sine_from(8000, 55.0, 0);
393 let offset: Vec<f64> = mono.iter().map(|x| x + 0.5).collect();
394 let clean = StemLevels::new(&mono, rate).harmonic_level(80, a1);
395 let shifted = StemLevels::new(&offset, rate).harmonic_level(80, a1);
396 assert!(
397 (shifted / clean - 1.0).abs() < 0.05,
398 "{shifted} against {clean}"
399 );
400 }
401
402 #[test]
403 fn below_a2_the_harmonic_level_hears_a_fifth_as_the_note() {
404 let rate = SampleRate::new(8000).unwrap();
407 let a1 = MidiPitch::new(33.0).unwrap();
408 let level = |frequency| {
409 let mono = sine_from(8000, frequency, 0);
410 StemLevels::new(&mono, rate).harmonic_level(80, a1)
411 };
412 let ratio = level(82.407) / level(55.0);
413 assert!(
414 ratio > 10_f64.powf(-3.0 / 20.0),
415 "the fifth reads {ratio} of A1"
416 );
417 }
418
419 #[test]
420 fn the_harmonic_level_stops_below_the_nyquist_frequency() {
421 let nyquist: Vec<f64> = (0..7040)
423 .map(|n| if n % 2 == 0 { 1.0 } else { -1.0 })
424 .collect();
425 let levels = StemLevels::new(&nyquist, SampleRate::new(7040).unwrap());
426 assert!(levels.harmonic_level(80, MidiPitch::new(69.0).unwrap()) < 1.0);
427 let pitch = MidiPitch::new(69.0 + 12.0 * (1250.0_f64 / 440.0).log2()).unwrap();
430 let mono = sine_from(8000, 1750.0, 0);
431 let levels = StemLevels::new(&mono, SampleRate::new(8000).unwrap());
432 assert!(levels.harmonic_level(80, pitch) < 1.0);
433 }
434
435 #[test]
436 fn silence_drops_every_note() {
437 let rate = SampleRate::new(8000).unwrap();
438 let (kept, dropped) = drop_quiet(
439 &NoteList::new(vec![note(0.1, 0.2)]),
440 &vec![0.0; 8000],
441 rate,
442 QUIET_GATE,
443 );
444 assert!(kept.is_empty());
445 assert_eq!(dropped, 1);
446 }
447}