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
//! [Discrete Fourier transform][1]. //! //! The `Transform` trait is responsible for performing the transform. The trait //! is implemented for real and complex data, which are represented by `[f64]` //! and `[c64]`, respectively. There are three transformation 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`. //! //! ## Example //! //! ``` //! use dft::{Operation, Plan, c64}; //! //! let size = 512; //! let plan = Plan::new(Operation::Forward, size); //! let mut data = vec![c64::new(42.0, 69.0); size]; //! //! dft::transform(&mut data, &plan); //! ``` //! //! ## Real Data //! //! When applied to real data, `Transform::transform` works as follows. If the //! operation is `Operation::Forward`, 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 `self[0]` //! and `self[1]`, respectively. If the operation is `Operation::Backward` or //! `Operation::Inverse`, the function assumes that the data are packed in the //! format that has just been described. See the reference below for further //! information on the format. //! //! ## 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; /// A complex number with 64-bit parts. #[allow(non_camel_case_types)] pub type c64 = num::Complex<f64>; macro_rules! c64(($re:expr, $im:expr) => (::c64::new($re, $im))); mod complex; mod real; pub use real::unpack; /// A transformation operation. #[derive(Clone, Copy, Debug, Eq, PartialEq)] pub enum Operation { /// The forward transform. Forward, /// The backward transform. Backward, /// The inverse transform. Inverse, } /// A transformation plan. #[derive(Clone, Debug)] pub struct Plan { size: usize, factors: Vec<c64>, operation: Operation, } /// A type suitable for transformation. pub trait Transform { /// Perform the transform. fn transform(&mut self, &Plan); } impl Plan { /// Create a plan for a specific operation and number of points. /// /// The number of points should be a power of two. pub fn new(operation: Operation, size: usize) -> Plan { use std::f64::consts::PI; assert!(size.is_power_of_two(), "the number of points should be a power of two"); let mut factors = vec![]; let sign = if let Operation::Forward = operation { -1.0 } else { 1.0 }; let mut step = 1; while step < size { let (multiplier, mut factor) = { let theta = PI / step as f64; let sine = (0.5 * theta).sin(); (c64!(-2.0 * sine * sine, sign * theta.sin()), c64!(1.0, 0.0)) }; for _ in 0..step { factors.push(factor); factor = multiplier * factor + factor; } step <<= 1; } Plan { size: size, factors: factors, operation: operation } } } /// Perform the transform. /// /// The function is the same as `Transform::transform`. #[inline(always)] pub fn transform<T: Transform + ?Sized>(data: &mut T, plan: &Plan) { Transform::transform(data, plan); }