From 6f7dac7af63fe37187aa1f5a48815dd50412d196 Mon Sep 17 00:00:00 2001 From: Ivan Ukhov Date: Thu, 25 Jun 2015 10:57:58 -0400 Subject: [PATCH] Split the code into two modules --- src/complex.rs | 78 ++++++++++++++++++++++ src/lib.rs | 175 ++----------------------------------------------- src/real.rs | 91 +++++++++++++++++++++++++ tests/lib.rs | 37 ++++++----- 4 files changed, 193 insertions(+), 188 deletions(-) create mode 100644 src/complex.rs create mode 100644 src/real.rs diff --git a/src/complex.rs b/src/complex.rs new file mode 100644 index 0000000..c125ec5 --- /dev/null +++ b/src/complex.rs @@ -0,0 +1,78 @@ +//! Transformation of complex data. + +// The implementation is based on: +// http://www.librow.com/articles/article-10 + +use number::c64; + +/// Perform the forward transform. +/// +/// The number of points should be a power of two. +pub fn forward(data: &mut [c64]) { + let n = data.len(); + power_of_two!(n); + rearrange(data, n); + perform(data, n, false); +} + +/// Perform the inverse transform. +/// +/// The number of points should be a power of two. +pub fn inverse(data: &mut [c64], scaling: bool) { + let n = data.len(); + power_of_two!(n); + rearrange(data, n); + perform(data, n, true); + if scaling { + scale(data, n); + } +} + +fn rearrange(data: &mut [c64], n: usize) { + let mut target = 0; + for position in 0..n { + if target > position { + data.swap(position, target); + } + let mut mask = n >> 1; + while target & mask != 0 { + target &= !mask; + mask >>= 1; + } + target |= mask; + } +} + +fn perform(data: &mut [c64], n: usize, inverse: bool) { + use std::f64::consts::PI; + + let pi = if inverse { PI } else { -PI }; + let mut step = 1; + while step < n { + let jump = step << 1; + let (multiplier, mut factor) = { + let delta = pi / step as f64; + let sine = (0.5 * delta).sin(); + (c64(-2.0 * sine * sine, delta.sin()), c64(1.0, 0.0)) + }; + for group in 0..step { + let mut pair = group; + while pair < n { + let match_pair = pair + step; + let product = factor * data[match_pair]; + data[match_pair] = data[pair] - product; + data[pair] = data[pair] + product; + pair += jump; + } + factor = multiplier * factor + factor; + } + step <<= 1; + } +} + +fn scale(data: &mut [c64], n: usize) { + let factor = 1.0 / n as f64; + for position in 0..n { + data[position] = data[position] * factor; + } +} diff --git a/src/lib.rs b/src/lib.rs index 1db875d..caca504 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -4,12 +4,7 @@ //! [1]: https://en.wikipedia.org/wiki/Fast_Fourier_transform //! [2]: https://en.wikipedia.org/wiki/Discrete_Fourier_transform -// The implementation is based on: -// http://www.librow.com/articles/article-10 - -extern crate complex; - -pub use complex::c64; +extern crate complex as number; macro_rules! power_of_two( ($n:expr) => (if $n < 1 || $n & ($n - 1) != 0 { @@ -17,169 +12,7 @@ macro_rules! power_of_two( }); ); -/// Perform the forward transform. -/// -/// The number of points should be a power of two. -pub fn forward(data: &mut [c64]) { - let n = data.len(); - power_of_two!(n); - rearrange(data, n); - perform(data, n, false); -} - -/// Perform the forward transform of real data. -/// -/// The number of points should be a power of two. The data are replaced by the -/// positive frequency half of their complex Fourier transform. The real-valued -/// first and last components of the complex transform are returned as elements -/// data[0] and data[1], respectively. -/// -/// ## References -/// -/// 1. William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. -/// Flannery, “Numerical Recipes 3rd Edition: The Art of Scientific -/// Computing,” Cambridge University Press, 2007. -pub fn forward_real(data: &mut [f64]) { - use std::f64::consts::PI; - use std::slice::from_raw_parts_mut; - - let n = data.len(); - power_of_two!(n); - - forward(unsafe { - from_raw_parts_mut(data.as_mut_ptr() as *mut _, n / 2) - }); - - let (c1, c2) = (0.5, -0.5); - let (wpr, wpi) = { - let theta = -PI / (n / 2) as f64; - let temp = (0.5 * theta).sin(); - (-2.0 * temp * temp, theta.sin()) - }; - let mut wr = 1.0 + wpr; - let mut wi = wpi; - for i in 1..(n / 2 / 2) { - let i1 = i + i; - let i2 = 1 + i1; - let i3 = n - i1; - let i4 = 1 + i3; - - let h1r = c1 * (data[i1] + data[i3]); - let h1i = c1 * (data[i2] - data[i4]); - let h2r = -c2 * (data[i2] + data[i4]); - let h2i = c2 * (data[i1] - data[i3]); - - data[i1] = h1r + wr * h2r - wi * h2i; - data[i2] = h1i + wr * h2i + wi * h2r; - data[i3] = h1r - wr * h2r + wi * h2i; - data[i4] = -h1i + wr * h2i + wi * h2r; - - let temp = wr; - wr = wr * wpr - wi * wpi + wr; - wi = wi * wpr + temp * wpi + wi; - } - data[n / 2 + 1] *= -1.0; - let temp = data[0]; - data[0] = temp + data[1]; - data[1] = temp - data[1]; -} - -/// Perform the inverse transform. -/// -/// The number of points should be a power of two. -pub fn inverse(data: &mut [c64], scaling: bool) { - let n = data.len(); - power_of_two!(n); - rearrange(data, n); - perform(data, n, true); - if scaling { - scale(data, n); - } -} - -/// Expand a compressed representation produced by `forward_real`. -pub fn unpack_real(data: &[f64]) -> Vec { - let n = data.len(); - power_of_two!(n); - - let mut cdata = Vec::with_capacity(n); - unsafe { cdata.set_len(n) }; - - cdata[0] = c64(data[0], 0.0); - for i in 1..(n / 2) { - cdata[i] = c64(data[2 * i], data[2 * i + 1]); - } - cdata[n / 2] = c64(data[1], 0.0); - for i in (n / 2 + 1)..n { - cdata[i] = c64(data[2 * (n - i)], -data[2 * (n - i) + 1]); - } - - cdata -} - -fn rearrange(data: &mut [c64], n: usize) { - let mut target = 0; - for position in 0..n { - if target > position { - data.swap(position, target); - } - let mut mask = n >> 1; - while target & mask != 0 { - target &= !mask; - mask >>= 1; - } - target |= mask; - } -} - -fn perform(data: &mut [c64], n: usize, inverse: bool) { - use std::f64::consts::PI; - - let pi = if inverse { PI } else { -PI }; - let mut step = 1; - while step < n { - let jump = step << 1; - let delta = pi / step as f64; - let sine = (0.5 * delta).sin(); - let multiplier = c64(-2.0 * sine * sine, delta.sin()); - let mut factor = c64(1.0, 0.0); - for group in 0..step { - let mut pair = group; - while pair < n { - let match_pair = pair + step; - let product = factor * data[match_pair]; - data[match_pair] = data[pair] - product; - data[pair] = data[pair] + product; - pair += jump; - } - factor = multiplier * factor + factor; - } - step <<= 1; - } -} - -fn scale(data: &mut [c64], n: usize) { - let factor = 1.0 / n as f64; - for position in 0..n { - data[position] = data[position] * factor; - } -} - -#[cfg(test)] -mod tests { - use c64; - - #[test] - fn unpack_real() { - let data = (0..4).map(|i| (i + 1) as f64).collect::>(); - assert_eq!(::unpack_real(&data), vec![ - c64(1.0, 0.0), c64(3.0, 4.0), c64(2.0, 0.0), c64(3.0, -4.0), - ]); +pub mod complex; +pub mod real; - let data = (0..8).map(|i| (i + 1) as f64).collect::>(); - assert_eq!(::unpack_real(&data), vec![ - c64(1.0, 0.0), c64(3.0, 4.0), c64(5.0, 6.0), c64(7.0, 8.0), - c64(2.0, 0.0), c64(7.0, -8.0), c64(5.0, -6.0), c64(3.0, -4.0), - ]); - } -} +pub use complex::{forward, inverse}; diff --git a/src/real.rs b/src/real.rs new file mode 100644 index 0000000..921bfeb --- /dev/null +++ b/src/real.rs @@ -0,0 +1,91 @@ +//! Transformation of real data. + +use number::{Complex, c64}; + +/// Perform the forward transform. +/// +/// The number of points should be a power of two. The data are replaced by the +/// positive frequency half of their complex Fourier transform. The real-valued +/// first and last components of the complex transform are returned as elements +/// data[0] and data[1], respectively. +/// +/// ## References +/// +/// 1. William H. Press, Saul A. Teukolsky, William T. Vetterling, Brian P. +/// Flannery, “Numerical Recipes 3rd Edition: The Art of Scientific +/// Computing,” Cambridge University Press, 2007. +pub fn forward(data: &mut [f64]) { + use std::f64::consts::PI; + use std::slice::from_raw_parts_mut; + + let (data, n) = unsafe { + let n = data.len(); + power_of_two!(n); + (from_raw_parts_mut(data.as_mut_ptr() as *mut _, n / 2), n / 2) + }; + + ::complex::forward(data); + + let (multiplier, mut factor) = { + let delta = -PI / n as f64; + let sine = (0.5 * delta).sin(); + (c64(-2.0 * sine * sine, delta.sin()), c64(0.0, 1.0)) + }; + + data[0] = c64( + data[0].re() + data[0].im(), + data[0].re() - data[0].im(), + ); + + for i in 1..(n / 2) { + let j = n - i; + + factor = multiplier * factor + factor; + let part1 = (data[i] + data[j].conj()) * 0.5; + let part2 = -(data[i] - data[j].conj()) * 0.5; + + data[i] = part1 + factor * part2; + data[j] = (part1 - factor * part2).conj(); + } + data[n / 2] = data[n / 2].conj(); +} + +/// Expand a compressed representation produced by `real::forward`. +pub fn unpack(data: &[f64]) -> Vec { + let n = data.len(); + power_of_two!(n); + + let mut cdata = Vec::with_capacity(n); + unsafe { cdata.set_len(n) }; + + cdata[0] = c64(data[0], 0.0); + for i in 1..(n / 2) { + cdata[i] = c64(data[2 * i], data[2 * i + 1]); + } + cdata[n / 2] = c64(data[1], 0.0); + for i in (n / 2 + 1)..n { + let j = n - i; + cdata[i] = c64(data[2 * j], -data[2 * j + 1]); + } + + cdata +} + +#[cfg(test)] +mod tests { + use number::c64; + + #[test] + fn unpack() { + let data = (0..4).map(|i| (i + 1) as f64).collect::>(); + assert_eq!(super::unpack(&data), vec![ + c64(1.0, 0.0), c64(3.0, 4.0), c64(2.0, 0.0), c64(3.0, -4.0), + ]); + + let data = (0..8).map(|i| (i + 1) as f64).collect::>(); + assert_eq!(super::unpack(&data), vec![ + c64(1.0, 0.0), c64(3.0, 4.0), c64(5.0, 6.0), c64(7.0, 8.0), + c64(2.0, 0.0), c64(7.0, -8.0), c64(5.0, -6.0), c64(3.0, -4.0), + ]); + } +} diff --git a/tests/lib.rs b/tests/lib.rs index 8754697..bad40ed 100644 --- a/tests/lib.rs +++ b/tests/lib.rs @@ -1,49 +1,52 @@ extern crate assert; +extern crate complex; extern crate fft; +use complex::c64; + mod fixtures; #[test] -fn forward() { +fn complex_forward() { let mut data = fixtures::TIME_DATA.to_vec(); - fft::forward(as_c64_mut(&mut data)); + fft::complex::forward(as_c64_mut(&mut data)); assert::close(&data, &fixtures::FREQUENCY_DATA_FROM_COMPLEX[..], 1e-14); } #[test] -fn forward_real() { +fn complex_inverse() { + let mut data = fixtures::FREQUENCY_DATA_FROM_COMPLEX.to_vec(); + fft::complex::inverse(as_c64_mut(&mut data), true); + assert::close(&data, &fixtures::TIME_DATA[..], 1e-14); +} + +#[test] +fn real_forward() { let mut data = fixtures::TIME_DATA.to_vec(); { let mut data = to_c64(&data); - fft::forward(&mut data); + fft::complex::forward(&mut data); assert::close(as_f64(&data), &fixtures::FREQUENCY_DATA_FROM_REAL[..], 1e-13); } { - fft::forward_real(&mut data); - let data = fft::unpack_real(&data); + fft::real::forward(&mut data); + let data = fft::real::unpack(&data); assert::close(as_f64(&data), &fixtures::FREQUENCY_DATA_FROM_REAL[..], 1e-13); } } -#[test] -fn inverse() { - let mut data = fixtures::FREQUENCY_DATA_FROM_COMPLEX.to_vec(); - fft::inverse(as_c64_mut(&mut data), true); - assert::close(&data, &fixtures::TIME_DATA[..], 1e-14); -} - -fn as_f64<'l>(slice: &'l [fft::c64]) -> &'l [f64] { +fn as_f64<'l>(slice: &'l [c64]) -> &'l [f64] { unsafe { std::slice::from_raw_parts(slice.as_ptr() as *const _, 2 * slice.len()) } } -fn as_c64_mut<'l>(slice: &'l mut [f64]) -> &'l mut [fft::c64] { +fn as_c64_mut<'l>(slice: &'l mut [f64]) -> &'l mut [c64] { unsafe { std::slice::from_raw_parts_mut(slice.as_mut_ptr() as *mut _, slice.len() / 2) } } -fn to_c64(slice: &[f64]) -> Vec { - slice.iter().map(|&re| fft::c64(re, 0.0)).collect() +fn to_c64(slice: &[f64]) -> Vec { + slice.iter().map(|&re| c64(re, 0.0)).collect() }