Fast Fourier Transform module

#include <contrib/Eigen/FFT>

This module provides Fast Fourier transformation, with a configurable backend implementation.

The default implementation is based on kissfft. It is a small, free, and reasonably efficient default.

There are currently four implementation backend:

• kissfft(https://github.com/mborgerding/kissfft) : Simple and not so fast, BSD-3-Clause. It is a mixed-radix Fast Fourier Transform based up on the principle, "Keep It Simple, Stupid." Notice that:kissfft fails to handle "atypically-sized" inputs(i.e., sizes with large factors),a workaround is using fftw or pocketfft.
• fftw (http://www.fftw.org) : faster, GPL – incompatible with Eigen in LGPL form, bigger code size.
• MKL (https://www.intel.com/content/www/us/en/developer/tools/oneapi/onemkl-download.html) : fastest, free – may be incompatible with Eigen in GPL form.
• PocketFFT/DUCC (https://gitlab.mpcdf.mpg.de/mtr/pocketfft, https://gitlab.mpcdf.mpg.de/mtr/ducc) : faster than kissfft, BSD 3-clause. A heavily modified implementation of FFTPack with the following advantages:
  1. strictly C++14 compliant
  2. more accurate twiddle factor computation
  3. very fast plan generation
  4. worst case complexity for transform sizes with large prime factors is N*log(N), because Bluestein's algorithm is used for these cases.
  According to the author, DUCC contains the "evolution" of pocketfft, though the interface is very similar.

Thread safety

An FFT object is not thread-safe: it caches a plan, so concurrent transforms need one object per thread. Distinct objects may be used concurrently.

That is not by itself enough for the FFTW backend, whose planner keeps process-wide state and whose only concurrent entry points are fftw_execute and its new-array variants. Eigen therefore serializes plan creation and destruction on a mutex of its own, which covers every planner call Eigen makes from within one process image.

Two cases fall outside that. FFTW planner calls made directly, not through Eigen, share nothing with Eigen's mutex. And a program that loads several shared libraries which each carry their own copy of Eigen may end up with one mutex per library instead of one per process: Eigen gives the mutex default visibility so that the dynamic linker collapses those copies, which covers ELF libraries linked normally or loaded with RTLD_GLOBAL, but not separate Windows DLLs, and not a library loaded with RTLD_LOCAL unless the compiler emitted a unique symbol for the definition. In either case, call FFTW's own fftw_make_planner_thread_safe() once before any transform: it puts the lock inside the planner, where it applies to every caller. It needs FFTW 3.3.5 or later, built with thread support.

Design

The following design decisions were made concerning scaling and half-spectrum for real FFT.

The intent is to facilitate generic programming and ease migrating code from Matlab/octave. We think the default behavior of Eigen/FFT should favor correctness and generality over speed. Of course, the caller should be able to "opt-out" from this behavior and get the speed increase if they want it.

1) Scaling: Other libraries (FFTW,IMKL,KISSFFT) 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 worked correctly in Matlab/octave don't behave the same way once implemented in C++.

How Eigen/FFT differs: invertible scaling is performed so IFFT( FFT(x) ) = x.

2) Real FFT half-spectrum 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 them a copy and some memory. The downside is the caller needs to have special logic for the number of bins in complex vs real.

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 output type is real.