From bad50ff5e200aab957ea1ddfa11f9e2ced41bff4 Mon Sep 17 00:00:00 2001 From: Jasper Date: Fri, 12 Mar 2021 16:54:18 +0800 Subject: [PATCH 01/13] migrate libtest to criterion; migrate rustfft to 3.0.1 --- Cargo.toml | 9 +++++- benches/lib.rs | 74 ++++++++++++++++++++++++++------------------------ src/lib.rs | 2 +- 3 files changed, 47 insertions(+), 38 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 3bc28e6..e855c45 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -13,5 +13,12 @@ version = "0.2.0" [dependencies] apodize = "0.3.1" num = "0.2.0" -rustfft = "3.0.0" +rustfft = "3.0.1" strider = "0.1.3" + +[dev-dependencies] +criterion = "0.3" + +[[bench]] +name = "lib" +harness = false \ No newline at end of file diff --git a/benches/lib.rs b/benches/lib.rs index f893b8e..67dd5ec 100644 --- a/benches/lib.rs +++ b/benches/lib.rs @@ -1,73 +1,73 @@ -#![feature(test)] -extern crate test; +extern crate criterion; +use criterion::{criterion_group, criterion_main, Criterion}; extern crate num; use num::complex::Complex; extern crate rustfft; -use rustfft::FFT; +use rustfft::FFTplanner; extern crate stft; use stft::{STFT, WindowType}; macro_rules! bench_fft_process { - ($bencher:expr, $window_size:expr, $float:ty) => {{ + ($c:expr, $window_size:expr, $float:ty) => {{ let inverse = false; - let window_size = $window_size; - let mut fft = FFT::<$float>::new(window_size, inverse); - let input = std::iter::repeat(Complex::new(0., 0.)) - .take(window_size) + let mut planner = FFTplanner::new(inverse); + let fft = planner.plan_fft($window_size); + let mut input = std::iter::repeat(Complex::new(0., 0.)) + .take($window_size) .collect::>>(); let mut output = std::iter::repeat(Complex::new(0., 0.)) - .take(window_size) + .take($window_size) .collect::>>(); - $bencher.iter(|| { - fft.process(&input[..], &mut output[..]) - }); + $c.bench_function(concat!("bench_fft_process_", stringify!($window_size), "_", stringify!($float)), |b| b.iter(|| { + fft.process(&mut input[..], &mut output[..]) + })); }} } -#[bench] -fn bench_fft_process_1024_f32(bencher: &mut test::Bencher) { - bench_fft_process!(bencher, 1024, f32); +fn bench_fft_process_1024_f32(c: &mut Criterion) { + bench_fft_process!(c, 1024, f32); } -#[bench] -fn bench_fft_process_1024_f64(bencher: &mut test::Bencher) { - bench_fft_process!(bencher, 1024, f64); +fn bench_fft_process_1024_f64(c: &mut Criterion) { + bench_fft_process!(c, 1024, f64); } +criterion_group!(benches_fft_process, bench_fft_process_1024_f32, bench_fft_process_1024_f64); + macro_rules! bench_stft_compute { - ($bencher:expr, $window_size:expr, $float:ty) => {{ + ($c:expr, $window_size:expr, $float:ty) => {{ let mut stft = STFT::<$float>::new(WindowType::Hanning, $window_size, 0); let input = std::iter::repeat(1.).take($window_size).collect::>(); let mut output = std::iter::repeat(0.).take(stft.output_size()).collect::>(); stft.append_samples(&input[..]); - $bencher.iter(|| { + $c.bench_function(concat!("bench_stft_compute_", stringify!($window_size), "_", stringify!($float)), |b| b.iter(|| { stft.compute_column(&mut output[..]) - }); + })); }} } -#[bench] -fn bench_stft_compute_1024_f32(bencher: &mut test::Bencher) { - bench_stft_compute!(bencher, 1024, f32); +fn bench_stft_compute_1024_f32(c: &mut Criterion) { + bench_stft_compute!(c, 1024, f32); } -#[bench] -fn bench_stft_compute_1024_f64(bencher: &mut test::Bencher) { - bench_stft_compute!(bencher, 1024, f64); +fn bench_stft_compute_1024_f64(c: &mut Criterion) { + bench_stft_compute!(c, 1024, f64); } +criterion_group!(benches_stft_compute, bench_stft_compute_1024_f32, bench_stft_compute_1024_f64); + macro_rules! bench_stft_audio { - ($bencher:expr, $seconds:expr, $float:ty) => {{ + ($c:expr, $seconds:expr, $float:ty) => {{ // let's generate some fake audio let sample_rate: usize = 44100; let seconds: usize = $seconds; let sample_count = sample_rate * seconds; let all_samples = (0..sample_count).map(|x| x as $float).collect::>(); - $bencher.iter(|| { + $c.bench_function(concat!("bench_stft_audio_", stringify!($windowsize), "_", stringify!($float)), |b| b.iter(|| { // let's initialize our short-time fourier transform let window_type: WindowType = WindowType::Hanning; let window_size: usize = 1024; @@ -83,16 +83,18 @@ macro_rules! bench_stft_audio { stft.move_to_next_column(); } } - }); + })); }} } -#[bench] -fn bench_stft_10_seconds_audio_f32(bencher: &mut test::Bencher) { - bench_stft_audio!(bencher, 10, f32); +fn bench_stft_10_seconds_audio_f32(c: &mut Criterion) { + bench_stft_audio!(c, 10, f32); } -#[bench] -fn bench_stft_10_seconds_audio_f64(bencher: &mut test::Bencher) { - bench_stft_audio!(bencher, 10, f64); +fn bench_stft_10_seconds_audio_f64(c: &mut Criterion) { + bench_stft_audio!(c, 10, f64); } + +criterion_group!(benches_stft_audio, bench_stft_10_seconds_audio_f32, bench_stft_10_seconds_audio_f64); + +criterion_main!(benches_fft_process, benches_stft_compute, benches_stft_audio); \ No newline at end of file diff --git a/src/lib.rs b/src/lib.rs index 39c1c2b..603da34 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -146,7 +146,7 @@ pub struct STFT { pub window_size: usize, pub step_size: usize, - pub fft: Arc>, + pub fft: Arc>, pub window: Option>, /// internal ringbuffer used to store samples pub sample_ring: SliceRingImpl, From 3695eb86ebd7e7c7cd5949faf939c52d25fbd7d7 Mon Sep 17 00:00:00 2001 From: Jasper Date: Sat, 13 Mar 2021 14:14:20 +0800 Subject: [PATCH 02/13] add zero-padding FFT --- Cargo.toml | 2 +- src/lib.rs | 65 +++++++++++++++++++++++++++++++++++++++--------------- 2 files changed, 48 insertions(+), 19 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index e855c45..4a63425 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -8,7 +8,7 @@ license = "MIT OR Apache-2.0" name = "stft" readme = "README.md" repository = "https://github.com/snd/stft.git" -version = "0.2.0" +version = "0.3.0" [dependencies] apodize = "0.3.1" diff --git a/src/lib.rs b/src/lib.rs index 603da34..02ab366 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -145,6 +145,7 @@ pub struct STFT where T: FFTnum + FromF64 + num::Float { pub window_size: usize, + pub fft_size: usize, pub step_size: usize, pub fft: Arc>, pub window: Option>, @@ -174,34 +175,42 @@ impl STFT let window = Self::window_type_to_window_vec(window_type, window_size); Self::new_with_window_vec(window, window_size, step_size) } + + pub fn new_with_zero_padding(window_type: WindowType, window_size: usize, fft_size: usize, step_size: usize) -> Self { + let window = Self::window_type_to_window_vec(window_type, window_size); + Self::new_with_window_vec_and_zero_padding(window, window_size, fft_size, step_size) + } // TODO this should ideally take an iterator and not a vec - pub fn new_with_window_vec(window: Option>, - window_size: usize, - step_size: usize) - -> Self { - // TODO more assertions: - // window_size is power of two - // step_size > 0 - assert!(step_size <= window_size); - let inverse = false; - let mut planner = FFTplanner::new(inverse); + pub fn new_with_window_vec_and_zero_padding(window:Option>, + window_size: usize, + fft_size: usize, + step_size: usize) -> Self { + assert!(step_size > 0 && step_size < window_size); + let mut planner = FFTplanner::new(false); STFT { window_size: window_size, + fft_size: fft_size, step_size: step_size, - fft: planner.plan_fft(window_size), + fft: planner.plan_fft(fft_size), sample_ring: SliceRingImpl::new(), window: window, real_input: std::iter::repeat(T::zero()) - .take(window_size) - .collect(), + .take(window_size).collect(), + // zero-padded complex_input, so the size is fft_size, not window_size complex_input: std::iter::repeat(Complex::::zero()) - .take(window_size) - .collect(), + .take(fft_size).collect(), + // same size as complex_output complex_output: std::iter::repeat(Complex::::zero()) - .take(window_size) - .collect(), + .take(fft_size).collect(), } + + } + + pub fn new_with_window_vec(window: Option>, + window_size: usize, + step_size: usize) -> Self { + Self::new_with_window_vec_and_zero_padding(window, window_size, window_size, step_size) } #[inline] @@ -237,7 +246,8 @@ impl STFT } // copy windowed real_input as real parts into complex_input - for (dst, src) in self.complex_input.iter_mut().zip(self.real_input.iter()) { + // only copy `window_size` size, leave the rest in `complex_input` be zero + for (src, dst) in self.real_input.iter().zip(self.complex_input.iter_mut()) { dst.re = src.clone(); } @@ -287,6 +297,25 @@ impl STFT pub fn move_to_next_column(&mut self) { self.sample_ring.drop_many_front(self.step_size); } + + /// corresponding frequencies of a column of the spectrogram + /// # Arguments + /// `fs`: sampling frequency. + pub fn freqs(&self, fs: f64) -> Vec { + let n_freqs = self.output_size(); + (0..n_freqs).map(|f| (f as f64) / ((n_freqs - 1) as f64) * (fs / 2.)) + .collect() + } + + /// corresponding time of first columns of the spectrogram + pub fn first_time(&self, fs: f64) -> f64 { + (self.window_size as f64) / (fs * 2.) + } + + /// time interval between two adjacent columns of the spectrogram + pub fn time_interval(&self, fs: f64) -> f64 { + (self.step_size as f64) / fs + } } pub trait FromF64 { From a6a46f695b52c5ad5fad29baa455250d6da434ce Mon Sep 17 00:00:00 2001 From: Markus Mayer Date: Wed, 13 Oct 2021 20:09:03 +0200 Subject: [PATCH 03/13] Update project to Rust 2018 edition --- Cargo.toml | 3 +- benches/lib.rs | 122 +++++++++++++++++++++++++++++---------------- src/lib.rs | 131 +++++++++++++++++++++++++++++-------------------- tests/lib.rs | 25 +++++++--- 4 files changed, 178 insertions(+), 103 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 4a63425..983f72f 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -9,6 +9,7 @@ name = "stft" readme = "README.md" repository = "https://github.com/snd/stft.git" version = "0.3.0" +edition = "2018" [dependencies] apodize = "0.3.1" @@ -21,4 +22,4 @@ criterion = "0.3" [[bench]] name = "lib" -harness = false \ No newline at end of file +harness = false diff --git a/benches/lib.rs b/benches/lib.rs index 67dd5ec..373c7ea 100644 --- a/benches/lib.rs +++ b/benches/lib.rs @@ -1,14 +1,7 @@ -extern crate criterion; use criterion::{criterion_group, criterion_main, Criterion}; - -extern crate num; use num::complex::Complex; - -extern crate rustfft; use rustfft::FFTplanner; - -extern crate stft; -use stft::{STFT, WindowType}; +use stft::{WindowType, STFT}; macro_rules! bench_fft_process { ($c:expr, $window_size:expr, $float:ty) => {{ @@ -21,11 +14,16 @@ macro_rules! bench_fft_process { let mut output = std::iter::repeat(Complex::new(0., 0.)) .take($window_size) .collect::>>(); - $c.bench_function(concat!("bench_fft_process_", stringify!($window_size), "_", stringify!($float)), |b| b.iter(|| { - fft.process(&mut input[..], &mut output[..]) - })); - - }} + $c.bench_function( + concat!( + "bench_fft_process_", + stringify!($window_size), + "_", + stringify!($float) + ), + |b| b.iter(|| fft.process(&mut input[..], &mut output[..])), + ); + }}; } fn bench_fft_process_1024_f32(c: &mut Criterion) { @@ -36,18 +34,32 @@ fn bench_fft_process_1024_f64(c: &mut Criterion) { bench_fft_process!(c, 1024, f64); } -criterion_group!(benches_fft_process, bench_fft_process_1024_f32, bench_fft_process_1024_f64); +criterion_group!( + benches_fft_process, + bench_fft_process_1024_f32, + bench_fft_process_1024_f64 +); macro_rules! bench_stft_compute { ($c:expr, $window_size:expr, $float:ty) => {{ let mut stft = STFT::<$float>::new(WindowType::Hanning, $window_size, 0); - let input = std::iter::repeat(1.).take($window_size).collect::>(); - let mut output = std::iter::repeat(0.).take(stft.output_size()).collect::>(); + let input = std::iter::repeat(1.) + .take($window_size) + .collect::>(); + let mut output = std::iter::repeat(0.) + .take(stft.output_size()) + .collect::>(); stft.append_samples(&input[..]); - $c.bench_function(concat!("bench_stft_compute_", stringify!($window_size), "_", stringify!($float)), |b| b.iter(|| { - stft.compute_column(&mut output[..]) - })); - }} + $c.bench_function( + concat!( + "bench_stft_compute_", + stringify!($window_size), + "_", + stringify!($float) + ), + |b| b.iter(|| stft.compute_column(&mut output[..])), + ); + }}; } fn bench_stft_compute_1024_f32(c: &mut Criterion) { @@ -58,7 +70,11 @@ fn bench_stft_compute_1024_f64(c: &mut Criterion) { bench_stft_compute!(c, 1024, f64); } -criterion_group!(benches_stft_compute, bench_stft_compute_1024_f32, bench_stft_compute_1024_f64); +criterion_group!( + benches_stft_compute, + bench_stft_compute_1024_f32, + bench_stft_compute_1024_f64 +); macro_rules! bench_stft_audio { ($c:expr, $seconds:expr, $float:ty) => {{ @@ -66,25 +82,37 @@ macro_rules! bench_stft_audio { let sample_rate: usize = 44100; let seconds: usize = $seconds; let sample_count = sample_rate * seconds; - let all_samples = (0..sample_count).map(|x| x as $float).collect::>(); - $c.bench_function(concat!("bench_stft_audio_", stringify!($windowsize), "_", stringify!($float)), |b| b.iter(|| { - // let's initialize our short-time fourier transform - let window_type: WindowType = WindowType::Hanning; - let window_size: usize = 1024; - let step_size: usize = 512; - let mut stft = STFT::<$float>::new(window_type, window_size, step_size); - // we need a buffer to hold a computed column of the spectrogram - let mut spectrogram_column: Vec<$float> = - std::iter::repeat(0.).take(stft.output_size()).collect(); - for some_samples in (&all_samples[..]).chunks(3000) { - stft.append_samples(some_samples); - while stft.contains_enough_to_compute() { - stft.compute_column(&mut spectrogram_column[..]); - stft.move_to_next_column(); - } - } - })); - }} + let all_samples = (0..sample_count) + .map(|x| x as $float) + .collect::>(); + $c.bench_function( + concat!( + "bench_stft_audio_", + stringify!($windowsize), + "_", + stringify!($float) + ), + |b| { + b.iter(|| { + // let's initialize our short-time fourier transform + let window_type: WindowType = WindowType::Hanning; + let window_size: usize = 1024; + let step_size: usize = 512; + let mut stft = STFT::<$float>::new(window_type, window_size, step_size); + // we need a buffer to hold a computed column of the spectrogram + let mut spectrogram_column: Vec<$float> = + std::iter::repeat(0.).take(stft.output_size()).collect(); + for some_samples in (&all_samples[..]).chunks(3000) { + stft.append_samples(some_samples); + while stft.contains_enough_to_compute() { + stft.compute_column(&mut spectrogram_column[..]); + stft.move_to_next_column(); + } + } + }) + }, + ); + }}; } fn bench_stft_10_seconds_audio_f32(c: &mut Criterion) { @@ -95,6 +123,14 @@ fn bench_stft_10_seconds_audio_f64(c: &mut Criterion) { bench_stft_audio!(c, 10, f64); } -criterion_group!(benches_stft_audio, bench_stft_10_seconds_audio_f32, bench_stft_10_seconds_audio_f64); - -criterion_main!(benches_fft_process, benches_stft_compute, benches_stft_audio); \ No newline at end of file +criterion_group!( + benches_stft_audio, + bench_stft_10_seconds_audio_f32, + bench_stft_10_seconds_audio_f64 +); + +criterion_main!( + benches_fft_process, + benches_stft_compute, + benches_stft_audio +); diff --git a/src/lib.rs b/src/lib.rs index 02ab366..914fbfe 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -9,7 +9,6 @@ to the `[dependencies]` section of your `Cargo.toml` and call `extern crate stft ## example ``` -extern crate stft; use stft::{STFT, WindowType}; fn main() { @@ -59,21 +58,13 @@ fn main() { ``` */ -use std::str::FromStr; -use std::sync::Arc; - -extern crate num; use num::complex::Complex; use num::traits::{Float, Signed, Zero}; - -extern crate apodize; - -extern crate strider; +use rustfft::{FFTnum, FFTplanner, FFT}; +use std::str::FromStr; +use std::sync::Arc; use strider::{SliceRing, SliceRingImpl}; -extern crate rustfft; -use rustfft::{FFT,FFTnum,FFTplanner}; - /// returns `0` if `log10(value).is_negative()`. /// otherwise returns `log10(value)`. /// `log10` turns values in domain `0..1` into values @@ -109,7 +100,7 @@ impl FromStr for WindowType { fn from_str(s: &str) -> Result { let lower = s.to_lowercase(); - match &lower[..] { + match lower.as_str() { "hanning" => Ok(WindowType::Hanning), "hann" => Ok(WindowType::Hanning), "hamming" => Ok(WindowType::Hamming), @@ -129,11 +120,13 @@ impl std::fmt::Display for WindowType { } // TODO write a macro that does this automatically for any enum -static WINDOW_TYPES: [WindowType; 5] = [WindowType::Hanning, - WindowType::Hamming, - WindowType::Blackman, - WindowType::Nuttall, - WindowType::None]; +static WINDOW_TYPES: [WindowType; 5] = [ + WindowType::Hanning, + WindowType::Hamming, + WindowType::Blackman, + WindowType::Nuttall, + WindowType::None, +]; impl WindowType { pub fn values() -> [WindowType; 5] { @@ -142,7 +135,8 @@ impl WindowType { } pub struct STFT - where T: FFTnum + FromF64 + num::Float +where + T: FFTnum + FromF64 + num::Float, { pub window_size: usize, pub fft_size: usize, @@ -157,16 +151,34 @@ pub struct STFT } impl STFT - where T: FFTnum + FromF64 + num::Float +where + T: FFTnum + FromF64 + num::Float, { - pub fn window_type_to_window_vec(window_type: WindowType, - window_size: usize) - -> Option> { + pub fn window_type_to_window_vec( + window_type: WindowType, + window_size: usize, + ) -> Option> { match window_type { - WindowType::Hanning => Some(apodize::hanning_iter(window_size).map(FromF64::from_f64).collect()), - WindowType::Hamming => Some(apodize::hamming_iter(window_size).map(FromF64::from_f64).collect()), - WindowType::Blackman => Some(apodize::blackman_iter(window_size).map(FromF64::from_f64).collect()), - WindowType::Nuttall => Some(apodize::nuttall_iter(window_size).map(FromF64::from_f64).collect()), + WindowType::Hanning => Some( + apodize::hanning_iter(window_size) + .map(FromF64::from_f64) + .collect(), + ), + WindowType::Hamming => Some( + apodize::hamming_iter(window_size) + .map(FromF64::from_f64) + .collect(), + ), + WindowType::Blackman => Some( + apodize::blackman_iter(window_size) + .map(FromF64::from_f64) + .collect(), + ), + WindowType::Nuttall => Some( + apodize::nuttall_iter(window_size) + .map(FromF64::from_f64) + .collect(), + ), WindowType::None => None, } } @@ -175,41 +187,50 @@ impl STFT let window = Self::window_type_to_window_vec(window_type, window_size); Self::new_with_window_vec(window, window_size, step_size) } - - pub fn new_with_zero_padding(window_type: WindowType, window_size: usize, fft_size: usize, step_size: usize) -> Self { + + pub fn new_with_zero_padding( + window_type: WindowType, + window_size: usize, + fft_size: usize, + step_size: usize, + ) -> Self { let window = Self::window_type_to_window_vec(window_type, window_size); Self::new_with_window_vec_and_zero_padding(window, window_size, fft_size, step_size) } // TODO this should ideally take an iterator and not a vec - pub fn new_with_window_vec_and_zero_padding(window:Option>, - window_size: usize, - fft_size: usize, - step_size: usize) -> Self { + pub fn new_with_window_vec_and_zero_padding( + window: Option>, + window_size: usize, + fft_size: usize, + step_size: usize, + ) -> Self { assert!(step_size > 0 && step_size < window_size); let mut planner = FFTplanner::new(false); STFT { - window_size: window_size, - fft_size: fft_size, - step_size: step_size, + window_size, + fft_size, + step_size, fft: planner.plan_fft(fft_size), sample_ring: SliceRingImpl::new(), - window: window, - real_input: std::iter::repeat(T::zero()) - .take(window_size).collect(), + window, + real_input: std::iter::repeat(T::zero()).take(window_size).collect(), // zero-padded complex_input, so the size is fft_size, not window_size complex_input: std::iter::repeat(Complex::::zero()) - .take(fft_size).collect(), - // same size as complex_output + .take(fft_size) + .collect(), + // same size as complex_output complex_output: std::iter::repeat(Complex::::zero()) - .take(fft_size).collect(), + .take(fft_size) + .collect(), } - } - pub fn new_with_window_vec(window: Option>, - window_size: usize, - step_size: usize) -> Self { + pub fn new_with_window_vec( + window: Option>, + window_size: usize, + step_size: usize, + ) -> Self { Self::new_with_window_vec_and_zero_padding(window, window_size, window_size, step_size) } @@ -246,13 +267,14 @@ impl STFT } // copy windowed real_input as real parts into complex_input - // only copy `window_size` size, leave the rest in `complex_input` be zero + // only copy `window_size` size, leave the rest in `complex_input` be zero for (src, dst) in self.real_input.iter().zip(self.complex_input.iter_mut()) { dst.re = src.clone(); } // compute fft - self.fft.process(&mut self.complex_input, &mut self.complex_output); + self.fft + .process(&mut self.complex_input, &mut self.complex_output); } /// # Panics @@ -303,13 +325,14 @@ impl STFT /// `fs`: sampling frequency. pub fn freqs(&self, fs: f64) -> Vec { let n_freqs = self.output_size(); - (0..n_freqs).map(|f| (f as f64) / ((n_freqs - 1) as f64) * (fs / 2.)) + (0..n_freqs) + .map(|f| (f as f64) / ((n_freqs - 1) as f64) * (fs / 2.)) .collect() } /// corresponding time of first columns of the spectrogram pub fn first_time(&self, fs: f64) -> f64 { - (self.window_size as f64) / (fs * 2.) + (self.window_size as f64) / (fs * 2.) } /// time interval between two adjacent columns of the spectrogram @@ -323,9 +346,13 @@ pub trait FromF64 { } impl FromF64 for f64 { - fn from_f64(n: f64) -> Self { n } + fn from_f64(n: f64) -> Self { + n + } } impl FromF64 for f32 { - fn from_f64(n: f64) -> Self { n as f32 } + fn from_f64(n: f64) -> Self { + n as f32 + } } diff --git a/tests/lib.rs b/tests/lib.rs index 3fdc11b..ea5f393 100644 --- a/tests/lib.rs +++ b/tests/lib.rs @@ -1,14 +1,21 @@ use std::str::FromStr; - -extern crate stft; -use stft::{STFT, WindowType}; +use stft::{WindowType, STFT}; #[test] fn test_window_type_from_string() { - assert_eq!(WindowType::from_str("Hanning").unwrap(), WindowType::Hanning); - assert_eq!(WindowType::from_str("hanning").unwrap(), WindowType::Hanning); + assert_eq!( + WindowType::from_str("Hanning").unwrap(), + WindowType::Hanning + ); + assert_eq!( + WindowType::from_str("hanning").unwrap(), + WindowType::Hanning + ); assert_eq!(WindowType::from_str("hann").unwrap(), WindowType::Hanning); - assert_eq!(WindowType::from_str("blackman").unwrap(), WindowType::Blackman); + assert_eq!( + WindowType::from_str("blackman").unwrap(), + WindowType::Blackman + ); } #[test] @@ -20,7 +27,11 @@ fn test_window_type_to_string() { fn test_window_types_to_strings() { assert_eq!( vec!["Hanning", "Hamming", "Blackman", "Nuttall", "None"], - WindowType::values().iter().map(|x| x.to_string()).collect::>()); + WindowType::values() + .iter() + .map(|x| x.to_string()) + .collect::>() + ); } #[test] From 1d7e28d113c68bdb17812af23fe8ea42246c0f43 Mon Sep 17 00:00:00 2001 From: Markus Mayer Date: Wed, 13 Oct 2021 20:32:19 +0200 Subject: [PATCH 04/13] Fix an invalid STFT step size in benchmark --- benches/lib.rs | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/benches/lib.rs b/benches/lib.rs index 373c7ea..a77e7d7 100644 --- a/benches/lib.rs +++ b/benches/lib.rs @@ -42,7 +42,8 @@ criterion_group!( macro_rules! bench_stft_compute { ($c:expr, $window_size:expr, $float:ty) => {{ - let mut stft = STFT::<$float>::new(WindowType::Hanning, $window_size, 0); + let step_size: usize = 512; + let mut stft = STFT::<$float>::new(WindowType::Hanning, $window_size, step_size); let input = std::iter::repeat(1.) .take($window_size) .collect::>(); From 3407a9eac04f5855880aeefa576358436c78dc7c Mon Sep 17 00:00:00 2001 From: Markus Mayer Date: Wed, 13 Oct 2021 20:38:03 +0200 Subject: [PATCH 05/13] Update rustfft to 6.0.1 --- Cargo.toml | 6 +++--- benches/lib.rs | 13 +++++-------- src/lib.rs | 21 ++++++++++++++------- 3 files changed, 22 insertions(+), 18 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 983f72f..7a79d55 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -13,12 +13,12 @@ edition = "2018" [dependencies] apodize = "0.3.1" -num = "0.2.0" -rustfft = "3.0.1" +num = "0.4.0" +rustfft = "6.0.1" strider = "0.1.3" [dev-dependencies] -criterion = "0.3" +criterion = "0.3.4" [[bench]] name = "lib" diff --git a/benches/lib.rs b/benches/lib.rs index a77e7d7..6b6d752 100644 --- a/benches/lib.rs +++ b/benches/lib.rs @@ -1,16 +1,13 @@ use criterion::{criterion_group, criterion_main, Criterion}; use num::complex::Complex; -use rustfft::FFTplanner; +use rustfft::{FftDirection, FftPlanner}; use stft::{WindowType, STFT}; macro_rules! bench_fft_process { ($c:expr, $window_size:expr, $float:ty) => {{ - let inverse = false; - let mut planner = FFTplanner::new(inverse); - let fft = planner.plan_fft($window_size); - let mut input = std::iter::repeat(Complex::new(0., 0.)) - .take($window_size) - .collect::>>(); + let mut planner = FftPlanner::new(); + let fft = planner.plan_fft($window_size, FftDirection::Forward); + // input is processed in-place let mut output = std::iter::repeat(Complex::new(0., 0.)) .take($window_size) .collect::>>(); @@ -21,7 +18,7 @@ macro_rules! bench_fft_process { "_", stringify!($float) ), - |b| b.iter(|| fft.process(&mut input[..], &mut output[..])), + |b| b.iter(|| fft.process(&mut output[..])), ); }}; } diff --git a/src/lib.rs b/src/lib.rs index 914fbfe..af5180a 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -60,7 +60,7 @@ fn main() { use num::complex::Complex; use num::traits::{Float, Signed, Zero}; -use rustfft::{FFTnum, FFTplanner, FFT}; +use rustfft::{Fft, FftDirection, FftNum, FftPlanner}; use std::str::FromStr; use std::sync::Arc; use strider::{SliceRing, SliceRingImpl}; @@ -136,23 +136,24 @@ impl WindowType { pub struct STFT where - T: FFTnum + FromF64 + num::Float, + T: FftNum + FromF64 + num::Float, { pub window_size: usize, pub fft_size: usize, pub step_size: usize, - pub fft: Arc>, + pub fft: Arc>, pub window: Option>, /// internal ringbuffer used to store samples pub sample_ring: SliceRingImpl, pub real_input: Vec, pub complex_input: Vec>, pub complex_output: Vec>, + fft_scratch: Vec>, } impl STFT where - T: FFTnum + FromF64 + num::Float, + T: FftNum + FromF64 + num::Float, { pub fn window_type_to_window_vec( window_type: WindowType, @@ -206,12 +207,18 @@ where step_size: usize, ) -> Self { assert!(step_size > 0 && step_size < window_size); - let mut planner = FFTplanner::new(false); + let fft = FftPlanner::new().plan_fft(fft_size, FftDirection::Forward); + + // allocate a scratch buffer for the FFT + let scratch_len = fft.get_inplace_scratch_len(); + let fft_scratch = vec![Complex::::zero(); scratch_len]; + STFT { window_size, fft_size, step_size, - fft: planner.plan_fft(fft_size), + fft, + fft_scratch, sample_ring: SliceRingImpl::new(), window, real_input: std::iter::repeat(T::zero()).take(window_size).collect(), @@ -274,7 +281,7 @@ where // compute fft self.fft - .process(&mut self.complex_input, &mut self.complex_output); + .process_with_scratch(&mut self.complex_output, &mut self.fft_scratch) } /// # Panics From b96ee400b7ea3f2a1a963d686413eaea825f1846 Mon Sep 17 00:00:00 2001 From: Markus Mayer Date: Wed, 13 Oct 2021 20:38:16 +0200 Subject: [PATCH 06/13] Update apodize to 1.0.0 --- Cargo.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Cargo.toml b/Cargo.toml index 7a79d55..f800952 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -12,7 +12,7 @@ version = "0.3.0" edition = "2018" [dependencies] -apodize = "0.3.1" +apodize = "1.0.0" num = "0.4.0" rustfft = "6.0.1" strider = "0.1.3" From 8547443b9bb805a967e8d8c079c3cfe7be5ef04f Mon Sep 17 00:00:00 2001 From: Markus Mayer Date: Wed, 13 Oct 2021 20:39:06 +0200 Subject: [PATCH 07/13] Update .gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index 7883767..f1d8a0d 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ target **.DS_Store Cargo.lock +.idea/ From fef79d34783c05bb490e8e63c8a51871cd5bc331 Mon Sep 17 00:00:00 2001 From: Markus Mayer Date: Wed, 13 Oct 2021 21:04:19 +0200 Subject: [PATCH 08/13] Fix "needless `fn main` in doctest" clippy warning --- src/lib.rs | 82 ++++++++++++++++++++++++++---------------------------- 1 file changed, 40 insertions(+), 42 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index af5180a..113b8b0 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -11,48 +11,46 @@ to the `[dependencies]` section of your `Cargo.toml` and call `extern crate stft ``` use stft::{STFT, WindowType}; -fn main() { - // let's generate ten seconds of fake audio - let sample_rate: usize = 44100; - let seconds: usize = 10; - let sample_count = sample_rate * seconds; - let all_samples = (0..sample_count).map(|x| x as f64).collect::>(); - - // let's initialize our short-time fourier transform - let window_type: WindowType = WindowType::Hanning; - let window_size: usize = 1024; - let step_size: usize = 512; - let mut stft = STFT::new(window_type, window_size, step_size); - - // we need a buffer to hold a computed column of the spectrogram - let mut spectrogram_column: Vec = - std::iter::repeat(0.).take(stft.output_size()).collect(); - - // iterate over all the samples in chunks of 3000 samples. - // in a real program you would probably read from something instead. - for some_samples in (&all_samples[..]).chunks(3000) { - // append the samples to the internal ringbuffer of the stft - stft.append_samples(some_samples); - - // as long as there remain window_size samples in the internal - // ringbuffer of the stft - while stft.contains_enough_to_compute() { - // compute one column of the stft by - // taking the first window_size samples of the internal ringbuffer, - // multiplying them with the window, - // computing the fast fourier transform, - // taking half of the symetric complex outputs, - // computing the norm of the complex outputs and - // taking the log10 - stft.compute_column(&mut spectrogram_column[..]); - - // here's where you would do something with the - // spectrogram_column... - - // drop step_size samples from the internal ringbuffer of the stft - // making a step of size step_size - stft.move_to_next_column(); - } +// let's generate ten seconds of fake audio +let sample_rate: usize = 44100; +let seconds: usize = 10; +let sample_count = sample_rate * seconds; +let all_samples = (0..sample_count).map(|x| x as f64).collect::>(); + +// let's initialize our short-time fourier transform +let window_type: WindowType = WindowType::Hanning; +let window_size: usize = 1024; +let step_size: usize = 512; +let mut stft = STFT::new(window_type, window_size, step_size); + +// we need a buffer to hold a computed column of the spectrogram +let mut spectrogram_column: Vec = + std::iter::repeat(0.).take(stft.output_size()).collect(); + +// iterate over all the samples in chunks of 3000 samples. +// in a real program you would probably read from something instead. +for some_samples in (&all_samples[..]).chunks(3000) { + // append the samples to the internal ringbuffer of the stft + stft.append_samples(some_samples); + + // as long as there remain window_size samples in the internal + // ringbuffer of the stft + while stft.contains_enough_to_compute() { + // compute one column of the stft by + // taking the first window_size samples of the internal ringbuffer, + // multiplying them with the window, + // computing the fast fourier transform, + // taking half of the symetric complex outputs, + // computing the norm of the complex outputs and + // taking the log10 + stft.compute_column(&mut spectrogram_column[..]); + + // here's where you would do something with the + // spectrogram_column... + + // drop step_size samples from the internal ringbuffer of the stft + // making a step of size step_size + stft.move_to_next_column(); } } ``` From f22c999c7591c6d08a0be7d59b9b7afe86b407b6 Mon Sep 17 00:00:00 2001 From: Markus Mayer Date: Wed, 13 Oct 2021 21:04:37 +0200 Subject: [PATCH 09/13] Fix "len without is_empty" clippy warning --- src/lib.rs | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 113b8b0..0ac8d8a 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -249,6 +249,11 @@ where self.sample_ring.len() } + #[inline] + pub fn is_empty(&self) -> bool { + self.len() == 0 + } + pub fn append_samples(&mut self, input: &[T]) { self.sample_ring.push_many_back(input); } @@ -274,7 +279,7 @@ where // copy windowed real_input as real parts into complex_input // only copy `window_size` size, leave the rest in `complex_input` be zero for (src, dst) in self.real_input.iter().zip(self.complex_input.iter_mut()) { - dst.re = src.clone(); + dst.re = *src; } // compute fft @@ -290,7 +295,7 @@ where self.compute_into_complex_output(); for (dst, src) in output.iter_mut().zip(self.complex_output.iter()) { - *dst = src.clone(); + *dst = *src; } } From 1ec2791a9881bb7a9b96faf852c308f63b472d13 Mon Sep 17 00:00:00 2001 From: Markus Mayer Date: Wed, 20 Oct 2021 19:06:30 +0200 Subject: [PATCH 10/13] Fix wrong output_size, add tests for padded STFT --- src/lib.rs | 2 +- tests/lib.rs | 21 +++++++++++++++++++++ 2 files changed, 22 insertions(+), 1 deletion(-) diff --git a/src/lib.rs b/src/lib.rs index 0ac8d8a..5dde823 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -241,7 +241,7 @@ where #[inline] pub fn output_size(&self) -> usize { - self.window_size / 2 + self.fft_size / 2 } #[inline] diff --git a/tests/lib.rs b/tests/lib.rs index ea5f393..f2cfe95 100644 --- a/tests/lib.rs +++ b/tests/lib.rs @@ -64,3 +64,24 @@ fn test_stft() { stft.compute_column(&mut output[..]); println!("{:?}", output); } + +#[test] +fn test_stft_padded() { + let mut stft = STFT::new_with_zero_padding(WindowType::Hanning, 8, 32, 4); + assert!(!stft.contains_enough_to_compute()); + assert_eq!(stft.output_size(), 16); + assert_eq!(stft.len(), 0); + stft.append_samples(&vec![500., 0., 100.][..]); + assert_eq!(stft.len(), 3); + assert!(!stft.contains_enough_to_compute()); + stft.append_samples(&vec![500., 0., 100., 0.][..]); + assert_eq!(stft.len(), 7); + assert!(!stft.contains_enough_to_compute()); + + stft.append_samples(&vec![500.][..]); + assert!(stft.contains_enough_to_compute()); + + let mut output: Vec = vec![0.; 16]; + stft.compute_column(&mut output[..]); + println!("{:?}", output); +} From 932c96632875abdb4244ed2439238aeb2d46a256 Mon Sep 17 00:00:00 2001 From: Markus Mayer Date: Wed, 20 Oct 2021 19:08:55 +0200 Subject: [PATCH 11/13] Simplify passing of vectors in tests --- tests/lib.rs | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/tests/lib.rs b/tests/lib.rs index f2cfe95..54edd67 100644 --- a/tests/lib.rs +++ b/tests/lib.rs @@ -50,18 +50,18 @@ fn test_stft() { assert!(!stft.contains_enough_to_compute()); assert_eq!(stft.output_size(), 4); assert_eq!(stft.len(), 0); - stft.append_samples(&vec![500., 0., 100.][..]); + stft.append_samples(&[500., 0., 100.]); assert_eq!(stft.len(), 3); assert!(!stft.contains_enough_to_compute()); - stft.append_samples(&vec![500., 0., 100., 0.][..]); + stft.append_samples(&[500., 0., 100., 0.]); assert_eq!(stft.len(), 7); assert!(!stft.contains_enough_to_compute()); - stft.append_samples(&vec![500.][..]); + stft.append_samples(&[500.]); assert!(stft.contains_enough_to_compute()); let mut output: Vec = vec![0.; 4]; - stft.compute_column(&mut output[..]); + stft.compute_column(&mut output); println!("{:?}", output); } @@ -71,17 +71,17 @@ fn test_stft_padded() { assert!(!stft.contains_enough_to_compute()); assert_eq!(stft.output_size(), 16); assert_eq!(stft.len(), 0); - stft.append_samples(&vec![500., 0., 100.][..]); + stft.append_samples(&[500., 0., 100.]); assert_eq!(stft.len(), 3); assert!(!stft.contains_enough_to_compute()); - stft.append_samples(&vec![500., 0., 100., 0.][..]); + stft.append_samples(&[500., 0., 100., 0.]); assert_eq!(stft.len(), 7); assert!(!stft.contains_enough_to_compute()); - stft.append_samples(&vec![500.][..]); + stft.append_samples(&[500.]); assert!(stft.contains_enough_to_compute()); let mut output: Vec = vec![0.; 16]; - stft.compute_column(&mut output[..]); + stft.compute_column(&mut output); println!("{:?}", output); } From be2be92497c5dd63eae64b7a06f04349dab561cb Mon Sep 17 00:00:00 2001 From: Markus Mayer Date: Wed, 20 Oct 2021 19:30:09 +0200 Subject: [PATCH 12/13] Reuse complex input/output buffer This also adds explicit zero-padding of the input/output buffer to ensure that previous calculations don't "spill over". --- src/lib.rs | 33 +++++++++++++++++++++------------ 1 file changed, 21 insertions(+), 12 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 5dde823..f0065cf 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -144,8 +144,7 @@ where /// internal ringbuffer used to store samples pub sample_ring: SliceRingImpl, pub real_input: Vec, - pub complex_input: Vec>, - pub complex_output: Vec>, + pub complex_input_output: Vec>, fft_scratch: Vec>, } @@ -221,13 +220,10 @@ where window, real_input: std::iter::repeat(T::zero()).take(window_size).collect(), // zero-padded complex_input, so the size is fft_size, not window_size - complex_input: std::iter::repeat(Complex::::zero()) + complex_input_output: std::iter::repeat(Complex::::zero()) .take(fft_size) .collect(), // same size as complex_output - complex_output: std::iter::repeat(Complex::::zero()) - .take(fft_size) - .collect(), } } @@ -267,7 +263,7 @@ where assert!(self.contains_enough_to_compute()); // read into real_input - self.sample_ring.read_many_front(&mut self.real_input[..]); + self.sample_ring.read_many_front(&mut self.real_input); // multiply real_input with window if let Some(ref window) = self.window { @@ -278,13 +274,26 @@ where // copy windowed real_input as real parts into complex_input // only copy `window_size` size, leave the rest in `complex_input` be zero - for (src, dst) in self.real_input.iter().zip(self.complex_input.iter_mut()) { + for (src, dst) in self + .real_input + .iter() + .zip(self.complex_input_output.iter_mut()) + { dst.re = *src; + dst.im = T::zero(); + } + + // ensure the buffer is indeed zero-padded when needed. + if self.window_size < self.fft_size { + for dst in self.complex_input_output.iter_mut().skip(self.window_size) { + dst.re = T::zero(); + dst.im = T::zero(); + } } // compute fft self.fft - .process_with_scratch(&mut self.complex_output, &mut self.fft_scratch) + .process_with_scratch(&mut self.complex_input_output, &mut self.fft_scratch) } /// # Panics @@ -294,7 +303,7 @@ where self.compute_into_complex_output(); - for (dst, src) in output.iter_mut().zip(self.complex_output.iter()) { + for (dst, src) in output.iter_mut().zip(self.complex_input_output.iter()) { *dst = *src; } } @@ -306,7 +315,7 @@ where self.compute_into_complex_output(); - for (dst, src) in output.iter_mut().zip(self.complex_output.iter()) { + for (dst, src) in output.iter_mut().zip(self.complex_input_output.iter()) { *dst = src.norm(); } } @@ -319,7 +328,7 @@ where self.compute_into_complex_output(); - for (dst, src) in output.iter_mut().zip(self.complex_output.iter()) { + for (dst, src) in output.iter_mut().zip(self.complex_input_output.iter()) { *dst = log10_positive(src.norm()); } } From bc091034cf968cab5f3145ff457a5d7b9873d049 Mon Sep 17 00:00:00 2001 From: Markus Mayer Date: Wed, 20 Oct 2021 19:30:25 +0200 Subject: [PATCH 13/13] Add missing result assertion to integration tests In addition, the FFT calculation is now performed twice to ensure that the shared input/output buffer does not affect future calculations. --- Cargo.toml | 3 ++- tests/lib.rs | 39 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 41 insertions(+), 1 deletion(-) diff --git a/Cargo.toml b/Cargo.toml index f800952..bca2040 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -8,7 +8,7 @@ license = "MIT OR Apache-2.0" name = "stft" readme = "README.md" repository = "https://github.com/snd/stft.git" -version = "0.3.0" +version = "0.2.0" edition = "2018" [dependencies] @@ -19,6 +19,7 @@ strider = "0.1.3" [dev-dependencies] criterion = "0.3.4" +approx = "0.5.0" [[bench]] name = "lib" diff --git a/tests/lib.rs b/tests/lib.rs index 54edd67..ff0a3ec 100644 --- a/tests/lib.rs +++ b/tests/lib.rs @@ -1,3 +1,4 @@ +use approx::assert_ulps_eq; use std::str::FromStr; use stft::{WindowType, STFT}; @@ -63,6 +64,19 @@ fn test_stft() { let mut output: Vec = vec![0.; 4]; stft.compute_column(&mut output); println!("{:?}", output); + + let expected = vec![ + 2.7763337740785166, + 2.7149781042402594, + 2.6218024907053796, + 2.647816050270838, + ]; + assert_ulps_eq!(output.as_slice(), expected.as_slice(), max_ulps = 10); + + // repeat the calculation to ensure results are independent of the internal buffer + let mut output2: Vec = vec![0.; 4]; + stft.compute_column(&mut output2); + assert_ulps_eq!(output.as_slice(), output2.as_slice(), max_ulps = 10); } #[test] @@ -84,4 +98,29 @@ fn test_stft_padded() { let mut output: Vec = vec![0.; 16]; stft.compute_column(&mut output); println!("{:?}", output); + + let expected = vec![ + 2.7763337740785166, + 2.772158781619449, + 2.7598791705720664, + 2.740299218211912, + 2.7149781042402594, + 2.686495897766628, + 2.6585877421915676, + 2.635728083951981, + 2.6218024907053796, + 2.6183544930578027, + 2.6238833073831658, + 2.634925941918913, + 2.647816050270838, + 2.65977332745612, + 2.6691025866822033, + 2.6749381613735683, + ]; + assert_ulps_eq!(output.as_slice(), expected.as_slice(), max_ulps = 10); + + // repeat the calculation to ensure results are independent of the internal buffer + let mut output2: Vec = vec![0.; 16]; + stft.compute_column(&mut output2); + assert_ulps_eq!(output.as_slice(), output2.as_slice(), max_ulps = 10); }