1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113

//! [Discrete Fourier transform][1].
//!
//! The `Transform` trait is responsible for performing the transform. The trait
//! is implemented for both real and complex data. There are three transform
//! operations available: forward, backward, and inverse. The desired operation
//! is specified by the `Operation` enumeration passed to the `Plan::new`
//! function, which precomputes auxiliary information needed for
//! `Transform::transform`. All the operations are preformed in place.
//!
//! When applied to real data, the transform works as follows. If the operation
//! is forward, the data are replaced by the positive frequency half of their
//! complex transform. The first and last components of the complex transform,
//! which are real, are stored in `self[0]` and `self[1]`, respectively. If the
//! operation is backward or inverse, the data are assumed to be stored
//! according to the above convention. See the reference below for further
//! details.
//!
//! ## Example
//!
//! ```
//! use dft::{Operation, Plan, c64};
//!
//! let plan = Plan::new(Operation::Forward, 512);
//! let mut data = vec![c64::new(42.0, 69.0); 512];
//! dft::transform(&mut data, &plan);
//! ```
//!
//! ## References
//!
//! 1. W. Press, S. Teukolsky, W. Vetterling, and B. Flannery, “Numerical
//! Recipes 3rd Edition: The Art of Scientific Computing,” Cambridge University
//! Press, 2007.
//!
//! [1]: https://en.wikipedia.org/wiki/Discrete_Fourier_transform

extern crate num_complex;
extern crate num_traits;

use num_complex::Complex;
use num_traits::{Float, FloatConst, One};

/// A complex number with 32-bit parts.
#[allow(non_camel_case_types)]
pub type c32 = Complex<f32>;

/// A complex number with 64-bit parts.
#[allow(non_camel_case_types)]
pub type c64 = Complex<f64>;

mod complex;
mod real;

pub use real::unpack;

/// A transform operation.
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
pub enum Operation {
    /// The forward transform.
    Forward,
    /// The backward transform.
    Backward,
    /// The inverse transform.
    Inverse,
}

/// A transform plan.
#[derive(Clone, Debug)]
pub struct Plan<T> {
    n: usize,
    factors: Vec<Complex<T>>,
    operation: Operation,
}

/// The transform.
pub trait Transform<T> {
    /// Perform the transform.
    fn transform(&mut self, &Plan<T>);
}

impl<T> Plan<T> where T: Float + FloatConst {
    /// Create a plan for a specific operation and specific number of points.
    ///
    /// The number of points should be a power of two.
    pub fn new(operation: Operation, n: usize) -> Self {
        assert!(n.is_power_of_two());
        let one = T::one();
        let two = one + one;
        let mut factors = vec![];
        let sign = if let Operation::Forward = operation { -one } else { one };
        let mut step = 1;
        while step < n {
            let (multiplier, mut factor) = {
                let theta = T::PI() / T::from(step).unwrap();
                let sine = (theta / two).sin();
                (Complex::new(-two * sine * sine, sign * theta.sin()), Complex::one())
            };
            for _ in 0..step {
                factors.push(factor);
                factor = multiplier * factor + factor;
            }
            step <<= 1;
        }
        Plan { n: n, factors: factors, operation: operation }
    }
}

/// Perform the transform.
///
/// The function is a shortcut for `Transform::transform`.
#[inline(always)]
pub fn transform<D: ?Sized, T>(data: &mut D, plan: &Plan<T>) where D: Transform<T> {
    Transform::transform(data, plan);
}