Warning, /include/eigen3/unsupported/Eigen/FFT is written in an unsupported language. File is not indexed.
0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2009 Mark Borgerding mark a borgerding net
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_FFT_H
0011 #define EIGEN_FFT_H
0012
0013 #include <complex>
0014 #include <vector>
0015 #include <map>
0016 #include "../../Eigen/Core"
0017
0018
0019 /**
0020 * \defgroup FFT_Module Fast Fourier Transform module
0021 *
0022 * \code
0023 * #include <unsupported/Eigen/FFT>
0024 * \endcode
0025 *
0026 * This module provides Fast Fourier transformation, with a configurable backend
0027 * implementation.
0028 *
0029 * The default implementation is based on kissfft. It is a small, free, and
0030 * reasonably efficient default.
0031 *
0032 * There are currently two implementation backend:
0033 *
0034 * - fftw (http://www.fftw.org) : faster, GPL -- incompatible with Eigen in LGPL form, bigger code size.
0035 * - MKL (http://en.wikipedia.org/wiki/Math_Kernel_Library) : fastest, commercial -- may be incompatible with Eigen in GPL form.
0036 *
0037 * \section FFTDesign Design
0038 *
0039 * The following design decisions were made concerning scaling and
0040 * half-spectrum for real FFT.
0041 *
0042 * The intent is to facilitate generic programming and ease migrating code
0043 * from Matlab/octave.
0044 * We think the default behavior of Eigen/FFT should favor correctness and
0045 * generality over speed. Of course, the caller should be able to "opt-out" from this
0046 * behavior and get the speed increase if they want it.
0047 *
0048 * 1) %Scaling:
0049 * Other libraries (FFTW,IMKL,KISSFFT) do not perform scaling, so there
0050 * is a constant gain incurred after the forward&inverse transforms , so
0051 * IFFT(FFT(x)) = Kx; this is done to avoid a vector-by-value multiply.
0052 * The downside is that algorithms that worked correctly in Matlab/octave
0053 * don't behave the same way once implemented in C++.
0054 *
0055 * How Eigen/FFT differs: invertible scaling is performed so IFFT( FFT(x) ) = x.
0056 *
0057 * 2) Real FFT half-spectrum
0058 * Other libraries use only half the frequency spectrum (plus one extra
0059 * sample for the Nyquist bin) for a real FFT, the other half is the
0060 * conjugate-symmetric of the first half. This saves them a copy and some
0061 * memory. The downside is the caller needs to have special logic for the
0062 * number of bins in complex vs real.
0063 *
0064 * How Eigen/FFT differs: The full spectrum is returned from the forward
0065 * transform. This facilitates generic template programming by obviating
0066 * separate specializations for real vs complex. On the inverse
0067 * transform, only half the spectrum is actually used if the output type is real.
0068 */
0069
0070
0071 #include "../../Eigen/src/Core/util/DisableStupidWarnings.h"
0072
0073 #ifdef EIGEN_FFTW_DEFAULT
0074 // FFTW: faster, GPL -- incompatible with Eigen in LGPL form, bigger code size
0075 # include <fftw3.h>
0076 # include "src/FFT/ei_fftw_impl.h"
0077 namespace Eigen {
0078 //template <typename T> typedef struct internal::fftw_impl default_fft_impl; this does not work
0079 template <typename T> struct default_fft_impl : public internal::fftw_impl<T> {};
0080 }
0081 #elif defined EIGEN_MKL_DEFAULT
0082 // TODO
0083 // intel Math Kernel Library: fastest, commercial -- may be incompatible with Eigen in GPL form
0084 # include "src/FFT/ei_imklfft_impl.h"
0085 namespace Eigen {
0086 template <typename T> struct default_fft_impl : public internal::imklfft_impl {};
0087 }
0088 #else
0089 // internal::kissfft_impl: small, free, reasonably efficient default, derived from kissfft
0090 //
0091 # include "src/FFT/ei_kissfft_impl.h"
0092 namespace Eigen {
0093 template <typename T>
0094 struct default_fft_impl : public internal::kissfft_impl<T> {};
0095 }
0096 #endif
0097
0098 namespace Eigen {
0099
0100
0101 //
0102 template<typename T_SrcMat,typename T_FftIfc> struct fft_fwd_proxy;
0103 template<typename T_SrcMat,typename T_FftIfc> struct fft_inv_proxy;
0104
0105 namespace internal {
0106 template<typename T_SrcMat,typename T_FftIfc>
0107 struct traits< fft_fwd_proxy<T_SrcMat,T_FftIfc> >
0108 {
0109 typedef typename T_SrcMat::PlainObject ReturnType;
0110 };
0111 template<typename T_SrcMat,typename T_FftIfc>
0112 struct traits< fft_inv_proxy<T_SrcMat,T_FftIfc> >
0113 {
0114 typedef typename T_SrcMat::PlainObject ReturnType;
0115 };
0116 }
0117
0118 template<typename T_SrcMat,typename T_FftIfc>
0119 struct fft_fwd_proxy
0120 : public ReturnByValue<fft_fwd_proxy<T_SrcMat,T_FftIfc> >
0121 {
0122 typedef DenseIndex Index;
0123
0124 fft_fwd_proxy(const T_SrcMat& src,T_FftIfc & fft, Index nfft) : m_src(src),m_ifc(fft), m_nfft(nfft) {}
0125
0126 template<typename T_DestMat> void evalTo(T_DestMat& dst) const;
0127
0128 Index rows() const { return m_src.rows(); }
0129 Index cols() const { return m_src.cols(); }
0130 protected:
0131 const T_SrcMat & m_src;
0132 T_FftIfc & m_ifc;
0133 Index m_nfft;
0134 };
0135
0136 template<typename T_SrcMat,typename T_FftIfc>
0137 struct fft_inv_proxy
0138 : public ReturnByValue<fft_inv_proxy<T_SrcMat,T_FftIfc> >
0139 {
0140 typedef DenseIndex Index;
0141
0142 fft_inv_proxy(const T_SrcMat& src,T_FftIfc & fft, Index nfft) : m_src(src),m_ifc(fft), m_nfft(nfft) {}
0143
0144 template<typename T_DestMat> void evalTo(T_DestMat& dst) const;
0145
0146 Index rows() const { return m_src.rows(); }
0147 Index cols() const { return m_src.cols(); }
0148 protected:
0149 const T_SrcMat & m_src;
0150 T_FftIfc & m_ifc;
0151 Index m_nfft;
0152 };
0153
0154
0155 template <typename T_Scalar,
0156 typename T_Impl=default_fft_impl<T_Scalar> >
0157 class FFT
0158 {
0159 public:
0160 typedef T_Impl impl_type;
0161 typedef DenseIndex Index;
0162 typedef typename impl_type::Scalar Scalar;
0163 typedef typename impl_type::Complex Complex;
0164
0165 enum Flag {
0166 Default=0, // goof proof
0167 Unscaled=1,
0168 HalfSpectrum=2,
0169 // SomeOtherSpeedOptimization=4
0170 Speedy=32767
0171 };
0172
0173 FFT( const impl_type & impl=impl_type() , Flag flags=Default ) :m_impl(impl),m_flag(flags) { }
0174
0175 inline
0176 bool HasFlag(Flag f) const { return (m_flag & (int)f) == f;}
0177
0178 inline
0179 void SetFlag(Flag f) { m_flag |= (int)f;}
0180
0181 inline
0182 void ClearFlag(Flag f) { m_flag &= (~(int)f);}
0183
0184 inline
0185 void fwd( Complex * dst, const Scalar * src, Index nfft)
0186 {
0187 m_impl.fwd(dst,src,static_cast<int>(nfft));
0188 if ( HasFlag(HalfSpectrum) == false)
0189 ReflectSpectrum(dst,nfft);
0190 }
0191
0192 inline
0193 void fwd( Complex * dst, const Complex * src, Index nfft)
0194 {
0195 m_impl.fwd(dst,src,static_cast<int>(nfft));
0196 }
0197
0198 /*
0199 inline
0200 void fwd2(Complex * dst, const Complex * src, int n0,int n1)
0201 {
0202 m_impl.fwd2(dst,src,n0,n1);
0203 }
0204 */
0205
0206 template <typename _Input>
0207 inline
0208 void fwd( std::vector<Complex> & dst, const std::vector<_Input> & src)
0209 {
0210 if ( NumTraits<_Input>::IsComplex == 0 && HasFlag(HalfSpectrum) )
0211 dst.resize( (src.size()>>1)+1); // half the bins + Nyquist bin
0212 else
0213 dst.resize(src.size());
0214 fwd(&dst[0],&src[0],src.size());
0215 }
0216
0217 template<typename InputDerived, typename ComplexDerived>
0218 inline
0219 void fwd( MatrixBase<ComplexDerived> & dst, const MatrixBase<InputDerived> & src, Index nfft=-1)
0220 {
0221 typedef typename ComplexDerived::Scalar dst_type;
0222 typedef typename InputDerived::Scalar src_type;
0223 EIGEN_STATIC_ASSERT_VECTOR_ONLY(InputDerived)
0224 EIGEN_STATIC_ASSERT_VECTOR_ONLY(ComplexDerived)
0225 EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(ComplexDerived,InputDerived) // size at compile-time
0226 EIGEN_STATIC_ASSERT((internal::is_same<dst_type, Complex>::value),
0227 YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
0228 EIGEN_STATIC_ASSERT(int(InputDerived::Flags)&int(ComplexDerived::Flags)&DirectAccessBit,
0229 THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_WITH_DIRECT_MEMORY_ACCESS_SUCH_AS_MAP_OR_PLAIN_MATRICES)
0230
0231 if (nfft<1)
0232 nfft = src.size();
0233
0234 if ( NumTraits< src_type >::IsComplex == 0 && HasFlag(HalfSpectrum) )
0235 dst.derived().resize( (nfft>>1)+1);
0236 else
0237 dst.derived().resize(nfft);
0238
0239 if ( src.innerStride() != 1 || src.size() < nfft ) {
0240 Matrix<src_type,1,Dynamic> tmp;
0241 if (src.size()<nfft) {
0242 tmp.setZero(nfft);
0243 tmp.block(0,0,src.size(),1 ) = src;
0244 }else{
0245 tmp = src;
0246 }
0247 fwd( &dst[0],&tmp[0],nfft );
0248 }else{
0249 fwd( &dst[0],&src[0],nfft );
0250 }
0251 }
0252
0253 template<typename InputDerived>
0254 inline
0255 fft_fwd_proxy< MatrixBase<InputDerived>, FFT<T_Scalar,T_Impl> >
0256 fwd( const MatrixBase<InputDerived> & src, Index nfft=-1)
0257 {
0258 return fft_fwd_proxy< MatrixBase<InputDerived> ,FFT<T_Scalar,T_Impl> >( src, *this,nfft );
0259 }
0260
0261 template<typename InputDerived>
0262 inline
0263 fft_inv_proxy< MatrixBase<InputDerived>, FFT<T_Scalar,T_Impl> >
0264 inv( const MatrixBase<InputDerived> & src, Index nfft=-1)
0265 {
0266 return fft_inv_proxy< MatrixBase<InputDerived> ,FFT<T_Scalar,T_Impl> >( src, *this,nfft );
0267 }
0268
0269 inline
0270 void inv( Complex * dst, const Complex * src, Index nfft)
0271 {
0272 m_impl.inv( dst,src,static_cast<int>(nfft) );
0273 if ( HasFlag( Unscaled ) == false)
0274 scale(dst,Scalar(1./nfft),nfft); // scale the time series
0275 }
0276
0277 inline
0278 void inv( Scalar * dst, const Complex * src, Index nfft)
0279 {
0280 m_impl.inv( dst,src,static_cast<int>(nfft) );
0281 if ( HasFlag( Unscaled ) == false)
0282 scale(dst,Scalar(1./nfft),nfft); // scale the time series
0283 }
0284
0285 template<typename OutputDerived, typename ComplexDerived>
0286 inline
0287 void inv( MatrixBase<OutputDerived> & dst, const MatrixBase<ComplexDerived> & src, Index nfft=-1)
0288 {
0289 typedef typename ComplexDerived::Scalar src_type;
0290 typedef typename ComplexDerived::RealScalar real_type;
0291 typedef typename OutputDerived::Scalar dst_type;
0292 const bool realfft= (NumTraits<dst_type>::IsComplex == 0);
0293 EIGEN_STATIC_ASSERT_VECTOR_ONLY(OutputDerived)
0294 EIGEN_STATIC_ASSERT_VECTOR_ONLY(ComplexDerived)
0295 EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(ComplexDerived,OutputDerived) // size at compile-time
0296 EIGEN_STATIC_ASSERT((internal::is_same<src_type, Complex>::value),
0297 YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
0298 EIGEN_STATIC_ASSERT(int(OutputDerived::Flags)&int(ComplexDerived::Flags)&DirectAccessBit,
0299 THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_WITH_DIRECT_MEMORY_ACCESS_SUCH_AS_MAP_OR_PLAIN_MATRICES)
0300
0301 if (nfft<1) { //automatic FFT size determination
0302 if ( realfft && HasFlag(HalfSpectrum) )
0303 nfft = 2*(src.size()-1); //assume even fft size
0304 else
0305 nfft = src.size();
0306 }
0307 dst.derived().resize( nfft );
0308
0309 // check for nfft that does not fit the input data size
0310 Index resize_input= ( realfft && HasFlag(HalfSpectrum) )
0311 ? ( (nfft/2+1) - src.size() )
0312 : ( nfft - src.size() );
0313
0314 if ( src.innerStride() != 1 || resize_input ) {
0315 // if the vector is strided, then we need to copy it to a packed temporary
0316 Matrix<src_type,1,Dynamic> tmp;
0317 if ( resize_input ) {
0318 size_t ncopy = (std::min)(src.size(),src.size() + resize_input);
0319 tmp.setZero(src.size() + resize_input);
0320 if ( realfft && HasFlag(HalfSpectrum) ) {
0321 // pad at the Nyquist bin
0322 tmp.head(ncopy) = src.head(ncopy);
0323 tmp(ncopy-1) = real(tmp(ncopy-1)); // enforce real-only Nyquist bin
0324 }else{
0325 size_t nhead,ntail;
0326 nhead = 1+ncopy/2-1; // range [0:pi)
0327 ntail = ncopy/2-1; // range (-pi:0)
0328 tmp.head(nhead) = src.head(nhead);
0329 tmp.tail(ntail) = src.tail(ntail);
0330 if (resize_input<0) { //shrinking -- create the Nyquist bin as the average of the two bins that fold into it
0331 tmp(nhead) = ( src(nfft/2) + src( src.size() - nfft/2 ) )*real_type(.5);
0332 }else{ // expanding -- split the old Nyquist bin into two halves
0333 tmp(nhead) = src(nhead) * real_type(.5);
0334 tmp(tmp.size()-nhead) = tmp(nhead);
0335 }
0336 }
0337 }else{
0338 tmp = src;
0339 }
0340 inv( &dst[0],&tmp[0], nfft);
0341 }else{
0342 inv( &dst[0],&src[0], nfft);
0343 }
0344 }
0345
0346 template <typename _Output>
0347 inline
0348 void inv( std::vector<_Output> & dst, const std::vector<Complex> & src,Index nfft=-1)
0349 {
0350 if (nfft<1)
0351 nfft = ( NumTraits<_Output>::IsComplex == 0 && HasFlag(HalfSpectrum) ) ? 2*(src.size()-1) : src.size();
0352 dst.resize( nfft );
0353 inv( &dst[0],&src[0],nfft);
0354 }
0355
0356
0357 /*
0358 // TODO: multi-dimensional FFTs
0359 inline
0360 void inv2(Complex * dst, const Complex * src, int n0,int n1)
0361 {
0362 m_impl.inv2(dst,src,n0,n1);
0363 if ( HasFlag( Unscaled ) == false)
0364 scale(dst,1./(n0*n1),n0*n1);
0365 }
0366 */
0367
0368 inline
0369 impl_type & impl() {return m_impl;}
0370 private:
0371
0372 template <typename T_Data>
0373 inline
0374 void scale(T_Data * x,Scalar s,Index nx)
0375 {
0376 #if 1
0377 for (int k=0;k<nx;++k)
0378 *x++ *= s;
0379 #else
0380 if ( ((ptrdiff_t)x) & 15 )
0381 Matrix<T_Data, Dynamic, 1>::Map(x,nx) *= s;
0382 else
0383 Matrix<T_Data, Dynamic, 1>::MapAligned(x,nx) *= s;
0384 //Matrix<T_Data, Dynamic, Dynamic>::Map(x,nx) * s;
0385 #endif
0386 }
0387
0388 inline
0389 void ReflectSpectrum(Complex * freq, Index nfft)
0390 {
0391 // create the implicit right-half spectrum (conjugate-mirror of the left-half)
0392 Index nhbins=(nfft>>1)+1;
0393 for (Index k=nhbins;k < nfft; ++k )
0394 freq[k] = conj(freq[nfft-k]);
0395 }
0396
0397 impl_type m_impl;
0398 int m_flag;
0399 };
0400
0401 template<typename T_SrcMat,typename T_FftIfc>
0402 template<typename T_DestMat> inline
0403 void fft_fwd_proxy<T_SrcMat,T_FftIfc>::evalTo(T_DestMat& dst) const
0404 {
0405 m_ifc.fwd( dst, m_src, m_nfft);
0406 }
0407
0408 template<typename T_SrcMat,typename T_FftIfc>
0409 template<typename T_DestMat> inline
0410 void fft_inv_proxy<T_SrcMat,T_FftIfc>::evalTo(T_DestMat& dst) const
0411 {
0412 m_ifc.inv( dst, m_src, m_nfft);
0413 }
0414
0415 }
0416
0417 #include "../../Eigen/src/Core/util/ReenableStupidWarnings.h"
0418
0419 #endif
