fftw3 0.0.2

#![crate_type="lib"]

//! High-level bindings to [FFTW3](http://fftw.org/): Fastest Fourier
//! Transform in the West.
//!
//! At the moment, this only provides one-dimensional
//! complex-to-complex transforms via `c2c_1d`, with no explicit plan
//! reuse.
//!
//! There are some modules that provide helpers to make using the
//! low-level interface manually slightly nicer and safer, in
//! particular the `lock` module assists with keeping the use of FFTW3
//! threadsafe, and the `plan` module gives a nicer interface to
//! creating and using plans.
//!
//! # Licensing
//!
//! The FFTW3 library is licensed under GPLv2 (or later) with
//! non-free/commercial licenses available from MIT via the FFTW
//! website. On the other hand, these bindings are dual-licensed
//! MIT/Apache-2.0, however, any application using FFTW3 via these
//! bindings *must* still be distributed under the terms of the GPL
//! (unless a non-free license has been arranged). The explicit
//! difference in licenses is to enable code-reuse between these
//! bindings and bindings to other FFT libraries.

#![feature(std_misc, core)]

extern crate libc;
extern crate num;

extern crate "fftw3-sys" as ffi;

use std::{error, fmt, mem};
use num::Complex;

pub mod plan;
pub mod lock;

/// Errors that can occur while planning or running a fourier transform.
#[derive(Debug)]
pub enum FftError {
    /// The input slice was too long for the the output (elements are
    /// the lengths of those respective slices).
    InsufficientSpace(usize, usize),
    /// FFTW returned an error (a null plan) when planning.
    FailedPlan,
}

impl fmt::Display for FftError {
    fn fmt(&self, fmtr: &mut fmt::Formatter) -> fmt::Result {
        match *self {
            FftError::InsufficientSpace(in_, out) => {
                write!(fmtr, "insufficient output space: given {}, needed {}", out, in_)
            }
            FftError::FailedPlan => write!(fmtr, "planning failed"),
        }
    }
}

impl error::Error for FftError {
    fn description(&self) -> &str {
        match *self {
            FftError::InsufficientSpace(_, _) => "insufficient output space",
            FftError::FailedPlan => "planning failed"
        }
    }
}

/// Perform a complex-to-complex one-dimensional Fourier transform of
/// `in_`, writing to `out`.
///
/// The direction is controlled by `forward`: `true` for the forward
/// transform (that is, multiplying by `exp(-2πi/N)`) and `false` for
/// the inverse (`exp(2πi/N)`). This uses the ESTIMATE rigor for
/// making a plan to avoid overwriting `in_` and `out` while planning,
/// hence this definitely does not achieve the peak performance
/// possible with FFTW.
///
/// An error is returned if `out` is shorter than `in_`, since there
/// is insufficient space to perform the transform, or if the FFTW3
/// planning function does not succeed.
///
/// # Example
///
/// ```rust
/// extern crate fftw3;
/// extern crate num;
/// use num::Complex;
///
/// fn main() {
///     let input = [Complex::new(2.0, 0.0), Complex::new(1.0, 1.0),
///                  Complex::new(0.0, 3.0), Complex::new(2.0, 4.0)];
///     let mut output = [Complex::new(0.0, 0.0); 4];
///
///     fftw3::c2c_1d(&input, &mut output, true).unwrap(); // panic on error
///
///     println!("the transform of {:?} is {:?}", &input[..], &output[..]);
/// }
/// ```
pub fn c2c_1d(in_: &[Complex<f64>], out: &mut [Complex<f64>],
              forward: bool) -> Result<(), FftError> {
    assert!(in_.len() < (1 << (8*mem::size_of::<libc::c_int>() - 1)) - 1);

    if out.len() < in_.len() {
        return Err(FftError::InsufficientSpace(in_.len(), out.len()))
    }

    let sign = if forward { ffi::FFTW_FORWARD } else { ffi::FFTW_BACKWARD };

    unsafe {
        let plan = plan::RawPlan::new(|| {
            ffi::fftw_plan_dft_1d(in_.len() as libc::c_int,
                                  // this is only correct because we use ESTIMATE
                                  in_.as_ptr() as *mut ffi::fftw_complex,
                                  out.as_mut_ptr() as *mut ffi::fftw_complex,
                                  sign,
                                  ffi::FFTW_ESTIMATE)
        });

        match plan {
            None => Err(FftError::FailedPlan),
            Some(mut p) => {
                p.execute();
                Ok(())
            }
        }
    }
}


#[cfg(test)]
mod tests {
    use num::Complex;
    extern crate rand;
    #[test]
    fn c2c_1d_roundtrip() {
        let input = (0..1000)
            .map(|_| Complex::new(rand::random(),
                                  rand::random()))
            .collect::<Vec<_>>();
        let mut output = (0..1000).map(|_| Complex::new(0.0,0.0)).collect::<Vec<_>>();
        let mut output2 = output.clone();

        super::c2c_1d(&*input, &mut *output, true).unwrap();
        super::c2c_1d(&*output, &mut *output2, false).unwrap();

        for (x, raw_y) in input.iter().zip(output2.iter()) {
            let y = raw_y.unscale(1000.0);
            if (*x - y).norm() >= 1e-10 {
                panic!("roundtrip elements too far apart: {}, {}", *x, y)
            }
        }
    }
}