diff --git a/crates/dsp/src/offline.rs b/crates/dsp/src/offline.rs index fcda93b..8bbc80a 100644 --- a/crates/dsp/src/offline.rs +++ b/crates/dsp/src/offline.rs @@ -379,73 +379,143 @@ pub fn signal_generator(kind: Waveform, freq: f32, level_db: f32, sample_rate: f (0..channels).map(|_| mono.clone()).collect() } -/// Onset detection by half-wave-rectified spectral flux with an adaptive threshold. -/// `sensitivity` 0..1 (higher finds more transients). Returns sample positions, ascending. +/// Onset detection: how much a log-magnitude spectrum on 24 bands per octave rises against the frame +/// before, each band against the loudest of its neighbours there so vibrato doesn't count +/// (vibrato-suppressed spectral flux). A peak must stand out from the flux within 0.1 s by a fixed +/// amount, not relative to the loudest hit, so quiet hits count and the result doesn't depend on where +/// the analysed audio starts. Each is placed where its attack starts. `sensitivity` 0..1 (higher finds +/// more). Returns sample positions, ascending. pub fn detect_transients(ch: &[Vec], sample_rate: f32, sensitivity: f32) -> Vec { let len = max_len(ch); let sr = sane_sr(sample_rate); let sens = if sensitivity.is_finite() { sensitivity.clamp(0.0, 1.0) } else { 0.5 }; - let n = ((sr * 0.023) as usize).clamp(64, 8192).next_power_of_two(); + let n = ((sr * 0.012) as usize).clamp(64, 8192).next_power_of_two(); let hop = n / 4; if len < n || ch.is_empty() { return Vec::new(); } + let ch: Vec> = + ch.iter().take(crate::util::MAX_CHANNELS).map(|c| c.iter().map(|v| if v.is_finite() { *v } else { 0.0 }).collect()).collect(); let nch = ch.len() as f32; - let mono: Vec = (0..len).map(|i| ch.iter().map(|c| c.get(i).copied().filter(|v| v.is_finite()).unwrap_or(0.0)).sum::() / nch).collect(); - let mut planner = FftPlanner::::new(); - let fft = planner.plan_fft_forward(n); + let bands = log_bands(n, sr); + let fft = FftPlanner::::new().plan_fft_forward(n); let window: Vec = (0..n).map(|i| (0.5 - 0.5 * (2.0 * PI * i as f64 / n as f64).cos()) as f32).collect(); + // Magnitudes relative to a full-scale sine; below about -80 dB a band's changes barely count. + let gain = 1.0e4 / (n as f32 / 4.0); let frames = (len - n) / hop + 1; - let mut prev = vec![0.0f32; n / 2 + 1]; let mut buf = vec![Complex::new(0.0f32, 0.0); n]; + let mut power = vec![0.0f32; n / 2 + 1]; + let mut prev: Vec = Vec::new(); let mut flux = Vec::with_capacity(frames); + let mut energy: Vec = Vec::with_capacity(frames); for m in 0..frames { - let start = m * hop; - for (i, b) in buf.iter_mut().enumerate() { - *b = Complex::new(mono.get(start + i).copied().unwrap_or(0.0) * window.get(i).copied().unwrap_or(0.0), 0.0); - } - fft.process(&mut buf); - let mut f = 0.0; - for (b, p) in buf.iter().zip(prev.iter_mut()) { - let mag = (1.0 + 100.0 * b.norm()).ln(); - f += (mag - *p).max(0.0); - *p = mag; - } - flux.push(if m == 0 { 0.0 } else { f }); - } - let max_flux = flux.iter().copied().fold(0.0f32, f32::max); - if max_flux <= 1e-6 { - return Vec::new(); - } - let w = ((0.1 * sr) as usize / hop).max(2); - let mult = 1.0 + 2.5 * (1.0 - sens); - let delta = max_flux * (0.02 + 0.25 * (1.0 - sens)); - let min_gap = ((0.05 * sr) as usize / hop).max(1); - let mut out = Vec::new(); - let mut last: Option = None; - for m in 1..flux.len() { - let cur = flux.get(m).copied().unwrap_or(0.0); - let lo = m.saturating_sub(w); - let hi = (m + w + 1).min(flux.len()); - let local = flux.get(lo..hi).map(|s| s.iter().sum::() / s.len().max(1) as f32).unwrap_or(0.0); - let prev_v = flux.get(m - 1).copied().unwrap_or(0.0); - let next_v = flux.get(m + 1).copied().unwrap_or(0.0); - if cur > local * mult + delta && cur >= prev_v && cur > next_v && last.is_none_or(|l| m - l >= min_gap) { - last = Some(m); - // Refine: first sample in the analysis frame reaching a quarter of the frame's peak. - let start = m * hop; - let seg = mono.get(start..(start + n).min(len)).unwrap_or(&[]); - let pk = seg.iter().fold(0.0f32, |a, v| a.max(v.abs())); - let off = seg.iter().position(|v| v.abs() >= 0.25 * pk).unwrap_or(n / 2); - let pos = start + off; - if out.last().is_none_or(|&p: &usize| pos > p) { - out.push(pos); + power.iter_mut().for_each(|p| *p = 0.0); + for c in &ch { + for (i, b) in buf.iter_mut().enumerate() { + *b = Complex::new(c.get(m * hop + i).copied().unwrap_or(0.0) * window.get(i).copied().unwrap_or(0.0), 0.0); } + fft.process(&mut buf); + for (p, b) in power.iter_mut().zip(&buf) { + *p += b.norm_sqr() / nch; + } + } + energy.push(power.iter().sum()); + let cur: Vec = bands + .iter() + .map(|band| (1.0 + gain * band.iter().map(|&(k, w)| power.get(k).map_or(0.0, |p| p.sqrt()) * w).sum::()).log10()) + .collect(); + let rise = |b: usize, v: f32| { + (v - prev.get(b.saturating_sub(1)..(b + 2).min(prev.len())).map_or(0.0, |s| s.iter().copied().fold(0.0, f32::max))).max(0.0) + }; + flux.push(if prev.is_empty() { 0.0 } else { cur.iter().enumerate().map(|(b, v)| rise(b, *v)).sum() }); + prev = cur; + } + let frames_of = |secs: f32| ((secs * sr) as usize / hop).max(1); + let (near, before, after, gap) = (frames_of(0.03), frames_of(0.1), frames_of(0.07), frames_of(0.03)); + let delta = 2.5 * 2f32.powf((0.5 - sens) * 4.0); + let mut out: Vec = Vec::new(); + let mut last: Option = None; + for (m, &cur) in flux.iter().enumerate() { + let span = |a: usize, b: usize| flux.get(m.saturating_sub(a)..(m + b + 1).min(frames)).unwrap_or(&[]); + let peak = span(near, near).iter().copied().fold(0.0, f32::max); + let around = span(before, after); + let mean = around.iter().sum::() / around.len().max(1) as f32; + if cur <= 0.0 || cur < peak || cur < mean + delta || last.is_some_and(|l| m - l < gap) { + continue; + } + let (at, contrast) = attack_start(&ch, m * hop, n, sr); + // A fast release smears energy between a note's partials, which looks like a rise. Its + // attack edge is weak and the energy doesn't grow, unlike a real onset's. + let was = energy.get(m.saturating_sub(1)).copied().unwrap_or(0.0); + let grows = energy.get(m..(m + 2).min(frames)).is_some_and(|e| e.iter().any(|v| *v > was)); + if contrast < std::f64::consts::LN_2 && !grows { + continue; + } + last = Some(m); + if out.last().is_none_or(|&p| at > p) { + out.push(at); } } out } +/// Triangular filters 24 to the octave from 30 Hz to 20 kHz over the bins of an `n`-point FFT, without +/// duplicates where bins are coarser than that; each filter's weights sum to 1. +fn log_bands(n: usize, sr: f32) -> Vec> { + let bin_hz = sr / n as f32; + let top = 20_000.0f32.min(sr * 0.45); + let mut centres: Vec = Vec::new(); + let mut f = 30.0f32; + while f <= top { + let k = (f / bin_hz).round() as usize; + if centres.last().is_none_or(|&c| k > c) { + centres.push(k); + } + f *= 2f32.powf(1.0 / 24.0); + } + centres + .windows(3) + .filter_map(|w| { + let (&lo, &mid, &hi) = (w.first()?, w.get(1)?, w.get(2)?); + let tri: Vec<(usize, f32)> = (lo..=hi) + .map(|k| (k, if k <= mid { (k - lo) as f32 / (mid - lo).max(1) as f32 } else { (hi - k) as f32 / (hi - mid).max(1) as f32 })) + .filter(|(_, v)| *v > 0.0) + .collect(); + let sum: f32 = tri.iter().map(|(_, v)| v).sum(); + (sum > 0.0).then(|| tri.into_iter().map(|(k, v)| (k, v / sum)).collect()) + }) + .collect() +} + +/// Where the attack found in the frame starting at `frame` begins. The log spectrum jumps as soon as the +/// window reaches an attack, so it can sit up to half a window past the frame. Find the step edge of the +/// high-passed energy, the sample where the next 3 ms gain most over the 6 ms before, then go back to +/// where the energy first crosses the midpoint (in dB) between the two. +fn attack_start(ch: &[Vec], frame: usize, n: usize, sr: f32) -> (usize, f64) { + let (before, after, short) = (((sr * 0.006) as usize).max(1), ((sr * 0.003) as usize).max(1), ((sr * 0.0005) as usize).max(1)); + let (from, to) = (frame.saturating_sub(before), frame + n + n / 2 + after); + let mut cum = vec![0.0f64; to - from + 1]; + for (j, i) in (from..to).enumerate() { + let e: f32 = ch + .iter() + .map(|x| { + let d = x.get(i).copied().unwrap_or(0.0) - i.checked_sub(1).and_then(|p| x.get(p)).copied().unwrap_or(0.0); + d * d + }) + .sum(); + let sum = cum.get(j).copied().unwrap_or(0.0) + f64::from(e); + if let Some(c) = cum.get_mut(j + 1) { + *c = sum; + } + } + let at = |i: usize| cum.get(i.saturating_sub(from)).copied().unwrap_or(0.0); + let mean = |a: usize, b: usize| (at(b) - at(a)) / b.saturating_sub(a).max(1) as f64; + let contrast = |k: usize| ((mean(k, k + after) + 1e-12) / (mean(k.saturating_sub(before), k) + 1e-12)).ln(); + let k = (frame..frame + n + n / 2).max_by(|&a, &b| contrast(a).total_cmp(&contrast(b))).unwrap_or(frame); + let level = (mean(k.saturating_sub(before), k) * mean(k, k + after)).sqrt(); + ((k.saturating_sub(after)..=k).find(|&i| mean(i, i + short) > level).unwrap_or(k), contrast(k)) +} + /// Strip Silence: ranges `(start, end)` (end exclusive) of material whose peak across channels /// exceeds `threshold_db`. Silent gaps shorter than `min_len` samples do not split a region. /// Each region is padded by `pad_before`/`pad_after` (clamped to the buffer) and overlapping diff --git a/crates/dsp/tests/offline_meters.rs b/crates/dsp/tests/offline_meters.rs index 305ef17..1f9bd9b 100644 --- a/crates/dsp/tests/offline_meters.rs +++ b/crates/dsp/tests/offline_meters.rs @@ -254,7 +254,7 @@ fn transients_detected_near_onsets() { let mut x = vec![0.0f32; 48_000 * 2]; let onsets = [9_600usize, 33_600, 60_000, 81_000]; for &o in &onsets { - for i in 0..6000 { + for i in 0..15000 { x[o + i] += 0.8 * (-(i as f32) / 1500.0).exp() * (TAU * 180.0 * i as f32 / sr).sin(); } } diff --git a/crates/dsp/tests/transients.rs b/crates/dsp/tests/transients.rs new file mode 100644 index 0000000..893a048 --- /dev/null +++ b/crates/dsp/tests/transients.rs @@ -0,0 +1,149 @@ +use soundcraft_dsp::offline::detect_transients; +use std::f32::consts::TAU; + +const SR: f32 = 48_000.0; + +fn at(secs: f32) -> usize { + (secs * SR) as usize +} + +fn mix(len: f32, parts: &[(f32, Vec)]) -> Vec { + let mut x = vec![0.0; at(len)]; + for (start, sig) in parts { + let fade = at(0.03).min(sig.len()); + for (i, v) in sig.iter().enumerate() { + let left = sig.len() - i; + let g = if left < fade { 0.5 - 0.5 * (std::f32::consts::PI * left as f32 / fade as f32).cos() } else { 1.0 }; + if let Some(o) = x.get_mut(at(*start) + i) { + *o += v * g; + } + } + } + x +} + +fn noise(seed: &mut u32) -> f32 { + *seed ^= *seed << 13; + *seed ^= *seed >> 17; + *seed ^= *seed << 5; + *seed as f32 / u32::MAX as f32 * 2.0 - 1.0 +} + +fn hit(level: f32, decay: f32, seed: &mut u32) -> Vec { + let mut prev = 0.0; + (0..at(decay * 6.0)) + .map(|i| { + let w = noise(seed); + let hp = w - prev; + prev = w; + level * hp * (-(i as f32) / SR / decay).exp() + }) + .collect() +} + +fn note(f: f32, level: f32, secs: f32) -> Vec { + (0..at(secs)) + .map(|i| { + let t = i as f32 / SR; + let s: f32 = (1..=8).map(|k| (TAU * f * k as f32 * t).sin() * (-t * (0.6 + 0.4 * k as f32)).exp() / k as f32).sum(); + level * s * (t / 0.002).min(1.0) + }) + .collect() +} + +fn drums(len: f32) -> (Vec, Vec) { + let mut seed = 7; + let mut parts = Vec::new(); + let mut onsets = Vec::new(); + let mut t = 0.3; + while t < len - 0.5 { + parts.push((t, hit(0.9, 0.04, &mut seed))); + onsets.push(at(t)); + for k in 1..4 { + parts.push((t + 0.14 * k as f32, hit(0.09, 0.03, &mut seed))); + onsets.push(at(t + 0.14 * k as f32)); + } + t += 0.56; + } + (mix(len, &parts), onsets) +} + +fn near(found: &[usize], onsets: &[usize], tolerance: usize) -> bool { + found.len() == onsets.len() && found.iter().zip(onsets).all(|(f, o)| f.abs_diff(*o) <= tolerance) +} + +#[test] +fn notes_starting_over_a_ringing_note_are_placed_on_their_attacks() { + let starts = [0.3, 0.72, 1.14, 1.56, 1.98, 2.4]; + let freqs = [196.0, 261.6, 196.0, 329.6, 246.9, 392.0]; + let parts: Vec<(f32, Vec)> = + starts.iter().zip(freqs).enumerate().map(|(k, (s, f))| (*s, note(f, 0.2 + 0.1 * (k % 3) as f32, 1.1))).collect(); + let found = detect_transients(&[mix(4.0, &parts)], SR, 0.5); + let onsets: Vec = starts.iter().map(|s| at(*s)).collect(); + assert!(near(&found, &onsets, at(0.003)), "{found:?} vs {onsets:?}"); +} + +#[test] +fn quieter_hits_between_loud_ones_count() { + let (x, onsets) = drums(6.0); + let found = detect_transients(&[x], SR, 0.5); + assert!(near(&found, &onsets, at(0.001)), "{found:?} vs {onsets:?}"); +} + +#[test] +fn a_held_chord_with_vibrato_and_chorus_has_no_transients() { + let mut ph = [0.0f32; 6]; + let x: Vec = (0..at(6.0)) + .map(|i| { + let t = i as f32 / SR; + let mut s = 0.0; + for (j, f) in [220.0f32, 277.2, 329.6].iter().enumerate() { + for (d, sign) in [(0, 1.0f32), (1, -1.0)] { + let p = &mut ph[j * 2 + d]; + *p += TAU * f * (1.0 + sign * 0.004 * (TAU * (0.3 + 0.07 * j as f32) * t).sin()) / SR; + s += (1..=6).map(|k| (*p * k as f32).sin() / k as f32).sum::(); + } + } + 0.04 * s * (t / 0.4).min(1.0) * (1.0 + 0.3 * (TAU * 4.0 * t).sin()) + }) + .collect(); + let found: Vec = detect_transients(&[x], SR, 0.5).into_iter().filter(|p| *p > at(0.5)).collect(); + assert!(found.is_empty(), "{found:?}"); +} + +#[test] +fn transients_do_not_depend_on_where_the_analysed_audio_starts() { + let (mut x, _) = drums(8.0); + let from = at(3.3); + x.iter_mut().skip(from - at(0.4)).for_each(|v| *v *= 0.05); + let whole = detect_transients(&[x.clone()], SR, 0.5); + let part: Vec = detect_transients(&[x.get(from..).unwrap_or_default().to_vec()], SR, 0.5).into_iter().map(|p| p + from).collect(); + let settled = |v: &[usize]| v.iter().copied().filter(|p| *p > from + at(0.25)).collect::>(); + assert!(!settled(&whole).is_empty()); + assert_eq!(settled(&whole), settled(&part), "{whole:?} vs {part:?}"); +} + +#[test] +fn a_fast_release_is_not_a_transient() { + let starts: Vec = (0..8).map(|k| 0.3 + k as f32 * 0.5).collect(); + let note: Vec = (0..at(0.4)) + .map(|i| { + let t = i as f32 / SR; + let level = if t < 0.15 { (t / 0.01).min(1.0) } else { 0.001f32.powf((t - 0.15) / 0.04) }; + let saw: f32 = (1..=20).map(|k| (TAU * 659.0 * k as f32 * t).sin() / k as f32).sum(); + 0.2 * level * saw + }) + .collect(); + let parts: Vec<(f32, Vec)> = starts.iter().map(|s| (*s, note.clone())).collect(); + let found = detect_transients(&[mix(4.5, &parts)], SR, 0.5); + let onsets: Vec = starts.iter().map(|s| at(*s)).collect(); + assert!(near(&found, &onsets, at(0.002)), "{found:?} vs {onsets:?}"); +} + +#[test] +fn channels_that_cancel_when_mixed_still_show_their_transients() { + let (left, onsets) = drums(4.0); + let right: Vec = left.iter().map(|v| -v).collect(); + let found = detect_transients(&[left, right], SR, 0.5); + assert!(near(&found, &onsets, at(0.001)), "{found:?} vs {onsets:?}"); +}