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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130

//! [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);
}