|
|
@ -9,85 +9,20 @@ |
|
|
|
|
|
|
|
extern crate complex;
|
|
|
|
|
|
|
|
use complex::{Complex, c64};
|
|
|
|
use std::slice;
|
|
|
|
|
|
|
|
/// A means of obtaining a slice of mutable complex numbers.
|
|
|
|
pub trait AsMutComplex<'l> {
|
|
|
|
type Complex: Complex + 'l;
|
|
|
|
|
|
|
|
fn as_mut_complex(self) -> &'l mut [Self::Complex];
|
|
|
|
}
|
|
|
|
|
|
|
|
macro_rules! implement(
|
|
|
|
($complex:ty) => (
|
|
|
|
impl<'l> AsMutComplex<'l> for &'l mut [$complex] {
|
|
|
|
type Complex = $complex;
|
|
|
|
|
|
|
|
#[inline(always)]
|
|
|
|
fn as_mut_complex(self) -> &'l mut [Self::Complex] {
|
|
|
|
self
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
impl<'l> AsMutComplex<'l> for &'l mut Vec<$complex> {
|
|
|
|
type Complex = $complex;
|
|
|
|
|
|
|
|
#[inline]
|
|
|
|
fn as_mut_complex(self) -> &'l mut [Self::Complex] {
|
|
|
|
<&mut [$complex]>::as_mut_complex(&mut *self)
|
|
|
|
}
|
|
|
|
}
|
|
|
|
);
|
|
|
|
);
|
|
|
|
|
|
|
|
implement!(c64);
|
|
|
|
|
|
|
|
macro_rules! implement(
|
|
|
|
($complex:ty, $real:ty) => (
|
|
|
|
impl<'l> AsMutComplex<'l> for &'l mut [$real] {
|
|
|
|
type Complex = $complex;
|
|
|
|
|
|
|
|
/// Treat the slice as a collection of pairs of real and imaginary
|
|
|
|
/// parts and reinterpret it as a slice of complex numbers.
|
|
|
|
///
|
|
|
|
/// ## Panics
|
|
|
|
///
|
|
|
|
/// The function panics if the number of elements is not even.
|
|
|
|
#[inline]
|
|
|
|
fn as_mut_complex(self) -> &'l mut [Self::Complex] {
|
|
|
|
unsafe {
|
|
|
|
let length = self.len();
|
|
|
|
assert!(length % 2 == 0, "the number of elements should be even");
|
|
|
|
slice::from_raw_parts_mut(self.as_mut_ptr() as *mut _, length / 2)
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
impl<'l> AsMutComplex<'l> for &'l mut Vec<$real> {
|
|
|
|
type Complex = $complex;
|
|
|
|
|
|
|
|
#[inline]
|
|
|
|
fn as_mut_complex(self) -> &'l mut [Self::Complex] {
|
|
|
|
<&mut [$real]>::as_mut_complex(&mut *self)
|
|
|
|
}
|
|
|
|
}
|
|
|
|
);
|
|
|
|
);
|
|
|
|
|
|
|
|
implement!(c64, f64);
|
|
|
|
pub use complex::c64;
|
|
|
|
|
|
|
|
/// Perform the Fourier transform.
|
|
|
|
///
|
|
|
|
/// The number of points should be a power of two.
|
|
|
|
pub fn forward<'l, T>(data: T) where T: AsMutComplex<'l, Complex=c64> {
|
|
|
|
let data = data.as_mut_complex();
|
|
|
|
|
|
|
|
///
|
|
|
|
/// ## Panics
|
|
|
|
///
|
|
|
|
/// The function panics if the number of points is not a power of two.
|
|
|
|
pub fn forward(data: &mut [c64]) {
|
|
|
|
let n = data.len();
|
|
|
|
if n < 1 || n & (n - 1) != 0 {
|
|
|
|
panic!("expected the number of points to be a power of two");
|
|
|
|
}
|
|
|
|
|
|
|
|
rearrange(data, n);
|
|
|
|
perform(data, n, false);
|
|
|
|
}
|
|
|
@ -95,14 +30,15 @@ pub fn forward<'l, T>(data: T) where T: AsMutComplex<'l, Complex=c64> { |
|
|
|
/// Perform the inverse Fourier transform.
|
|
|
|
///
|
|
|
|
/// The number of points should be a power of two.
|
|
|
|
pub fn inverse<'l, T>(data: T, scaling: bool) where T: AsMutComplex<'l, Complex=c64> {
|
|
|
|
let data = data.as_mut_complex();
|
|
|
|
|
|
|
|
///
|
|
|
|
/// ## Panics
|
|
|
|
///
|
|
|
|
/// The function panics if the number of points is not a power of two.
|
|
|
|
pub fn inverse(data: &mut [c64], scaling: bool) {
|
|
|
|
let n = data.len();
|
|
|
|
if n < 1 || n & (n - 1) != 0 {
|
|
|
|
panic!("expected the number of points to be a power of two");
|
|
|
|
}
|
|
|
|
|
|
|
|
rearrange(data, n);
|
|
|
|
perform(data, n, true);
|
|
|
|
if scaling {
|
|
|
@ -131,7 +67,6 @@ 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;
|
|
|
|