// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Mark Borgerding mark a borgerding net
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
/**
* \defgroup FFT_Module Fast Fourier Transform module
*
* \code
* #include <unsupported/Eigen/FFT>
* \endcode
*
* 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 two implementation backend:
*
* - fftw (http://www.fftw.org) : faster, GPL -- incompatible with Eigen in LGPL form, bigger code size.
* - MKL (http://en.wikipedia.org/wiki/Math_Kernel_Library) : fastest, commercial -- may be incompatible with Eigen in GPL form.
*
* \section FFTDesign 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.
*/
#ifdef EIGEN_FFTW_DEFAULT
// FFTW: faster, GPL -- incompatible with Eigen in LGPL form, bigger code size
# include <fftw3.h>
# include "src/FFT/ei_fftw_impl.h"
namespace Eigen {
//template <typename T> typedef struct internal::fftw_impl default_fft_impl; this does not work
template <typename T> struct default_fft_impl : public internal::fftw_impl<T> {};
}
#elif defined EIGEN_MKL_DEFAULT
// TODO
// intel Math Kernel Library: fastest, commercial -- may be incompatible with Eigen in GPL form
# include "src/FFT/ei_imklfft_impl.h"
namespace Eigen {
template <typename T> struct default_fft_impl : public internal::imklfft_impl {};
}
#else
// internal::kissfft_impl: small, free, reasonably efficient default, derived from kissfft
//
# include "src/FFT/ei_kissfft_impl.h"
namespace Eigen {
template <typename T>
struct default_fft_impl : public internal::kissfft_impl<T> {};
}
#endif
if (nfft<1) { //automatic FFT size determination
if ( realfft && HasFlag(HalfSpectrum) )
nfft = 2*(src.size()-1); //assume even fft size
else
nfft = src.size();
}
dst.derived().resize( nfft );
// check for nfft that does not fit the input data size
Index resize_input= ( realfft && HasFlag(HalfSpectrum) )
? ( (nfft/2+1) - src.size() )
: ( nfft - src.size() );
if ( src.innerStride() != 1 || resize_input ) {
// if the vector is strided, then we need to copy it to a packed temporary
Matrix<src_type,1,Dynamic> tmp;
if ( resize_input ) {
size_t ncopy = (std::min)(src.size(),src.size() + resize_input);
tmp.setZero(src.size() + resize_input);
if ( realfft && HasFlag(HalfSpectrum) ) {
// pad at the Nyquist bin
tmp.head(ncopy) = src.head(ncopy);
tmp(ncopy-1) = real(tmp(ncopy-1)); // enforce real-only Nyquist bin
}else{
size_t nhead,ntail;
nhead = 1+ncopy/2-1; // range [0:pi)
ntail = ncopy/2-1; // range (-pi:0)
tmp.head(nhead) = src.head(nhead);
tmp.tail(ntail) = src.tail(ntail);
if (resize_input<0) { //shrinking -- create the Nyquist bin as the average of the two bins that fold into it
tmp(nhead) = ( src(nfft/2) + src( src.size() - nfft/2 ) )*real_type(.5);
}else{ // expanding -- split the old Nyquist bin into two halves
tmp(nhead) = src(nhead) * real_type(.5);
tmp(tmp.size()-nhead) = tmp(nhead);
}
}
}else{
tmp = src;
}
inv( &dst[0],&tmp[0], nfft);
}else{
inv( &dst[0],&src[0], nfft);
}
}
