Merge pull request #416 from grovebit/transient-detection
Transient detection: find quiet and legato onsets, place them on the attack
This commit is contained in:
Brandon Thomas
GitHub
commit
fe573bb31f
|
|
@ -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<f32>], sample_rate: f32, sensitivity: f32) -> Vec<usize> {
|
||||
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<Vec<f32>> =
|
||||
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<f32> = (0..len).map(|i| ch.iter().map(|c| c.get(i).copied().filter(|v| v.is_finite()).unwrap_or(0.0)).sum::<f32>() / nch).collect();
|
||||
let mut planner = FftPlanner::<f32>::new();
|
||||
let fft = planner.plan_fft_forward(n);
|
||||
let bands = log_bands(n, sr);
|
||||
let fft = FftPlanner::<f32>::new().plan_fft_forward(n);
|
||||
let window: Vec<f32> = (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<f32> = Vec::new();
|
||||
let mut flux = Vec::with_capacity(frames);
|
||||
let mut energy: Vec<f32> = 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<usize> = 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::<f32>() / 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<f32> = bands
|
||||
.iter()
|
||||
.map(|band| (1.0 + gain * band.iter().map(|&(k, w)| power.get(k).map_or(0.0, |p| p.sqrt()) * w).sum::<f32>()).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<usize> = Vec::new();
|
||||
let mut last: Option<usize> = 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::<f32>() / 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<Vec<(usize, f32)>> {
|
||||
let bin_hz = sr / n as f32;
|
||||
let top = 20_000.0f32.min(sr * 0.45);
|
||||
let mut centres: Vec<usize> = 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<f32>], 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
|
||||
|
|
|
|||
|
|
@ -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();
|
||||
}
|
||||
}
|
||||
|
|
|
|||
|
|
@ -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<f32>)]) -> Vec<f32> {
|
||||
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<f32> {
|
||||
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<f32> {
|
||||
(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<f32>, Vec<usize>) {
|
||||
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<f32>)> =
|
||||
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<usize> = 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<f32> = (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::<f32>();
|
||||
}
|
||||
}
|
||||
0.04 * s * (t / 0.4).min(1.0) * (1.0 + 0.3 * (TAU * 4.0 * t).sin())
|
||||
})
|
||||
.collect();
|
||||
let found: Vec<usize> = 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<usize> = 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::<Vec<_>>();
|
||||
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<f32> = (0..8).map(|k| 0.3 + k as f32 * 0.5).collect();
|
||||
let note: Vec<f32> = (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<f32>)> = starts.iter().map(|s| (*s, note.clone())).collect();
|
||||
let found = detect_transients(&[mix(4.5, &parts)], SR, 0.5);
|
||||
let onsets: Vec<usize> = 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<f32> = 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:?}");
|
||||
}
|
||||
Loading…
Reference in New Issue