From 5c57743678f5317ab8d55bd1141f9afd3aa6ad04 Mon Sep 17 00:00:00 2001 From: Sami Perttu Date: Tue, 26 Dec 2023 23:24:47 +0200 Subject: [PATCH] Frequency domain resynthesis. --- CHANGES.md | 1 + README.md | 1 + src/feedback.rs | 6 +- src/hacker.rs | 33 +++++- src/hacker32.rs | 33 +++++- src/prelude.rs | 34 +++++- src/resynth.rs | 292 +++++++++++++++++++++++++++++++++++++++++++++--- src/signal.rs | 10 +- tests/basic.rs | 25 +++++ 9 files changed, 408 insertions(+), 27 deletions(-) diff --git a/CHANGES.md b/CHANGES.md index ac3bfb0..d2a5d33 100644 --- a/CHANGES.md +++ b/CHANGES.md @@ -5,6 +5,7 @@ - `AudioNode` now requires `Send` and `Sync`. - Feedback units `Feedback64` and `Feedback32`. - `Shape::Atan` was contributed by Serdnad. +- New opcode `resynth` for frequency domain resynthesis. ### Version 0.15 diff --git a/README.md b/README.md index 6fb5e7e..dd95eab 100644 --- a/README.md +++ b/README.md @@ -903,6 +903,7 @@ The type parameters in the table refer to the hacker preludes. | `resample(node)` | 1 (speed) | `node` | Resample generator `node` using cubic interpolation at speed obtained from the input, where 1 is the original speed. | | `resonator()` | 3 (audio, frequency, bandwidth) | 1 | Constant-gain bandpass resonator (2nd order). | | `resonator_hz(f, bw)` | 1 | 1 | Constant-gain bandpass resonator (2nd order) with center frequency `f` Hz and bandwidth `bw` Hz. | +| `resynth::(w, f)` | `I` | `O` | Frequency domain resynthesis with window length `w` and processing function `f` | | `reverb_stereo(r, t)` | 2 | 2 | Stereo reverb with room size `r` meters (10 is average) and reverberation time `t` seconds. | | `reverse::()` | `N` | `N` | Reverse channel order, e.g., swap left and right channels. | | `rossler()` | 1 (frequency) | 1 | [Rössler dynamical system](https://en.wikipedia.org/wiki/R%C3%B6ssler_attractor) oscillator. | diff --git a/src/feedback.rs b/src/feedback.rs index d81d44b..aa074bc 100644 --- a/src/feedback.rs +++ b/src/feedback.rs @@ -158,7 +158,7 @@ where } fn route(&mut self, input: &SignalFrame, _frequency: f64) -> SignalFrame { - Routing::Arbitrary.propagate(input, self.outputs()) + Routing::Arbitrary(0.0).propagate(input, self.outputs()) } fn ping(&mut self, probe: bool, hash: AttoHash) -> AttoHash { @@ -264,7 +264,7 @@ where } fn route(&mut self, input: &SignalFrame, _frequency: f64) -> SignalFrame { - Routing::Arbitrary.propagate(input, self.outputs()) + Routing::Arbitrary(0.0).propagate(input, self.outputs()) } fn ping(&mut self, probe: bool, hash: AttoHash) -> AttoHash { @@ -437,7 +437,7 @@ impl AudioUnit48 for Feedback48 { } fn route(&mut self, input: &SignalFrame, _frequency: f64) -> SignalFrame { - Routing::Arbitrary.propagate(input, self.outputs()) + Routing::Arbitrary(0.0).propagate(input, self.outputs()) } fn get_id(&self) -> u64 { diff --git a/src/hacker.rs b/src/hacker.rs index 12edf78..55bd9b0 100644 --- a/src/hacker.rs +++ b/src/hacker.rs @@ -1051,7 +1051,7 @@ where O: ConstantFrame, O::Size: Size, { - An(Map::new(f, Routing::Arbitrary)) + An(Map::new(f, Routing::Arbitrary(0.0))) } /// Keeps a signal zero centered. @@ -2171,3 +2171,34 @@ pub fn snoop(capacity: usize) -> (Snoop, An>) { let (snoop, backend) = Snoop::new(capacity); (snoop, An(backend)) } + +/// Frequency domain resynthesizer. The number of inputs is `I` and the number of outputs is `O`. +/// The window length (in samples) must be a power of two and at least four. +/// Processes windows of input samples transformed into the frequency domain. +/// The processing function is issued arguments (time, window). +/// It processes frequency domain inputs into frequency domain outputs. +/// The latency in samples is equal to window length. +/// If any output is a copy of an input, then the input will be reconstructed exactly +/// once all windows are overlapping, which takes `window_length` samples. +/// -Input(s): input signals. +/// -Output(s): processed signals. +/// +/// ### Example: FFT Lowpass Filter +/// ``` +/// use fundsp::hacker::*; +/// let cutoff = 1000.0; +/// let resynth = resynth::(1024, |_time, fft| +/// for i in 0 .. fft.bins() { +/// if fft.frequency(i) <= cutoff { +/// fft.set(0, i, fft.at(0, i)); +/// } +/// }); +/// ``` +pub fn resynth(window_length: usize, processing: F) -> An> +where + I: Size, + O: Size, + F: FnMut(f32, &mut FftWindow) + Clone + Send + Sync, +{ + An(Resynth::new(window_length, processing)) +} diff --git a/src/hacker32.rs b/src/hacker32.rs index 55ad964..6ad2bd2 100644 --- a/src/hacker32.rs +++ b/src/hacker32.rs @@ -1051,7 +1051,7 @@ where O: ConstantFrame, O::Size: Size, { - An(Map::new(f, Routing::Arbitrary)) + An(Map::new(f, Routing::Arbitrary(0.0))) } /// Keeps a signal zero centered. @@ -2173,3 +2173,34 @@ pub fn snoop(capacity: usize) -> (Snoop, An>) { let (snoop, backend) = Snoop::new(capacity); (snoop, An(backend)) } + +/// Frequency domain resynthesizer. The number of inputs is `I` and the number of outputs is `O`. +/// The window length (in samples) must be a power of two and at least four. +/// Processes windows of input samples transformed into the frequency domain. +/// The processing function is issued arguments (time, window). +/// It processes frequency domain inputs into frequency domain outputs. +/// The latency in samples is equal to window length. +/// If any output is a copy of an input, then the input will be reconstructed exactly +/// once all windows are overlapping, which takes `window_length` samples. +/// -Input(s): input signals. +/// -Output(s): processed signals. +/// +/// ### Example: FFT Lowpass Filter +/// ``` +/// use fundsp::hacker32::*; +/// let cutoff = 1000.0; +/// let resynth = resynth::(1024, |_time, fft| +/// for i in 0 .. fft.bins() { +/// if fft.frequency(i) <= cutoff { +/// fft.set(0, i, fft.at(0, i)); +/// } +/// }); +/// ``` +pub fn resynth(window_length: usize, processing: F) -> An> +where + I: Size, + O: Size, + F: FnMut(f32, &mut FftWindow) + Clone + Send + Sync, +{ + An(Resynth::new(window_length, processing)) +} diff --git a/src/prelude.rs b/src/prelude.rs index 99250c6..3e08dd0 100644 --- a/src/prelude.rs +++ b/src/prelude.rs @@ -1091,7 +1091,7 @@ where O: ConstantFrame, O::Size: Size, { - An(Map::new(f, Routing::Arbitrary)) + An(Map::new(f, Routing::Arbitrary(0.0))) } /// Keeps a signal zero centered. @@ -2676,3 +2676,35 @@ pub fn snoop(capacity: usize) -> (Snoop, An>) { let (snoop, backend) = Snoop::new(capacity); (snoop, An(backend)) } + +/// Frequency domain resynthesizer. The number of inputs is `I` and the number of outputs is `O`. +/// The window length (in samples) must be a power of two and at least four. +/// Processes windows of input samples transformed into the frequency domain. +/// The processing function is issued arguments (time, window). +/// It processes frequency domain inputs into frequency domain outputs. +/// The latency in samples is equal to window length. +/// If any output is a copy of an input, then the input will be reconstructed exactly +/// once all windows are overlapping, which takes `window_length` samples. +/// -Input(s): input signals. +/// -Output(s): processed signals. +/// +/// ### Example: FFT Lowpass Filter +/// ``` +/// use fundsp::prelude::*; +/// let cutoff = 1000.0; +/// let resynth = resynth::(1024, |_time, fft| +/// for i in 0 .. fft.bins() { +/// if fft.frequency(i) <= cutoff { +/// fft.set(0, i, fft.at(0, i)); +/// } +/// }); +/// ``` +pub fn resynth(window_length: usize, processing: F) -> An> +where + I: Size, + O: Size, + T: Float, + F: FnMut(f32, &mut FftWindow) + Clone + Send + Sync, +{ + An(Resynth::new(window_length, processing)) +} diff --git a/src/resynth.rs b/src/resynth.rs index 701f0f5..f64bd63 100644 --- a/src/resynth.rs +++ b/src/resynth.rs @@ -1,16 +1,24 @@ -//! Frequency domain resynthesis. WIP. +//! Frequency domain resynthesis. + +// For more information on this technique, see +// "Fourier analysis and reconstruction of audio signals" at +// http://msp.ucsd.edu/techniques/v0.11/book-html/node172.html use super::audionode::*; +use super::math::*; +use super::signal::*; use super::*; use num_complex::Complex32; +use realfft::{ComplexToReal, RealFftPlanner, RealToComplex}; +use std::sync::Arc; -pub struct FftWindow -where - T: Float, - I: Size, - O: Size, -{ - /// Window length. Equals the length of each input and output channel vector. +/// Number of overlapping FFT windows. +const WINDOWS: usize = 4; + +#[derive(Clone)] +pub struct FftWindow { + /// Window length. Must be a power of two and at least four. + /// Equals the length of each input and output channel vector. length: usize, /// Input samples for each input channel. input: Vec>, @@ -22,15 +30,13 @@ where output: Vec>, /// Sample rate for convenience. sample_rate: f32, - _marker: std::marker::PhantomData<(T, I, O)>, + /// Current index into input and output vectors. + index: usize, + /// Total number of processed samples. + samples: u64, } -impl FftWindow -where - T: Float, - I: Size, - O: Size, -{ +impl FftWindow { /// Number of input channels. #[inline] pub fn inputs(&self) -> usize { @@ -43,7 +49,17 @@ where self.output.len() } - /// FFT window length. Must be divisible by four. + /// Get forward vectors for forward FFT. + pub fn forward_vectors(&mut self, channel: usize) -> (&mut Vec, &mut Vec) { + (&mut self.input[channel], &mut self.input_fft[channel]) + } + + /// Get inverse vectors for inverse FFT. + pub fn inverse_vectors(&mut self, channel: usize) -> (&mut Vec, &mut Vec) { + (&mut self.output_fft[channel], &mut self.output[channel]) + } + + /// FFT window length. This is a power of two and at least four. #[inline] pub fn length(&self) -> usize { self.length @@ -72,4 +88,248 @@ where pub fn frequency(&self, i: usize) -> f32 { self.sample_rate / self.length() as f32 * i as f32 } + + /// Create new window. + pub fn new(length: usize, index: usize, inputs: usize, outputs: usize) -> Self { + let mut window = Self { + length, + input: vec![vec!(0.0; length); inputs], + input_fft: Vec::new(), + output_fft: Vec::new(), + output: vec![vec!(0.0; length); outputs], + sample_rate: DEFAULT_SR as f32, + index, + samples: 0, + }; + window + .input_fft + .resize(inputs, vec![Complex32::default(); window.bins()]); + window + .output_fft + .resize(outputs, vec![Complex32::default(); window.bins()]); + window + } + + /// Set the sample rate. + pub fn set_sample_rate(&mut self, sample_rate: f32) { + self.sample_rate = sample_rate; + } + + /// Write input for current index. + #[inline] + pub fn write>(&mut self, input: &Frame, window_value: f32) { + for (channel, item) in input.iter().enumerate() { + self.input[channel][self.index] = item.to_f32() * window_value; + } + } + + /// Read output for current index. + #[inline] + pub fn read>(&self, window_value: f32) -> Frame { + Frame::generate(|channel| convert(self.output[channel][self.index] * window_value)) + } + + /// Set outputs to all zeros. + pub fn clear_output(&mut self) { + for i in 0..self.outputs() { + self.output_fft[i].fill(Complex32::default()); + } + } + + /// Current read and write index. + #[inline] + pub fn index(&self) -> usize { + self.index + } + + /// Reset the window to an empty state. + pub fn reset(&mut self, start_index: usize) { + self.index = start_index; + for channel in 0..self.inputs() { + self.input[channel].fill(0.0); + } + for channel in 0..self.outputs() { + self.output[channel].fill(0.0); + } + } + + /// Advance index to the next sample. + #[inline] + pub fn advance(&mut self) { + self.samples += 1; + self.index = (self.index + 1) & (self.length - 1); + } + + /// Return whether we should do FFT processing right now. + #[inline] + pub fn is_fft_time(&self) -> bool { + self.index == 0 && self.samples >= self.length as u64 + } +} + +/// Frequency domain resynthesizer. Processes windows of input samples with an overlap of four. +/// Each window is Fourier transformed and then processed into output spectra by the processing function. +/// The output windows are then inverse transformed into the outputs. +/// The latency is equal to the window length. +/// If any output is a copy of an input, then the input will be reconstructed exactly once +/// the windows are all overlapping, which takes one window length. +#[derive(Clone)] +pub struct Resynth +where + I: Size, + O: Size, + T: Float, + F: FnMut(f32, &mut FftWindow) + Clone + Send + Sync, +{ + _marker: std::marker::PhantomData<(T, I, O)>, + /// FFT windows. + window: [FftWindow; WINDOWS], + /// Window length. + window_length: usize, + /// Hann window function. + window_function: Vec, + /// Processing function is a function of (time, window). + processing: F, + /// Sample rate. + sample_rate: f64, + /// Forward transform. + forward: Arc>, + /// Inverse transform. + inverse: Arc>, + /// Temporary vector for FFT. + scratch: Vec, + /// Number of processed samples. + samples: u64, +} + +impl Resynth +where + I: Size, + O: Size, + T: Float, + F: FnMut(f32, &mut FftWindow) + Clone + Send + Sync, +{ + /// Create new resynthesizer. + pub fn new(window_length: usize, processing: F) -> Self { + assert!(window_length >= 4 && window_length.is_power_of_two()); + + let mut planner = RealFftPlanner::::new(); + let forward = planner.plan_fft_forward(window_length); + let inverse = planner.plan_fft_inverse(window_length); + + let mut window_function = Vec::new(); + + for i in 0..window_length { + let hann = 0.5 + + 0.5 + * cos((i as i32 - (window_length >> 1) as i32) as f32 * TAU as f32 + / window_length as f32); + window_function.push(hann); + } + + let window = [ + FftWindow::new(window_length, 0, I::USIZE, O::USIZE), + FftWindow::new(window_length, window_length >> 2, I::USIZE, O::USIZE), + FftWindow::new(window_length, window_length >> 1, I::USIZE, O::USIZE), + FftWindow::new( + window_length, + (window_length >> 1) + (window_length >> 2), + I::USIZE, + O::USIZE, + ), + ]; + + let scratch = + vec![Complex32::default(); max(forward.get_scratch_len(), inverse.get_scratch_len())]; + + Self { + _marker: std::marker::PhantomData, + window, + window_length, + window_function, + processing, + sample_rate: DEFAULT_SR, + forward, + inverse, + scratch, + samples: 0, + } + } +} + +impl AudioNode for Resynth +where + I: Size, + O: Size, + T: Float, + F: FnMut(f32, &mut FftWindow) + Clone + Send + Sync, +{ + const ID: u64 = 80; + type Sample = T; + type Inputs = I; + type Outputs = O; + type Setting = (); + + fn set_sample_rate(&mut self, sample_rate: f64) { + self.sample_rate = sample_rate; + for i in 0..WINDOWS { + self.window[i].set_sample_rate(sample_rate as f32); + } + } + + fn reset(&mut self) { + self.samples = 0; + for i in 0..WINDOWS { + self.window[i].reset(i * (self.window_length >> 2)); + } + } + + fn tick( + &mut self, + input: &Frame, + ) -> Frame { + let mut output = Frame::default(); + + let z = 2.0 / (3.0 * self.window_length as f32); + + for i in 0..WINDOWS { + let window_value = self.window_function[self.window[i].index()]; + self.window[i].write(input, window_value); + output += self.window[i].read(window_value * z); + self.window[i].advance(); + } + + self.samples += 1; + + if self.samples & ((self.window_length as u64 >> 2) - 1) == 0 { + for i in 0..WINDOWS { + if self.window[i].is_fft_time() { + for channel in 0..I::USIZE { + let (input, input_fft) = self.window[i].forward_vectors(channel); + self.forward + .process_with_scratch(input, input_fft, &mut self.scratch) + .expect("Internal error"); + } + + self.window[i].clear_output(); + (self.processing)( + (self.samples as f64 / self.sample_rate) as f32, + &mut self.window[i], + ); + + for channel in 0..O::USIZE { + let (output_fft, output) = self.window[i].inverse_vectors(channel); + self.inverse + .process_with_scratch(output_fft, output, &mut self.scratch) + .expect("Internal error"); + } + } + } + } + output + } + + fn route(&mut self, input: &SignalFrame, _frequency: f64) -> SignalFrame { + Routing::Arbitrary(self.window_length as f64).propagate(input, self.outputs()) + } } diff --git a/src/signal.rs b/src/signal.rs index 1c7fa2b..6f57be1 100644 --- a/src/signal.rs +++ b/src/signal.rs @@ -129,8 +129,8 @@ pub fn copy_signal_frame(source: &SignalFrame, i: usize, n: usize) -> SignalFram /// functionality. #[derive(Clone)] pub enum Routing { - /// Conservative routing: every input influences every output nonlinearly. - Arbitrary, + /// Conservative routing: every input influences every output nonlinearly with extra latency. + Arbitrary(f64), /// Split or multisplit semantics. Split, /// Join or multijoin semantics. @@ -147,10 +147,10 @@ impl Routing { return output; } match self { - Routing::Arbitrary => { - let mut combo = input[0].distort(0.0); + Routing::Arbitrary(latency) => { + let mut combo = input[0].distort(*latency); for i in 1..input.len() { - combo = combo.combine_nonlinear(input[i], 0.0); + combo = combo.combine_nonlinear(input[i], *latency); } output.fill(combo); } diff --git a/tests/basic.rs b/tests/basic.rs index 24f7abe..a050210 100644 --- a/tests/basic.rs +++ b/tests/basic.rs @@ -547,3 +547,28 @@ fn test_basic() { assert_eq!(inouts(!zero()), (0, 0)); // A null unit. Stacking it with a graph modifies its sound subtly, as the hash is altered. assert_eq!(inouts(!-!!!--!!!-!!--!zero()), (0, 0)); // Hot-rodded null unit with a custom hash. Uses more electricity. } + +#[test] +/// Test a pass-through resynthesizer. +fn test_resynth() { + let window = 8; + let mut synth: An> = resynth(window, |_time, fft| { + for i in 0..fft.bins() { + fft.set(0, i, fft.at(0, i)); + } + }); + + let duration = 16.0 / DEFAULT_SR; + + let input = Wave64::render(DEFAULT_SR, duration, &mut noise()); + let output = input.filter_latency(duration, &mut synth); + + for i in window..input.length() { + //println!("input {}, output {}", input.at(0, i), output.at(0, i)); + let tolerance = 1.0e-6; + assert!( + input.at(0, i) - tolerance <= output.at(0, i) + && input.at(0, i) + tolerance >= output.at(0, i) + ); + } +}