diff --git a/Cargo.lock b/Cargo.lock index 3ea28e1..6394e4d 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -54,7 +54,7 @@ dependencies = [ [[package]] name = "audio-pattern-detector" -version = "0.4.2" +version = "0.4.3" dependencies = [ "audio-pattern-detector-core", "clap", @@ -64,9 +64,10 @@ dependencies = [ [[package]] name = "audio-pattern-detector-core" -version = "0.4.2" +version = "0.4.3" dependencies = [ "base64", + "fft-correlation", "pyo3", "realfft", "rustfft", @@ -168,6 +169,14 @@ version = "2.5.0" source = "registry+https://github.com/rust-lang/crates.io-index" checksum = "da7c62ceae207dd37ea5b845da6a0696c799f85e97da1ab5b7910be3c1c80223" +[[package]] +name = "fft-correlation" +version = "0.0.6" +source = "git+https://github.com/andrewtheguy/fft-correlation?tag=v0.0.6#7d31bffd13e0e86361251f190fdd6f64e8b1441d" +dependencies = [ + "realfft", +] + [[package]] name = "getrandom" version = "0.4.3" diff --git a/Cargo.toml b/Cargo.toml index b8bd3db..9be667a 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -3,7 +3,7 @@ resolver = "2" members = ["crates/core", "crates/cli"] [workspace.package] -version = "0.4.2" +version = "0.4.3" edition = "2021" authors = ["Andrew Chen "] diff --git a/README.md b/README.md index 4b99e09..e674748 100644 --- a/README.md +++ b/README.md @@ -25,7 +25,7 @@ cargo install --git https://github.com/andrewtheguy/audio_pattern_detector.git - `ghcr.io/andrewtheguy/audio_pattern_detector:` is published with each release (`v`, plus `latest` for releases from `main`). It contains only the binary at `/usr/local/bin/audio-pattern-detector` (no ffmpeg), so the usual use is copying it into another image: ```dockerfile -COPY --from=ghcr.io/andrewtheguy/audio_pattern_detector:v0.4.2 /usr/local/bin/audio-pattern-detector /usr/local/bin/audio-pattern-detector +COPY --from=ghcr.io/andrewtheguy/audio_pattern_detector:v0.4.3 /usr/local/bin/audio-pattern-detector /usr/local/bin/audio-pattern-detector ``` The release workflow builds the image for amd64 and arm64 from the binaries it has already compiled (`runtime-prebuilt` target). To build it locally from source for the host architecture: `docker build --target runtime -t audio-pattern-detector .` diff --git a/crates/core/Cargo.toml b/crates/core/Cargo.toml index 533c373..e33263f 100644 --- a/crates/core/Cargo.toml +++ b/crates/core/Cargo.toml @@ -17,6 +17,7 @@ python = ["dep:pyo3", "pyo3/extension-module"] [dependencies] base64 = "0.23" +fft-correlation = { git = "https://github.com/andrewtheguy/fft-correlation", tag = "v0.0.6" } pyo3 = { version = "0.28", optional = true, features = ["abi3-py39"] } realfft = "3.5" rustfft = "6.4" diff --git a/crates/core/src/detector.rs b/crates/core/src/detector.rs index 37c579c..111fec6 100644 --- a/crates/core/src/detector.rs +++ b/crates/core/src/detector.rs @@ -3,14 +3,15 @@ use std::path::PathBuf; use serde::ser::SerializeMap; use serde::{Serialize, Serializer}; +use fft_correlation::{fft_correlate_1d, CorrelationTemplate, CorrelationWorkspace, Mode}; use serde_json::json; use crate::audio_clip::{ validate_sample_rate, AudioClip, MarkerToneThresholds, Strategy, DEFAULT_TARGET_SAMPLE_RATE, }; use crate::dsp::{ - fft_correlate_full, find_peaks_1d, integrated_loudness, loudness_normalize, pearson_correlation_1d, - resample_preserve_maxima_1d, CorrelationTemplate, CorrelationWorkspace, FindPeaksOptions, + find_peaks_1d, integrated_loudness, loudness_normalize, pearson_correlation_1d, resample_preserve_maxima_1d, + FindPeaksOptions, }; use crate::error::{Error, Result}; use crate::stream::AudioStream; @@ -192,11 +193,11 @@ fn absolute_in_place(values: &mut [f32]) -> f32 { } /// Absolute self-correlation normalized to a peak of 1, plus the original peak. -fn clip_correlation(clip: &[f32]) -> (Vec, f32) { - let mut correlation = fft_correlate_full(clip, clip); +fn clip_correlation(clip: &[f32]) -> Result<(Vec, f32)> { + let mut correlation = fft_correlate_1d(clip, clip, Mode::Full)?; let absolute_max = absolute_in_place(&mut correlation); correlation.iter_mut().for_each(|v| *v /= absolute_max); - (correlation, absolute_max) + Ok((correlation, absolute_max)) } /// Clips whose sections are identical (same lookback), so one @@ -321,7 +322,7 @@ impl AudioPatternDetector { let clips: Vec = audio_clips .into_iter() .map(|audio_clip| Self::prepare_clip(audio_clip, sr, debug_mode)) - .collect(); + .collect::>()?; let section_groups = section_groups(&clips); Ok(Self { @@ -337,7 +338,7 @@ impl AudioPatternDetector { } /// Pre-compute everything about a clip that doesn't depend on the audio stream. - fn prepare_clip(audio_clip: AudioClip, sr: u32, debug_mode: bool) -> ClipData { + fn prepare_clip(audio_clip: AudioClip, sr: u32, debug_mode: bool) -> Result { let clip_seconds = audio_clip.clip_length_seconds(); let sliding_window = clip_seconds.ceil() as u32; if sliding_window as f64 != clip_seconds { @@ -349,7 +350,7 @@ impl AudioPatternDetector { let mut clip = audio_clip.audio; normalize_loudness(&mut clip, sr); - let (correlation_clip, correlation_clip_absolute_max) = clip_correlation(&clip); + let (correlation_clip, correlation_clip_absolute_max) = clip_correlation(&clip)?; if debug_mode { eprintln!("clip_length {} {}", audio_clip.name, clip.len()); @@ -377,7 +378,7 @@ impl AudioPatternDetector { .map(|&window| downsample_window(&correlation_clip, window)) .collect(); - ClipData { + Ok(ClipData { name: audio_clip.name, clip, sliding_window, @@ -385,7 +386,7 @@ impl AudioPatternDetector { correlation_clip_absolute_max, tone, pearson_windows, - } + }) } /// Seconds per chunk actually used (after auto-computation). @@ -571,7 +572,7 @@ impl AudioPatternDetector { let lookback_seconds = lookback_samples as f64 / sr as f64; let RunBuffers { workspace, templates, audio_section, correlation } = buffers; - workspace.correlate(&mut templates[clip_index], correlation); + workspace.correlate(&mut templates[clip_index], Mode::Full, correlation)?; let peaks = self.correlation_method(clip_data, audio_section, correlation, index)?; let section_start = index as i64 * self.seconds_per_chunk as i64 * sr as i64 - lookback_samples as i64; @@ -979,7 +980,7 @@ mod tests { fn test_clip_correlation_is_normalized_with_original_peak() { let clip = [0.5_f32, -1.0, 0.25, 0.75]; let energy: f32 = clip.iter().map(|v| v * v).sum(); - let (correlation, absolute_max) = clip_correlation(&clip); + let (correlation, absolute_max) = clip_correlation(&clip).unwrap(); assert_eq!(correlation.len(), 2 * clip.len() - 1); assert!((absolute_max - energy).abs() < 1e-5, "{absolute_max} != {energy}"); assert_eq!(correlation[clip.len() - 1], 1.0); diff --git a/crates/core/src/dsp/correlate.rs b/crates/core/src/dsp/correlate.rs deleted file mode 100644 index 800ac99..0000000 --- a/crates/core/src/dsp/correlate.rs +++ /dev/null @@ -1,428 +0,0 @@ -// ── FFT cross-correlation ──────────────────────────────────────────── - -use realfft::num_complex::Complex; -use realfft::{ComplexToReal, RealFftPlanner, RealToComplex}; -use std::cell::RefCell; -use std::collections::VecDeque; -use std::sync::Arc; - -const FFT_PLAN_CACHE_CAPACITY: usize = 8; - -type FftPlans = (Arc>, Arc>); - -// Bounded per-thread cache of FFT plans keyed by size, most recently used first. -thread_local! { - static FFT_PLAN_CACHE: RefCell> = - RefCell::new(VecDeque::with_capacity(FFT_PLAN_CACHE_CAPACITY)); -} - -fn get_fft_plans(fft_size: usize) -> FftPlans { - FFT_PLAN_CACHE.with(|cache_cell| { - let mut cache = cache_cell.borrow_mut(); - if let Some(index) = cache.iter().position(|(size, _)| *size == fft_size) { - let entry = cache.remove(index).expect("cache entry exists at located index"); - let plans = entry.1.clone(); - cache.push_front(entry); - return plans; - } - - // Use a fresh planner per cache miss so planner-internal maps do not grow unbounded. - let mut planner = RealFftPlanner::::new(); - let plans: FftPlans = ( - planner.plan_fft_forward(fft_size), - planner.plan_fft_inverse(fft_size), - ); - if cache.len() == FFT_PLAN_CACHE_CAPACITY { - cache.pop_back(); - } - cache.push_front((fft_size, plans.clone())); - plans - }) -} - -/// Number of FFT sizes a [`CorrelationTemplate`] keeps spectra for. A run -/// sees at most three: the first chunk (no lookback), full chunks and the -/// final short chunk. -const TEMPLATE_SPECTRUM_CACHE_CAPACITY: usize = 4; - -/// A correlation template with its time-reversed spectrum cached per FFT -/// size, so repeated correlations against new signals only transform the -/// signal. -pub struct CorrelationTemplate { - reversed: Vec, - spectra: Vec<(usize, Vec>)>, -} - -impl CorrelationTemplate { - pub fn new(template: &[f32]) -> Self { - Self { - // Time-reversing the template turns convolution into correlation. - reversed: template.iter().rev().copied().collect(), - spectra: Vec::new(), - } - } - - pub fn len(&self) -> usize { - self.reversed.len() - } - - pub fn is_empty(&self) -> bool { - self.reversed.is_empty() - } - - /// Spectrum of the zero-padded, reversed template at `fft_size`. - fn spectrum(&mut self, fft_size: usize) -> &[Complex] { - if let Some(index) = self.spectra.iter().position(|(size, _)| *size == fft_size) { - return &self.spectra[index].1; - } - assert!(fft_size >= self.reversed.len(), "FFT size must cover the template"); - - let (r2c, _) = get_fft_plans(fft_size); - let mut padded = vec![0.0_f32; fft_size]; - padded[..self.reversed.len()].copy_from_slice(&self.reversed); - let mut spectrum = r2c.make_output_vec(); - let mut scratch = r2c.make_scratch_vec(); - r2c.process_with_scratch(&mut padded, &mut spectrum, &mut scratch) - .expect("forward FFT buffers are sized by the plan"); - - if self.spectra.len() == TEMPLATE_SPECTRUM_CACHE_CAPACITY { - self.spectra.remove(0); - } - self.spectra.push((fft_size, spectrum)); - &self.spectra.last().expect("spectrum was just pushed").1 - } -} - -/// Reusable buffers for FFT cross-correlation of one signal against any -/// number of templates: the signal is set by [`load_signal`] and transformed -/// once per FFT size, so each [`correlate`] call with a template needing the -/// same size costs a spectrum product and one inverse FFT. -/// -/// The FFT size of a correlation depends only on the signal and that -/// template, never on the other templates, so the result is bit-identical -/// to correlating the pair on its own. -/// -/// [`load_signal`]: CorrelationWorkspace::load_signal -/// [`correlate`]: CorrelationWorkspace::correlate -#[derive(Default)] -pub struct CorrelationWorkspace { - signal: Vec, - /// FFT size the buffers are allocated for (0 before the first transform). - fft_size: usize, - /// Whether `signal_spectrum` is the transform of `signal` at `fft_size`. - spectrum_loaded: bool, - plans: Option, - padded_signal: Vec, - signal_spectrum: Vec>, - product: Vec>, - forward_scratch: Vec>, - inverse_scratch: Vec>, - time_domain: Vec, -} - -impl CorrelationWorkspace { - pub fn new() -> Self { - Self::default() - } - - /// Set the signal for the following [`correlate`](Self::correlate) calls. - pub fn load_signal(&mut self, signal: &[f32]) { - self.signal.clear(); - self.signal.extend_from_slice(signal); - self.spectrum_loaded = false; - } - - /// Transform the loaded signal at `fft_size` unless that is already done. - fn transform_signal(&mut self, fft_size: usize) { - if self.spectrum_loaded && self.fft_size == fft_size { - return; - } - if self.fft_size != fft_size || self.plans.is_none() { - let plans = get_fft_plans(fft_size); - let (r2c, c2r) = &plans; - self.signal_spectrum = r2c.make_output_vec(); - self.product = r2c.make_output_vec(); - self.forward_scratch = r2c.make_scratch_vec(); - self.inverse_scratch = c2r.make_scratch_vec(); - self.time_domain = c2r.make_output_vec(); - self.plans = Some(plans); - self.fft_size = fft_size; - } - - self.padded_signal.clear(); - self.padded_signal.resize(fft_size, 0.0); - self.padded_signal[..self.signal.len()].copy_from_slice(&self.signal); - - let (r2c, _) = self.plans.as_ref().expect("plans were just set"); - r2c.process_with_scratch(&mut self.padded_signal, &mut self.signal_spectrum, &mut self.forward_scratch) - .expect("forward FFT buffers are sized by the plan"); - self.spectrum_loaded = true; - } - - /// Full cross-correlation of the loaded signal with `template`, written - /// to `output` (length `signal.len() + template.len() - 1`). - /// - /// Matches `scipy.signal.correlate(signal, template, mode="full")`: - /// output index `k` is the lag where `template[template.len() - 1]` - /// aligns with `signal[k]`. `output` is empty if either is empty. - pub fn correlate(&mut self, template: &mut CorrelationTemplate, output: &mut Vec) { - output.clear(); - if self.signal.is_empty() || template.is_empty() { - return; - } - - // Covers the full correlation without wrap-around. - let output_len = self.signal.len() + template.len() - 1; - self.transform_signal(output_len.next_power_of_two()); - - let template_spectrum = template.spectrum(self.fft_size); - for ((p, s), t) in self.product.iter_mut().zip(&self.signal_spectrum).zip(template_spectrum) { - *p = *s * *t; - } - - let (_, c2r) = self.plans.as_ref().expect("the signal was just transformed"); - c2r.process_with_scratch(&mut self.product, &mut self.time_domain, &mut self.inverse_scratch) - .expect("inverse FFT buffers are sized by the plan"); - - let normalization = 1.0 / self.fft_size as f32; - output.extend(self.time_domain[..output_len].iter().map(|x| x * normalization)); - } -} - -/// Cross-correlate two real 1-D signals using FFT, returning the full -/// correlation (length `signal.len() + template.len() - 1`). -/// -/// One-shot convenience over [`CorrelationWorkspace`]; see -/// [`CorrelationWorkspace::correlate`] for the semantics. -pub fn fft_correlate_full(signal: &[f32], template: &[f32]) -> Vec { - let mut workspace = CorrelationWorkspace::new(); - let mut template = CorrelationTemplate::new(template); - let mut output = Vec::new(); - workspace.load_signal(signal); - workspace.correlate(&mut template, &mut output); - output -} - -#[cfg(test)] -mod tests { - use super::*; - - // Naive correlation used for correctness checks. - fn naive_full_correlation(signal: &[f32], template: &[f32]) -> Vec { - let output_len = signal.len() + template.len() - 1; - let mut result = vec![0.0; output_len]; - for (lag, out) in result.iter_mut().enumerate() { - for (i, &t) in template.iter().enumerate() { - let signal_idx = lag as isize - (template.len() as isize - 1) + i as isize; - if (0..signal.len() as isize).contains(&signal_idx) { - *out += signal[signal_idx as usize] * t; - } - } - } - result - } - - #[test] - fn test_full_length() { - let result = fft_correlate_full(&[1.0, 2.0, 3.0, 4.0, 5.0], &[1.0, 0.0, 0.0]); - assert_eq!(result.len(), 7); - assert_eq!(fft_correlate_full(&[1.0; 100], &[1.0; 10]).len(), 109); - } - - #[test] - fn test_empty_inputs() { - assert_eq!(fft_correlate_full(&[], &[1.0]), Vec::::new()); - assert_eq!(fft_correlate_full(&[1.0], &[]), Vec::::new()); - assert_eq!(fft_correlate_full(&[1.0, 2.0], &[]), Vec::::new()); - assert_eq!(fft_correlate_full(&[1.0, 2.0, 3.0], &[]), Vec::::new()); - } - - #[test] - fn test_matches_naive() { - let signal: Vec = (0..37).map(|i| ((i * 7 % 11) as f32 - 5.0) / 5.0).collect(); - let template: Vec = (0..9).map(|i| ((i * 5 % 7) as f32 - 3.0) / 3.0).collect(); - let fft = fft_correlate_full(&signal, &template); - let naive = naive_full_correlation(&signal, &template); - assert_eq!(fft.len(), naive.len()); - for (a, b) in fft.iter().zip(naive.iter()) { - assert!((a - b).abs() < 1e-4, "{a} != {b}"); - } - } - - #[test] - fn test_workspace_reuse_matches_one_shot() { - let signal_a: Vec = (0..50).map(|i| ((i * 13 % 17) as f32 - 8.0) / 8.0).collect(); - let signal_b: Vec = (0..70).map(|i| ((i * 3 % 19) as f32 - 9.0) / 9.0).collect(); - let template_a: Vec = (0..7).map(|i| ((i * 5 % 7) as f32 - 3.0) / 3.0).collect(); - let template_b: Vec = (0..12).map(|i| ((i * 11 % 13) as f32 - 6.0) / 6.0).collect(); - - let mut workspace = CorrelationWorkspace::new(); - let mut cached_a = CorrelationTemplate::new(&template_a); - let mut cached_b = CorrelationTemplate::new(&template_b); - let mut output = Vec::new(); - - // Two templates against one signal, then a different signal size. - for signal in [&signal_a, &signal_b, &signal_a] { - workspace.load_signal(signal); - for (cached, template) in [(&mut cached_a, &template_a), (&mut cached_b, &template_b)] { - workspace.correlate(cached, &mut output); - let naive = naive_full_correlation(signal, template); - assert_eq!(output.len(), naive.len()); - for (a, b) in output.iter().zip(naive.iter()) { - assert!((a - b).abs() < 1e-4, "{a} != {b}"); - } - } - } - assert_eq!(cached_a.spectra.len(), 2, "one cached spectrum per FFT size"); - } - - #[test] - fn test_workspace_empty_inputs() { - let mut workspace = CorrelationWorkspace::new(); - let mut template = CorrelationTemplate::new(&[1.0, 2.0]); - let mut output = vec![1.0]; - workspace.load_signal(&[]); - workspace.correlate(&mut template, &mut output); - assert!(output.is_empty()); - - let mut empty = CorrelationTemplate::new(&[]); - workspace.load_signal(&[1.0, 2.0, 3.0]); - workspace.correlate(&mut empty, &mut output); - assert!(output.is_empty()); - } - - fn assert_close(actual: &[f32], expected: &[f32]) { - assert_eq!(actual.len(), expected.len()); - for (a, b) in actual.iter().zip(expected) { - assert!((a - b).abs() < 1e-4, "{a} != {b}"); - } - } - - #[test] - fn test_template_len_and_is_empty() { - let template = CorrelationTemplate::new(&[1.0, 2.0, 3.0]); - assert_eq!(template.len(), 3); - assert!(!template.is_empty()); - assert_eq!(template.reversed, vec![3.0, 2.0, 1.0]); - - let empty = CorrelationTemplate::new(&[]); - assert_eq!(empty.len(), 0); - assert!(empty.is_empty()); - } - - #[test] - fn test_cached_spectrum_gives_identical_output() { - let signal: Vec = (0..90).map(|i| ((i * 29 % 31) as f32 - 15.0) / 15.0).collect(); - let template: Vec = (0..11).map(|i| ((i * 7 % 13) as f32 - 6.0) / 6.0).collect(); - - let mut workspace = CorrelationWorkspace::new(); - let mut cached = CorrelationTemplate::new(&template); - let mut first = Vec::new(); - let mut second = Vec::new(); - - workspace.load_signal(&signal); - workspace.correlate(&mut cached, &mut first); - assert_eq!(cached.spectra.len(), 1); - // Same signal again: the cached spectrum is used (no new entry) and - // the result is bit-identical to the first pass and to the one-shot. - workspace.load_signal(&signal); - workspace.correlate(&mut cached, &mut second); - assert_eq!(cached.spectra.len(), 1); - assert_eq!(first, second); - assert_eq!(first, fft_correlate_full(&signal, &template)); - } - - #[test] - fn test_template_spectrum_cache_evicts_oldest() { - let template: Vec = (0..5).map(|i| ((i * 3 % 5) as f32 - 2.0) / 2.0).collect(); - let mut cached = CorrelationTemplate::new(&template); - let mut workspace = CorrelationWorkspace::new(); - let mut output = Vec::new(); - - // Six signal lengths, each needing a different power-of-two FFT. - let lengths = [4, 12, 28, 60, 124, 252]; - let mut fft_sizes = Vec::new(); - for &len in &lengths { - let signal: Vec = (0..len).map(|i| ((i * 17 % 23) as f32 - 11.0) / 11.0).collect(); - workspace.load_signal(&signal); - workspace.correlate(&mut cached, &mut output); - assert_close(&output, &naive_full_correlation(&signal, &template)); - fft_sizes.push(workspace.fft_size); - assert!(cached.spectra.len() <= TEMPLATE_SPECTRUM_CACHE_CAPACITY); - } - assert_eq!(fft_sizes, vec![8, 16, 32, 64, 128, 256]); - - // The oldest two sizes were evicted, the newest four kept in order. - let kept: Vec = cached.spectra.iter().map(|(size, _)| *size).collect(); - assert_eq!(kept, vec![32, 64, 128, 256]); - - // An evicted size is recomputed and still correct. - let signal: Vec = (0..4).map(|i| i as f32 - 1.5).collect(); - workspace.load_signal(&signal); - workspace.correlate(&mut cached, &mut output); - assert_close(&output, &naive_full_correlation(&signal, &template)); - let kept: Vec = cached.spectra.iter().map(|(size, _)| *size).collect(); - assert_eq!(kept, vec![64, 128, 256, 8]); - } - - #[test] - fn test_fft_size_is_independent_of_other_templates() { - let signal: Vec = (0..40).map(|i| ((i * 19 % 29) as f32 - 14.0) / 14.0).collect(); - let short: Vec = (0..6).map(|i| ((i * 5 % 7) as f32 - 3.0) / 3.0).collect(); - let long: Vec = (0..30).map(|i| ((i * 11 % 17) as f32 - 8.0) / 8.0).collect(); - let expected_short = fft_correlate_full(&signal, &short); - let expected_long = fft_correlate_full(&signal, &long); - - let mut workspace = CorrelationWorkspace::new(); - let mut output = Vec::new(); - workspace.load_signal(&signal); - // Each template gets the FFT size it needs on its own, in any order, - // so the output is bit-identical to the one-shot correlation. - let mut fft_sizes = Vec::new(); - for (template, expected) in [(&long, &expected_long), (&short, &expected_short), (&long, &expected_long)] { - workspace.correlate(&mut CorrelationTemplate::new(template), &mut output); - fft_sizes.push(workspace.fft_size); - assert_eq!(&output, expected); - assert_close(&output, &naive_full_correlation(&signal, template)); - } - assert_eq!(fft_sizes, vec![128, 64, 128]); - } - - #[test] - fn test_template_longer_than_signal_matches_naive() { - let signal = [1.0_f32, -2.0, 0.5]; - let template: Vec = (0..9).map(|i| ((i * 4 % 9) as f32 - 4.0) / 4.0).collect(); - let result = fft_correlate_full(&signal, &template); - assert_eq!(result.len(), 11); - assert_close(&result, &naive_full_correlation(&signal, &template)); - } - - #[test] - fn test_exact_values_match_scipy() { - // scipy.signal.correlate([1, 2, 3, 4], [1, 0.5], mode="full") - // = [0.5, 2.0, 3.5, 5.0, 4.0] - let result = fft_correlate_full(&[1.0, 2.0, 3.0, 4.0], &[1.0, 0.5]); - assert_close(&result, &[0.5, 2.0, 3.5, 5.0, 4.0]); - // Single-sample template scales the signal. - let result = fft_correlate_full(&[1.0, -2.0, 3.0], &[2.0]); - assert_close(&result, &[2.0, -4.0, 6.0]); - } - - #[test] - fn test_output_buffer_is_replaced_not_appended() { - let mut workspace = CorrelationWorkspace::new(); - let mut template = CorrelationTemplate::new(&[1.0, 0.5]); - let mut output = vec![9.0; 20]; - workspace.load_signal(&[1.0, 2.0, 3.0, 4.0]); - workspace.correlate(&mut template, &mut output); - assert_close(&output, &[0.5, 2.0, 3.5, 5.0, 4.0]); - } - - #[test] - fn test_autocorrelation_peak_is_centered() { - let clip = [0.5_f32, -1.0, 0.25, 0.75]; - let corr = fft_correlate_full(&clip, &clip); - let energy: f32 = clip.iter().map(|v| v * v).sum(); - assert!((corr[clip.len() - 1] - energy).abs() < 1e-5); - } -} diff --git a/crates/core/src/dsp/mod.rs b/crates/core/src/dsp/mod.rs index 266d0e8..1af4b40 100644 --- a/crates/core/src/dsp/mod.rs +++ b/crates/core/src/dsp/mod.rs @@ -1,14 +1,13 @@ -//! Numerical routines used by the detector: FFT cross-correlation, BS.1770 -//! loudness, peak finding, resampling, Pearson correlation and real spectra. +//! Numerical routines used by the detector: BS.1770 loudness, peak finding, +//! resampling, Pearson correlation and real spectra. FFT cross-correlation +//! comes from the `fft-correlation` crate. -pub mod correlate; pub mod loudness; pub mod pearson; pub mod peaks; pub mod resample; pub mod spectrum; -pub use correlate::{fft_correlate_full, CorrelationTemplate, CorrelationWorkspace}; pub use loudness::{integrated_loudness, loudness_normalize}; pub use pearson::pearson_correlation_1d; pub use peaks::{find_peaks_1d, FindPeaksOptions}; diff --git a/crates/core/src/error.rs b/crates/core/src/error.rs index 1288fb2..96d6d1e 100644 --- a/crates/core/src/error.rs +++ b/crates/core/src/error.rs @@ -7,6 +7,8 @@ pub enum Error { Invalid(String), /// Underlying I/O failure. Io(std::io::Error), + /// FFT cross-correlation failure. + Correlation(fft_correlation::FftCorrelationError), } pub type Result = std::result::Result; @@ -22,6 +24,7 @@ impl fmt::Display for Error { match self { Error::Invalid(message) => f.write_str(message), Error::Io(e) => write!(f, "I/O error: {e}"), + Error::Correlation(e) => write!(f, "correlation error: {e}"), } } } @@ -31,6 +34,7 @@ impl std::error::Error for Error { match self { Error::Invalid(_) => None, Error::Io(e) => Some(e), + Error::Correlation(e) => Some(e), } } } @@ -40,3 +44,9 @@ impl From for Error { Error::Io(e) } } + +impl From for Error { + fn from(e: fft_correlation::FftCorrelationError) -> Self { + Error::Correlation(e) + } +} diff --git a/crates/core/src/python.rs b/crates/core/src/python.rs index b7b8413..71b3637 100644 --- a/crates/core/src/python.rs +++ b/crates/core/src/python.rs @@ -4,7 +4,7 @@ use std::io::{self, Read}; use std::path::PathBuf; use std::sync::Mutex; -use pyo3::exceptions::{PyTypeError, PyValueError}; +use pyo3::exceptions::{PyRuntimeError, PyTypeError, PyValueError}; use pyo3::prelude::*; use pyo3::types::{PyBytes, PyDict}; @@ -23,6 +23,7 @@ fn to_py_err(error: Error) -> PyErr { match error { Error::Invalid(message) => PyValueError::new_err(message), Error::Io(e) => e.into(), + Error::Correlation(e) => PyRuntimeError::new_err(e.to_string()), } } diff --git a/docs/development.md b/docs/development.md index 8dadc59..1e83ac2 100644 --- a/docs/development.md +++ b/docs/development.md @@ -49,9 +49,9 @@ The version is shared through `[workspace.package]` in the root `Cargo.toml`. | `crates/core/src/wav.rs` | WAV reading and writing | | `crates/core/src/ffmpeg.rs` | ffmpeg subprocess for non-WAV files | | `crates/core/src/python.rs` | Python bindings (PyO3), compiled only with the `python` feature; see [python.md](python.md) | -| `crates/core/src/dsp/` | FFT cross-correlation, BS.1770 loudness, peak finding, resampling, Pearson correlation, spectra | +| `crates/core/src/dsp/` | BS.1770 loudness, peak finding, resampling, Pearson correlation, spectra | -The numerical routines in `crates/core/src/dsp/` are implemented in the crate rather than pulled in as dependencies; they follow the semantics of the scipy/numpy/pyloudnorm functions the algorithm was originally tuned against (`scipy.signal.find_peaks`, `scipy.signal.resample`, `scipy.signal.correlate`, `numpy.hanning`, BS.1770 integrated loudness), so thresholds carry over unchanged. +FFT cross-correlation comes from the [`fft-correlation`](https://github.com/andrewtheguy/fft-correlation) crate; fixes and optimizations to it belong there. The numerical routines in `crates/core/src/dsp/` and that crate follow the semantics of the scipy/numpy/pyloudnorm functions the algorithm was originally tuned against (`scipy.signal.find_peaks`, `scipy.signal.resample`, `scipy.signal.correlate`, `numpy.hanning`, BS.1770 integrated loudness), so thresholds carry over unchanged. ## Debug output diff --git a/docs/pattern-matching.md b/docs/pattern-matching.md index cdaa48d..c9d2728 100644 --- a/docs/pattern-matching.md +++ b/docs/pattern-matching.md @@ -49,7 +49,7 @@ Raw audio stream (float32 PCM) Before processing audio, each clip is prepared once: 1. **Loudness normalization** - clip audio normalized to -16 dB LUFS. -2. **Self-correlation** - FFT cross-correlation of the clip with itself (`fft_correlate_full(clip, clip)`), producing a reference correlation curve. The absolute max is stored for normalization later. +2. **Self-correlation** - FFT cross-correlation of the clip with itself (`fft_correlate_1d(clip, clip, Mode::Full)` from the `fft-correlation` crate), producing a reference correlation curve. The absolute max is stored for normalization later. 3. **Tone-strategy setup** - for clips whose strategy is `marker_tone`, the dominant frequency is recorded as clip metadata. All other clips have no tone verifier. This produces a `ClipData` struct per clip containing the normalized audio, clip name, sliding window, self-correlation curve, its absolute max, the downsampled Pearson windows of that curve, and the tone verifier (set only for tone-strategy clips). @@ -62,15 +62,15 @@ The raw chunks themselves are not read with overlap. Instead, the detector build The same prepend is applied to every chunk after the first, including a final short chunk. -Each `audio_section` is loudness-normalized to -16 dB LUFS before correlation, so a clip's detections depend only on its own section, never on which other clips are loaded. Clips with the same `sliding_window` have identical sections and share one: it is normalized once and transformed once per FFT size, and each clip's cross-correlation then costs a spectrum product with the clip's cached template spectrum and one inverse FFT. The FFT size is the one the section and that clip need on their own (next power of two of `len(audio_section) + len(clip) - 1`), so clips in a group usually share the forward FFT but never change each other's results. +Each `audio_section` is loudness-normalized to -16 dB LUFS before correlation, so a clip's detections depend only on its own section, never on which other clips are loaded. Clips with the same `sliding_window` have identical sections and share one: it is normalized once and transformed once per FFT size, and each clip's cross-correlation then costs a spectrum product with the clip's cached template spectrum and one inverse FFT. The FFT size is the one the section and that clip need on their own (the smallest `m * 2^k` with `m` in 1, 3, 5, 9, 15 that covers `len(audio_section) + len(clip) - 1`, so at most 25% above it), so clips in a group usually share the forward FFT but never change each other's results. ## Step 1: Candidate Detection (FFT Cross-Correlation) This step always runs first for every clip type. Its job is to locate and center potential match positions; it does not decide whether a candidate is a true match — that is Step 2. -For each chunk, each sliding-window group's `audio_section` is loaded into a `CorrelationWorkspace` (`load_signal`). Then for each clip in the group: +For each chunk, each sliding-window group's `audio_section` is loaded into a `CorrelationWorkspace` (`load_signal`) from the `fft-correlation` crate. Then for each clip in the group: -1. Compute the full cross-correlation with `CorrelationWorkspace::correlate`, using the clip's `CorrelationTemplate` (its time-reversed spectrum, cached per FFT size for the run), and take the absolute value. The result is the same as `fft_correlate_full(audio_section, clip)`; the workspace only avoids repeating the section's forward FFT. +1. Compute the full cross-correlation with `CorrelationWorkspace::correlate`, using the clip's `CorrelationTemplate` (its time-reversed spectrum, cached per FFT size for the run), and take the absolute value. The result is the same as `fft_correlate_1d(audio_section, clip, Mode::Full)`; the workspace only avoids repeating the section's forward FFT. 2. Normalize by `max(self_correlation_max, cross_correlation_max)` so the correlation curve is in [0, 1]. 3. Run peak detection with `height >= 0.25` and `distance >= clip_length` (prevents duplicate detections within one clip duration). diff --git a/docs/python.md b/docs/python.md index fb01857..4351d65 100644 --- a/docs/python.md +++ b/docs/python.md @@ -9,7 +9,7 @@ Wheels (Linux x86_64 and arm64, macOS Apple Silicon, Windows x86_64) are attache ```toml # pyproject.toml (uv) [project] -dependencies = ["audio-pattern-detector==0.4.2"] +dependencies = ["audio-pattern-detector==0.4.3"] [[tool.uv.index]] name = "audio-pattern-detector" @@ -85,7 +85,7 @@ Both match methods return a `MatchResult`: ### Behaviour -- Invalid input (missing files, bad WAV data, bad pattern configs, a chunk size too small for a clip) raises `ValueError`; I/O failures raise `OSError`. +- Invalid input (missing files, bad WAV data, bad pattern configs, a chunk size too small for a clip) raises `ValueError`; I/O failures raise `OSError`; an FFT correlation failure raises `RuntimeError`. - An exception raised by `on_detected` or by `stream.read()` is re-raised from the match call. A stream match stops reading as soon as that happens; a file match finishes the file first. - The GIL is released while matching, so detectors can run in threads. A `Detector` is immutable and can be shared between threads. - Diagnostics go to the process's stderr, as with the CLI.