EigenFFT

From Eigen
Jump to: navigation, search

Eigen/FFT is a WORK IN PROGRESS

A new effort was started in May 2009 to port/adapt the kissfft library to Eigen. The work is currently in the unsupported section.

Help is welcome. In particular, help is needed in the area of

  • testing and feedback on how to make the interface better
  • creating IMKL backend

Please send an email on the mailing list or to Mark directly in order to coordinate efforts.

The major development milestones are:

  1. DONE: one-dimensional FFTs, forward and inverse for float, double,long double
  2. DONE: real-only optimization
  3. TODO: multi-dimensional FFTs ( more cannibalizing from C version)
  4. DONE: backend for FFTW ( using this requires GPL-compatible application, or the commercial MIT license for FFTW)
  5. TODO: backend for intel Math Kernel Library (for closed source applications)

Example usage

#include <unsupported/Eigen/FFT>

...
Eigen::FFT<float> fft;

std::vector<float> timevec = MakeMyData();
std::vector<std::complex<float> > freqvec;

fft.fwd( freqvec,timevec);
// manipulate freqvec
fft.inv( timevec,freqvec);

// stored "plans" get destroyed with fft destructor


FFT backend options

Allowing different backends makes sure everybody gets what they need based on their individual requirements regarding licensing, finances, code size, and execution speed. The default backend is kissfft since it offers a good balance between easy licensing, size, and speed. The default FFT backend can be chosen at compile time.

kissfft

The FFT module is derived from the kissfft project, but has undergone considerable adaptation for use in C++ and Eigen. The licensing for the derived work has been changed by Mark Borgerding, the kissfft copyright holder, to the Eigen dual LGPL/GPL license.

  • default backend
  • small
  • no external dependencies
  • reasonably fast
  • Compatible with either choice of Eigen license : GPL or LGPL


FFTW

  • relatively large code size
  • very fast
  • compatible with GPL license (NOT LGPL compatible)
  • can be made the default backend by
    • compiling code with preprocessor definition EIGEN_FFTW_DEFAULT
    • linking with FFTW libraries e.g. -lfftw3 -lfftw3f -lfftw3l

Intel Math Kernel Library

TODO This backend has not yet been written.

  • relatively large code size
  • very fast
  • compatible with LGPL license (NOT GPL compatible)

Common features

All three backends have a common set of features

  • arbitrary FFT sizes
  • optimization for real-only FFTs (i.e. real x <=> half spectrum X)
  • multi-dimensional FFTs
  • the concept of a "plan"

Behavior Differences vs other FFT libraries

The Eigen::FFT module has a couple of behaviors that stray away from that of most FFT libraries (e.g. FFTW,IMKL,KISSFFT). These differences are intended to facilitate generic programming and ease migrating algorithms from Matlab/octave. The default behavior of Eigen::FFT favors correctness and generality over speed. The caller can "opt-out" from the behavior and get the speed increase if they want it.

  1. ) Scaling:
    1. Other libraries: do not perform scaling, so there is a constant gain incurred after the forward&inverse transforms , so IFFT(FFT(x)) = Kx; this is done to avoid a vector-by-value multiply. The downside is that algorithms that work correctly in Matlab/octave don't behave the same way once implemented in C++.
    2. How Eigen::FFT differs: invertible scaling is performed so IFFT( FFT(x) ) = x.
    3. Use the Eigen::FFT::Unscaled flag to change the default behavior
  2. ) Real FFT half-spectrum
    1. Other libraries: use only half the frequency spectrum (plus one extra sample for the Nyquist bin) for a real FFT, the other half is the conjugate-symmetric of the first half. This saves a copy and some memory. The downside is the caller needs to have special logic for the number of bins in complex vs real.
    2. How Eigen::FFT differs: The full spectrum is returned from the forward transform. This facilitates generic template programming by obviating separate specializations for real vs complex. On the inverse transform, only half the spectrum is actually used if the expected output type is real.
    3. Use the Eigen::FFT::HalfSpectrum flag to change the default behavior