From 51fc59782963ff66e2f75f29589fa1c8bfe5f150 Mon Sep 17 00:00:00 2001 From: "Sergey \"Shnatsel\" Davidoff" Date: Mon, 29 Apr 2024 22:01:49 +0100 Subject: [PATCH 1/7] advertise runtime feature selection --- README.md | 1 + 1 file changed, 1 insertion(+) diff --git a/README.md b/README.md index e885f3b..54c505a 100644 --- a/README.md +++ b/README.md @@ -15,6 +15,7 @@ Transform (FFT) library written in pure Rust. - Performance on par with other Rust FFT implementations - Zero `unsafe` code - Takes advantage of latest CPU features up to and including `AVX-512`, but performs well even without them +- Selects the fastest implementation at runtime. No need for `-C target-cpu=native`! - Optional parallelization of some steps to 2 threads (with even more planned) - 2x lower memory usage than [RustFFT](https://crates.io/crates/rustfft/) - Python bindings (via [PyO3](https://github.com/PyO3/pyo3)) From e150db3a065f02cf016ba2bf2da7faaef5df2bdf Mon Sep 17 00:00:00 2001 From: Saveliy Yusufov Date: Fri, 3 May 2024 12:09:29 -0400 Subject: [PATCH 2/7] Fixes inverse FFT ouput order issue --- src/lib.rs | 51 +++++++++++++++++++++++++++++++++++++++++++++++++- src/planner.rs | 5 +++-- 2 files changed, 53 insertions(+), 3 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index c65a74f..534afa0 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -47,14 +47,35 @@ macro_rules! impl_fft_for { imags.len() ); - let mut planner = <$planner>::new(reals.len(), direction); + let mut planner = <$planner>::new(reals.len(), Direction::Forward); assert!( planner.num_twiddles().is_power_of_two() && planner.num_twiddles() == reals.len() / 2 ); let opts = Options::guess_options(reals.len()); + + match direction { + Direction::Reverse => { + for z_im in imags.iter_mut() { + *z_im = -*z_im; + } + } + _ => (), + } + $opts_and_plan(reals, imags, &opts, &mut planner); + + match direction { + Direction::Reverse => { + let scaling_factor = (reals.len() as $precision).recip(); + for (z_re, z_im) in reals.iter_mut().zip(imags.iter_mut()) { + *z_re = (*z_re) * scaling_factor; + *z_im = -(*z_im) * -scaling_factor; + } + } + _ => (), + } } }; } @@ -256,4 +277,32 @@ mod tests { test_fft_correctness!(fft_correctness_32, f32, fft_32, 4, 9); test_fft_correctness!(fft_correctness_64, f64, fft_64, 4, 17); + + #[test] + fn ifft_using_fft() { + let n = 4; + let big_n = 1 << n; + + let mut reals: Vec<_> = (1..=big_n).map(|i| i as f64).collect(); + let mut imags: Vec<_> = vec![0.0; big_n]; + println!("{:?}", reals); + println!("{:?}\n", imags); + + fft_64(&mut reals, &mut imags, Direction::Forward); + println!("{:?}", reals); + println!("{:?}\n", imags); + + fft_64(&mut reals, &mut imags, Direction::Reverse); + println!("{:?}", reals); + println!("{:?}", imags); + + let mut signal_re = 1.0; + + // Now check that the identity is indeed the original signal we generated above + for (z_re, z_im) in reals.into_iter().zip(imags.into_iter()) { + assert_float_closeness(z_re, signal_re, 1e-4); + assert_float_closeness(z_im, 0.0, 1e-4); + signal_re += 1.0; + } + } } diff --git a/src/planner.rs b/src/planner.rs index 548e972..9564b89 100644 --- a/src/planner.rs +++ b/src/planner.rs @@ -6,6 +6,7 @@ use crate::twiddles::{generate_twiddles, generate_twiddles_simd_32, generate_twi /// Reverse is for running the Inverse Fast Fourier Transform (IFFT) /// Forward is for running the regular FFT +#[derive(Copy, Clone)] pub enum Direction { /// Leave the exponent term in the twiddle factor alone Forward = 1, @@ -46,9 +47,9 @@ macro_rules! impl_planner_for { let dist = num_points >> 1; let (twiddles_re, twiddles_im) = if dist >= 8 * 2 { - $generate_twiddles_simd_fn(dist, direction) + $generate_twiddles_simd_fn(dist, Direction::Forward) } else { - generate_twiddles(dist, direction) + generate_twiddles(dist, Direction::Forward) }; assert_eq!(twiddles_re.len(), twiddles_im.len()); From 44b80a651618cccd7523bd5a9cc1360b2b53c6b9 Mon Sep 17 00:00:00 2001 From: Saveliy Yusufov Date: Fri, 3 May 2024 12:37:00 -0400 Subject: [PATCH 3/7] Move twiddles fwd * rev = identity test - Testing that all the values in the forward twiddle factors pointwise multiplied by all the values in the inverse twiddle factors is more of a unit test for twiddles. It made no sense to keep it in the planner. - The previous commit modifies the planner in a way that this test would no longer pass. We can get the same, good coverage by keeping this test under twiddles. --- src/planner.rs | 40 +--------------------------------------- src/twiddles.rs | 45 +++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 46 insertions(+), 39 deletions(-) diff --git a/src/planner.rs b/src/planner.rs index 9564b89..b174567 100644 --- a/src/planner.rs +++ b/src/planner.rs @@ -35,7 +35,7 @@ macro_rules! impl_planner_for { /// # Panics /// /// Panics if `num_points < 1` or if `num_points` is __not__ a power of 2. - pub fn new(num_points: usize, direction: Direction) -> Self { + pub fn new(num_points: usize, _direction: Direction) -> Self { assert!(num_points > 0 && num_points.is_power_of_two()); if num_points <= 4 { return Self { @@ -73,8 +73,6 @@ impl_planner_for!(Planner32, f32, generate_twiddles_simd_32); #[cfg(test)] mod tests { - use utilities::assert_float_closeness; - use super::*; macro_rules! test_no_twiddles { @@ -91,40 +89,4 @@ mod tests { test_no_twiddles!(no_twiddles_64, Planner64); test_no_twiddles!(no_twiddles_32, Planner32); - - macro_rules! forward_mul_inverse_eq_identity { - ($test_name:ident, $planner:ty) => { - #[test] - fn $test_name() { - for i in 3..25 { - let num_points = 1 << i; - let planner_forward = <$planner>::new(num_points, Direction::Forward); - let planner_reverse = <$planner>::new(num_points, Direction::Reverse); - - assert_eq!( - planner_reverse.num_twiddles(), - planner_forward.num_twiddles() - ); - - // (a + ib) (c + id) = ac + iad + ibc - bd - // = ac - bd + i(ad + bc) - planner_forward - .twiddles_re - .iter() - .zip(planner_forward.twiddles_im.iter()) - .zip(planner_reverse.twiddles_re.iter()) - .zip(planner_reverse.twiddles_im) - .for_each(|(((a, b), c), d)| { - let temp_re = a * c - b * d; - let temp_im = a * d + b * c; - assert_float_closeness(temp_re, 1.0, 1e-2); - assert_float_closeness(temp_im, 0.0, 1e-2); - }); - } - } - }; - } - - forward_mul_inverse_eq_identity!(forward_reverse_eq_identity_64, Planner64); - forward_mul_inverse_eq_identity!(forward_reverse_eq_identity_32, Planner32); } diff --git a/src/twiddles.rs b/src/twiddles.rs index b754a11..13dedba 100644 --- a/src/twiddles.rs +++ b/src/twiddles.rs @@ -300,4 +300,49 @@ mod tests { } } } + + macro_rules! forward_mul_inverse_eq_identity { + ($test_name:ident, $generate_twiddles_simd_fn:ident) => { + #[test] + fn $test_name() { + for i in 3..25 { + let num_points = 1 << i; + let dist = num_points >> 1; + + let (fwd_twiddles_re, fwd_twiddles_im) = if dist >= 8 * 2 { + $generate_twiddles_simd_fn(dist, Direction::Forward) + } else { + generate_twiddles(dist, Direction::Forward) + }; + + assert_eq!(fwd_twiddles_re.len(), fwd_twiddles_im.len()); + + let (rev_twiddles_re, rev_twiddles_im) = if dist >= 8 * 2 { + $generate_twiddles_simd_fn(dist, Direction::Reverse) + } else { + generate_twiddles(dist, Direction::Reverse) + }; + + assert_eq!(rev_twiddles_re.len(), rev_twiddles_im.len()); + + // (a + ib) (c + id) = ac + iad + ibc - bd + // = ac - bd + i(ad + bc) + fwd_twiddles_re + .iter() + .zip(fwd_twiddles_im.iter()) + .zip(rev_twiddles_re.iter()) + .zip(rev_twiddles_im.iter()) + .for_each(|(((a, b), c), d)| { + let temp_re = a * c - b * d; + let temp_im = a * d + b * c; + assert_float_closeness(temp_re, 1.0, 1e-2); + assert_float_closeness(temp_im, 0.0, 1e-2); + }); + } + } + }; + } + + forward_mul_inverse_eq_identity!(forward_reverse_eq_identity_64, generate_twiddles_simd_64); + forward_mul_inverse_eq_identity!(forward_reverse_eq_identity_32, generate_twiddles_simd_32); } From aaa6a89bf1c76fbdebc8a6d5c63d883947bd6fbd Mon Sep 17 00:00:00 2001 From: Saveliy Yusufov Date: Fri, 3 May 2024 12:50:16 -0400 Subject: [PATCH 4/7] Placate clippy's complaints about assign op --- src/lib.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 534afa0..f729925 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -70,8 +70,8 @@ macro_rules! impl_fft_for { Direction::Reverse => { let scaling_factor = (reals.len() as $precision).recip(); for (z_re, z_im) in reals.iter_mut().zip(imags.iter_mut()) { - *z_re = (*z_re) * scaling_factor; - *z_im = -(*z_im) * -scaling_factor; + *z_re *= scaling_factor; + *z_im *= -scaling_factor; } } _ => (), From 899ce9cc4e70c1894b4a14049e87bf0700d45b7c Mon Sep 17 00:00:00 2001 From: Saveliy Yusufov Date: Fri, 3 May 2024 14:37:31 -0400 Subject: [PATCH 5/7] Bump version to 0.2.1 --- Cargo.toml | 2 +- pyphastft/Cargo.toml | 2 +- utilities/Cargo.toml | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index 98e51eb..12d50ef 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "phastft" -version = "0.2.0" +version = "0.2.1" edition = "2021" authors = ["Saveliy Yusufov", "Shnatsel"] license = "MIT OR Apache-2.0" diff --git a/pyphastft/Cargo.toml b/pyphastft/Cargo.toml index eaf4221..a1d234c 100644 --- a/pyphastft/Cargo.toml +++ b/pyphastft/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "pyphastft" -version = "0.2.0" +version = "0.2.1" edition = "2021" authors = ["Saveliy Yusufov", "Shnatsel"] license = "MIT OR Apache-2.0" diff --git a/utilities/Cargo.toml b/utilities/Cargo.toml index e4faada..4bd855b 100644 --- a/utilities/Cargo.toml +++ b/utilities/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "utilities" -version = "0.1.0" +version = "0.2.1" edition = "2021" # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html From d3254ffaad75f40ac77d935801c558d7d8b6177d Mon Sep 17 00:00:00 2001 From: Saveliy Yusufov Date: Fri, 10 May 2024 11:35:30 -0400 Subject: [PATCH 6/7] Add a more robust round-trip FFT test --- src/lib.rs | 54 ++++++++++++++++++++++++++++-------------------------- 1 file changed, 28 insertions(+), 26 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index f729925..18c7291 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -176,9 +176,9 @@ impl_fft_with_opts_and_plan_for!( mod tests { use std::ops::Range; - use utilities::assert_float_closeness; - use utilities::rustfft::num_complex::Complex; + use utilities::{assert_float_closeness, gen_random_signal}; use utilities::rustfft::FftPlanner; + use utilities::rustfft::num_complex::Complex; use super::*; @@ -279,30 +279,32 @@ mod tests { test_fft_correctness!(fft_correctness_64, f64, fft_64, 4, 17); #[test] - fn ifft_using_fft() { - let n = 4; - let big_n = 1 << n; - - let mut reals: Vec<_> = (1..=big_n).map(|i| i as f64).collect(); - let mut imags: Vec<_> = vec![0.0; big_n]; - println!("{:?}", reals); - println!("{:?}\n", imags); - - fft_64(&mut reals, &mut imags, Direction::Forward); - println!("{:?}", reals); - println!("{:?}\n", imags); - - fft_64(&mut reals, &mut imags, Direction::Reverse); - println!("{:?}", reals); - println!("{:?}", imags); - - let mut signal_re = 1.0; - - // Now check that the identity is indeed the original signal we generated above - for (z_re, z_im) in reals.into_iter().zip(imags.into_iter()) { - assert_float_closeness(z_re, signal_re, 1e-4); - assert_float_closeness(z_im, 0.0, 1e-4); - signal_re += 1.0; + fn fft_round_trip() { + for i in 4..23 { + let big_n = 1 << i; + let mut reals = vec![0.0; big_n]; + let mut imags = vec![0.0; big_n]; + + gen_random_signal(&mut reals, &mut imags); + + let original_reals = reals.clone(); + let original_imags = imags.clone(); + + // Forward FFT + fft_64(&mut reals, &mut imags, Direction::Forward); + + // Inverse FFT + fft_64(&mut reals, &mut imags, Direction::Reverse); + + // Ensure we get back the original signal within some tolerance + for ((orig_re, orig_im), (res_re, res_im)) in original_reals + .into_iter() + .zip(original_imags.into_iter()) + .zip(reals.into_iter().zip(imags.into_iter())) + { + assert_float_closeness(res_re, orig_re, 1e-6); + assert_float_closeness(res_im, orig_im, 1e-6); + } } } } From 7e8176a128e706a759ae74e92d11801908d3610d Mon Sep 17 00:00:00 2001 From: Saveliy Yusufov Date: Fri, 10 May 2024 11:36:04 -0400 Subject: [PATCH 7/7] Fix formatting --- src/lib.rs | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/lib.rs b/src/lib.rs index 18c7291..fd8e4b9 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -176,9 +176,9 @@ impl_fft_with_opts_and_plan_for!( mod tests { use std::ops::Range; - use utilities::{assert_float_closeness, gen_random_signal}; - use utilities::rustfft::FftPlanner; use utilities::rustfft::num_complex::Complex; + use utilities::rustfft::FftPlanner; + use utilities::{assert_float_closeness, gen_random_signal}; use super::*;