dft/
lib.rs

1//! [Discrete Fourier transform][1].
2//!
3//! The `Transform` trait is responsible for performing the transform. The trait
4//! is implemented for both real and complex data. There are three transform
5//! operations available: forward, backward, and inverse. The desired operation
6//! is specified by the `Operation` enumeration passed to the `Plan::new`
7//! function, which precomputes auxiliary information needed for
8//! `Transform::transform`. All the operations are preformed in place.
9//!
10//! When applied to real data, the transform works as follows. If the operation
11//! is forward, the data are replaced by the positive frequency half of their
12//! complex transform. The first and last components of the complex transform,
13//! which are real, are stored in `self[0]` and `self[1]`, respectively. If the
14//! operation is backward or inverse, the data are assumed to be stored
15//! according to the above convention. See the reference below for further
16//! details.
17//!
18//! ## Example
19//!
20//! ```
21//! use dft::{Operation, Plan, c64};
22//!
23//! let plan = Plan::new(Operation::Forward, 512);
24//! let mut data = vec![c64::new(42.0, 69.0); 512];
25//! dft::transform(&mut data, &plan);
26//! ```
27//!
28//! ## References
29//!
30//! 1. W. Press, S. Teukolsky, W. Vetterling, and B. Flannery, “Numerical
31//! Recipes 3rd Edition: The Art of Scientific Computing,” Cambridge University
32//! Press, 2007.
33//!
34//! [1]: https://en.wikipedia.org/wiki/Discrete_Fourier_transform
35
36extern crate num_complex;
37extern crate num_traits;
38
39use num_complex::Complex;
40use num_traits::{Float, FloatConst, One};
41
42/// A complex number with 32-bit parts.
43#[allow(non_camel_case_types)]
44pub type c32 = Complex<f32>;
45
46/// A complex number with 64-bit parts.
47#[allow(non_camel_case_types)]
48pub type c64 = Complex<f64>;
49
50mod complex;
51mod real;
52
53pub use real::unpack;
54
55/// A transform operation.
56#[derive(Clone, Copy, Debug, Eq, PartialEq)]
57pub enum Operation {
58    /// The forward transform.
59    Forward,
60    /// The backward transform.
61    Backward,
62    /// The inverse transform.
63    Inverse,
64}
65
66/// A transform plan.
67#[derive(Clone, Debug)]
68pub struct Plan<T> {
69    n: usize,
70    factors: Vec<Complex<T>>,
71    operation: Operation,
72}
73
74/// The transform.
75pub trait Transform<T> {
76    /// Perform the transform.
77    fn transform(&mut self, &Plan<T>);
78}
79
80impl<T> Plan<T>
81where
82    T: Float + FloatConst,
83{
84    /// Create a plan for a specific operation and specific number of points.
85    ///
86    /// The number of points should be a power of two.
87    pub fn new(operation: Operation, n: usize) -> Self {
88        assert!(n.is_power_of_two());
89        let one = T::one();
90        let two = one + one;
91        let mut factors = vec![];
92        let sign = if let Operation::Forward = operation {
93            -one
94        } else {
95            one
96        };
97        let mut step = 1;
98        while step < n {
99            let (multiplier, mut factor) = {
100                let theta = T::PI() / T::from(step).unwrap();
101                let sine = (theta / two).sin();
102                (
103                    Complex::new(-two * sine * sine, sign * theta.sin()),
104                    Complex::one(),
105                )
106            };
107            for _ in 0..step {
108                factors.push(factor);
109                factor = multiplier * factor + factor;
110            }
111            step <<= 1;
112        }
113        Plan {
114            n: n,
115            factors: factors,
116            operation: operation,
117        }
118    }
119}
120
121/// Perform the transform.
122///
123/// The function is a shortcut for `Transform::transform`.
124#[inline(always)]
125pub fn transform<D: ?Sized, T>(data: &mut D, plan: &Plan<T>)
126where
127    D: Transform<T>,
128{
129    Transform::transform(data, plan);
130}
dft/lib.rs

dft/
lib.rs
