From ecfb5fe9d3cdf7b61be69bd4617fe587f47850e0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ju=CC=88rgen=20Hock?= Date: Sun, 18 Jun 2023 12:14:10 +0200 Subject: [PATCH] Complete QDFT plan initialization #7 --- cpp/src/qdft/qdft.h | 2 +- rust/examples/analysis.rs | 13 +++- rust/src/lib.rs | 140 +++++++++++++++++++++++++++++++++++--- rust/tests/test.rs | 15 +++- 4 files changed, 157 insertions(+), 13 deletions(-) diff --git a/cpp/src/qdft/qdft.h b/cpp/src/qdft/qdft.h index 58f6fb0..495ecaa 100644 --- a/cpp/src/qdft/qdft.h +++ b/cpp/src/qdft/qdft.h @@ -84,7 +84,7 @@ class QDFT { const std::complex fiddle = std::polar(F(1), F(-2) * pi * (config.quality + k)); - data.fiddles[1 + k] = fiddle; + data.fiddles[k + 1] = fiddle; for (size_t i = 0, j = 1; i < config.size; ++i, j+=3) { diff --git a/rust/examples/analysis.rs b/rust/examples/analysis.rs index 7a765d8..e24066e 100644 --- a/rust/examples/analysis.rs +++ b/rust/examples/analysis.rs @@ -3,7 +3,18 @@ use num::complex::Complex; use qdft::QDFT; fn main() { - let mut qdft = QDFT::new(44100.0, (100.0, 20000.0)); + let samplerate = 44100.0; + let bandwidth = (50.0, samplerate / 2.0); + let resolution = 24.0; + let latency = 0.0; + let window = Some((0.5, -0.5)); + + let mut qdft = QDFT::new( + samplerate, + bandwidth, + resolution, + latency, + window); let n = 1; let mut samples = vec![f32::zero(); n]; diff --git a/rust/src/lib.rs b/rust/src/lib.rs index 5bec4a4..4a6c184 100644 --- a/rust/src/lib.rs +++ b/rust/src/lib.rs @@ -3,15 +3,26 @@ use num::Float; use num::Zero; use num::complex::Complex; -use std::marker::PhantomData; +use std::collections::VecDeque; pub struct QDFT where T: Float, F: Float { samplerate: f64, bandwidth: (f64, f64), - t: PhantomData, - f: PhantomData, + resolution: f64, + latency: f64, + quality: f64, + size: usize, + window: Option<(f64, f64)>, + frequencies: Vec, + periods: Vec, + offsets: Vec, + weights: Vec, + fiddles: Vec>, + twiddles: Vec>, + inputs: VecDeque, + outputs: Vec>, } pub type QDFT32 = QDFT; @@ -20,15 +31,79 @@ pub type QDFT64 = QDFT; impl QDFT where T: Float, F: Float { - pub fn new(samplerate: f64, - bandwidth: (f64, f64) + bandwidth: (f64, f64), + resolution: f64, + latency: f64, + window: Option<(f64, f64)>, ) -> Self { + let quality = f64::powf(f64::powf(2.0, 1.0 / resolution) - 1.0, -1.0); + let size = f64::ceil(resolution * f64::log2(bandwidth.1 / bandwidth.0)) as usize; + + let mut frequencies = vec![f64::zero(); size]; + let mut periods = vec![usize::zero(); size]; + let mut offsets = vec![usize::zero(); size]; + let mut weights = vec![F::zero(); size]; + + for i in 0..size { + let frequency = bandwidth.0 * f64::powf(2.0, (i as f64) / resolution); + frequencies[i] = frequency; + + let period = f64::ceil(quality * samplerate / frequency); + periods[i] = period as usize; + + let offset = f64::ceil(((periods[0] as f64) - period) + * f64::clamp(latency * 0.5 + 0.5, 0.0, 1.0)); + offsets[i] = offset as usize; + + let weight = F::one() / F::from(period).unwrap(); + weights[i] = weight; + } + + let mut fiddles = vec![Complex::::zero(); 3]; + let mut twiddles = vec![Complex::::zero(); size * 3]; + + for k in [-1, 0, 1] { + let pi = std::f64::consts::PI; // acos(-1) + + let fiddle = Complex::::from_polar(F::one(), F::from( + -2.0 * pi * (quality + (k as f64))).unwrap()); + fiddles[(k + 1) as usize] = fiddle; + + let mut i = 0; + let mut j = 1; + + while i < size { + let period = periods[i] as f64; + + let twiddle = Complex::::from_polar(F::one(), F::from( + 2.0 * pi * (quality + (k as f64)) / period).unwrap()); + twiddles[(j + k) as usize] = twiddle; + + i += 1; + j += 3; + } + } + + let inputs = VecDeque::from(vec![T::zero(); periods[0] + 1]); + let outputs = vec![Complex::::zero(); size * 3]; + QDFT { - samplerate: samplerate, - bandwidth: bandwidth, - t: PhantomData, - f: PhantomData, + samplerate, + bandwidth, + resolution, + latency, + quality, + size, + window, + frequencies, + periods, + offsets, + weights, + fiddles, + twiddles, + inputs, + outputs, } } @@ -64,3 +139,50 @@ impl QDFT self.iqdft_vector(dfts, samples); } } + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_sanity_checks() { + let samplerate = 44100.0; + let bandwidth = (50.0, samplerate / 2.0); + let resolution = 24.0; + let latency = 0.0; + let window = Some((0.5, -0.5)); + + let qdft = QDFT::::new( + samplerate, + bandwidth, + resolution, + latency, + window); + + let mut i = 0; + let mut j = 1; + + while i < qdft.size { + assert!((qdft.offsets[i]) < qdft.inputs.len()); + assert!((qdft.offsets[i] + qdft.periods[i]) < qdft.inputs.len()); + + assert!((-1 + 1) < qdft.fiddles.len() as isize); + assert!(( 0 + 1) < qdft.fiddles.len() as isize); + assert!(( 1 + 1) < qdft.fiddles.len() as isize); + + assert!((-1 + j) < qdft.twiddles.len() as isize); + assert!(( 0 + j) < qdft.twiddles.len() as isize); + assert!(( 1 + j) < qdft.twiddles.len() as isize); + + assert!((-1 + j) < qdft.outputs.len() as isize); + assert!(( 0 + j) < qdft.outputs.len() as isize); + assert!(( 1 + j) < qdft.outputs.len() as isize); + + i += 1; + j += 3; + } + + assert_eq!(i, qdft.size); + assert_eq!(j, (qdft.size as isize) * 3 + 1); + } +} diff --git a/rust/tests/test.rs b/rust/tests/test.rs index 16f9017..0a88c94 100644 --- a/rust/tests/test.rs +++ b/rust/tests/test.rs @@ -7,8 +7,19 @@ mod tests { use qdft::QDFT; #[test] - fn test() { - let mut qdft = QDFT::new(44100.0, (100.0, 20000.0)); + fn test_one_in_zero_out() { + let samplerate = 44100.0; + let bandwidth = (50.0, samplerate / 2.0); + let resolution = 24.0; + let latency = 0.0; + let window = Some((0.5, -0.5)); + + let mut qdft = QDFT::new( + samplerate, + bandwidth, + resolution, + latency, + window); let n = 1; let mut samples = vec![f32::zero(); n];