1// This file is part of Eigen, a lightweight C++ template library
4// Copyright (C) 2009 Mark Borgerding mark a borgerding net
6// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
16#include "../../Eigen/Core"
20 * \defgroup FFT_Module Fast Fourier Transform module
23 * #include <unsupported/Eigen/FFT>
26 * This module provides Fast Fourier transformation, with a configurable backend
29 * The default implementation is based on kissfft. It is a small, free, and
30 * reasonably efficient default.
32 * There are currently two implementation backend:
34 * - fftw (http://www.fftw.org) : faster, GPL -- incompatible with Eigen in LGPL form, bigger code size.
35 * - MKL (http://en.wikipedia.org/wiki/Math_Kernel_Library) : fastest, commercial -- may be incompatible with Eigen in GPL form.
37 * \section FFTDesign Design
39 * The following design decisions were made concerning scaling and
40 * half-spectrum for real FFT.
42 * The intent is to facilitate generic programming and ease migrating code
44 * We think the default behavior of Eigen/FFT should favor correctness and
45 * generality over speed. Of course, the caller should be able to "opt-out" from this
46 * behavior and get the speed increase if they want it.
49 * Other libraries (FFTW,IMKL,KISSFFT) do not perform scaling, so there
50 * is a constant gain incurred after the forward&inverse transforms , so
51 * IFFT(FFT(x)) = Kx; this is done to avoid a vector-by-value multiply.
52 * The downside is that algorithms that worked correctly in Matlab/octave
53 * don't behave the same way once implemented in C++.
55 * How Eigen/FFT differs: invertible scaling is performed so IFFT( FFT(x) ) = x.
57 * 2) Real FFT half-spectrum
58 * Other libraries use only half the frequency spectrum (plus one extra
59 * sample for the Nyquist bin) for a real FFT, the other half is the
60 * conjugate-symmetric of the first half. This saves them a copy and some
61 * memory. The downside is the caller needs to have special logic for the
62 * number of bins in complex vs real.
64 * How Eigen/FFT differs: The full spectrum is returned from the forward
65 * transform. This facilitates generic template programming by obviating
66 * separate specializations for real vs complex. On the inverse
67 * transform, only half the spectrum is actually used if the output type is real.
71#include "../../Eigen/src/Core/util/DisableStupidWarnings.h"
73#ifdef EIGEN_FFTW_DEFAULT
74// FFTW: faster, GPL -- incompatible with Eigen in LGPL form, bigger code size
76# include "src/FFT/ei_fftw_impl.h"
78 //template <typename T> typedef struct internal::fftw_impl default_fft_impl; this does not work
79 template <typename T> struct default_fft_impl : public internal::fftw_impl<T> {};
81#elif defined EIGEN_MKL_DEFAULT
83// intel Math Kernel Library: fastest, commercial -- may be incompatible with Eigen in GPL form
84# include "src/FFT/ei_imklfft_impl.h"
86 template <typename T> struct default_fft_impl : public internal::imklfft_impl {};
89// internal::kissfft_impl: small, free, reasonably efficient default, derived from kissfft
91# include "src/FFT/ei_kissfft_impl.h"
94 struct default_fft_impl : public internal::kissfft_impl<T> {};
102template<typename T_SrcMat,typename T_FftIfc> struct fft_fwd_proxy;
103template<typename T_SrcMat,typename T_FftIfc> struct fft_inv_proxy;
106template<typename T_SrcMat,typename T_FftIfc>
107struct traits< fft_fwd_proxy<T_SrcMat,T_FftIfc> >
109 typedef typename T_SrcMat::PlainObject ReturnType;
111template<typename T_SrcMat,typename T_FftIfc>
112struct traits< fft_inv_proxy<T_SrcMat,T_FftIfc> >
114 typedef typename T_SrcMat::PlainObject ReturnType;
118template<typename T_SrcMat,typename T_FftIfc>
120 : public ReturnByValue<fft_fwd_proxy<T_SrcMat,T_FftIfc> >
122 typedef DenseIndex Index;
124 fft_fwd_proxy(const T_SrcMat& src,T_FftIfc & fft, Index nfft) : m_src(src),m_ifc(fft), m_nfft(nfft) {}
126 template<typename T_DestMat> void evalTo(T_DestMat& dst) const;
128 Index rows() const { return m_src.rows(); }
129 Index cols() const { return m_src.cols(); }
131 const T_SrcMat & m_src;
136template<typename T_SrcMat,typename T_FftIfc>
138 : public ReturnByValue<fft_inv_proxy<T_SrcMat,T_FftIfc> >
140 typedef DenseIndex Index;
142 fft_inv_proxy(const T_SrcMat& src,T_FftIfc & fft, Index nfft) : m_src(src),m_ifc(fft), m_nfft(nfft) {}
144 template<typename T_DestMat> void evalTo(T_DestMat& dst) const;
146 Index rows() const { return m_src.rows(); }
147 Index cols() const { return m_src.cols(); }
149 const T_SrcMat & m_src;
155template <typename T_Scalar,
156 typename T_Impl=default_fft_impl<T_Scalar> >
160 typedef T_Impl impl_type;
161 typedef DenseIndex Index;
162 typedef typename impl_type::Scalar Scalar;
163 typedef typename impl_type::Complex Complex;
166 Default=0, // goof proof
169 // SomeOtherSpeedOptimization=4
173 FFT( const impl_type & impl=impl_type() , Flag flags=Default ) :m_impl(impl),m_flag(flags) { }
176 bool HasFlag(Flag f) const { return (m_flag & (int)f) == f;}
179 void SetFlag(Flag f) { m_flag |= (int)f;}
182 void ClearFlag(Flag f) { m_flag &= (~(int)f);}
185 void fwd( Complex * dst, const Scalar * src, Index nfft)
187 m_impl.fwd(dst,src,static_cast<int>(nfft));
188 if ( HasFlag(HalfSpectrum) == false)
189 ReflectSpectrum(dst,nfft);
193 void fwd( Complex * dst, const Complex * src, Index nfft)
195 m_impl.fwd(dst,src,static_cast<int>(nfft));
200 void fwd2(Complex * dst, const Complex * src, int n0,int n1)
202 m_impl.fwd2(dst,src,n0,n1);
206 template <typename _Input>
208 void fwd( std::vector<Complex> & dst, const std::vector<_Input> & src)
210 if ( NumTraits<_Input>::IsComplex == 0 && HasFlag(HalfSpectrum) )
211 dst.resize( (src.size()>>1)+1); // half the bins + Nyquist bin
213 dst.resize(src.size());
214 fwd(&dst[0],&src[0],src.size());
217 template<typename InputDerived, typename ComplexDerived>
219 void fwd( MatrixBase<ComplexDerived> & dst, const MatrixBase<InputDerived> & src, Index nfft=-1)
221 typedef typename ComplexDerived::Scalar dst_type;
222 typedef typename InputDerived::Scalar src_type;
223 EIGEN_STATIC_ASSERT_VECTOR_ONLY(InputDerived)
224 EIGEN_STATIC_ASSERT_VECTOR_ONLY(ComplexDerived)
225 EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(ComplexDerived,InputDerived) // size at compile-time
226 EIGEN_STATIC_ASSERT((internal::is_same<dst_type, Complex>::value),
227 YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
228 EIGEN_STATIC_ASSERT(int(InputDerived::Flags)&int(ComplexDerived::Flags)&DirectAccessBit,
229 THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_WITH_DIRECT_MEMORY_ACCESS_SUCH_AS_MAP_OR_PLAIN_MATRICES)
234 if ( NumTraits< src_type >::IsComplex == 0 && HasFlag(HalfSpectrum) )
235 dst.derived().resize( (nfft>>1)+1);
237 dst.derived().resize(nfft);
239 if ( src.innerStride() != 1 || src.size() < nfft ) {
240 Matrix<src_type,1,Dynamic> tmp;
241 if (src.size()<nfft) {
243 tmp.block(0,0,src.size(),1 ) = src;
247 fwd( &dst[0],&tmp[0],nfft );
249 fwd( &dst[0],&src[0],nfft );
253 template<typename InputDerived>
255 fft_fwd_proxy< MatrixBase<InputDerived>, FFT<T_Scalar,T_Impl> >
256 fwd( const MatrixBase<InputDerived> & src, Index nfft=-1)
258 return fft_fwd_proxy< MatrixBase<InputDerived> ,FFT<T_Scalar,T_Impl> >( src, *this,nfft );
261 template<typename InputDerived>
263 fft_inv_proxy< MatrixBase<InputDerived>, FFT<T_Scalar,T_Impl> >
264 inv( const MatrixBase<InputDerived> & src, Index nfft=-1)
266 return fft_inv_proxy< MatrixBase<InputDerived> ,FFT<T_Scalar,T_Impl> >( src, *this,nfft );
270 void inv( Complex * dst, const Complex * src, Index nfft)
272 m_impl.inv( dst,src,static_cast<int>(nfft) );
273 if ( HasFlag( Unscaled ) == false)
274 scale(dst,Scalar(1./nfft),nfft); // scale the time series
278 void inv( Scalar * dst, const Complex * src, Index nfft)
280 m_impl.inv( dst,src,static_cast<int>(nfft) );
281 if ( HasFlag( Unscaled ) == false)
282 scale(dst,Scalar(1./nfft),nfft); // scale the time series
285 template<typename OutputDerived, typename ComplexDerived>
287 void inv( MatrixBase<OutputDerived> & dst, const MatrixBase<ComplexDerived> & src, Index nfft=-1)
289 typedef typename ComplexDerived::Scalar src_type;
290 typedef typename ComplexDerived::RealScalar real_type;
291 typedef typename OutputDerived::Scalar dst_type;
292 const bool realfft= (NumTraits<dst_type>::IsComplex == 0);
293 EIGEN_STATIC_ASSERT_VECTOR_ONLY(OutputDerived)
294 EIGEN_STATIC_ASSERT_VECTOR_ONLY(ComplexDerived)
295 EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(ComplexDerived,OutputDerived) // size at compile-time
296 EIGEN_STATIC_ASSERT((internal::is_same<src_type, Complex>::value),
297 YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
298 EIGEN_STATIC_ASSERT(int(OutputDerived::Flags)&int(ComplexDerived::Flags)&DirectAccessBit,
299 THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_WITH_DIRECT_MEMORY_ACCESS_SUCH_AS_MAP_OR_PLAIN_MATRICES)
301 if (nfft<1) { //automatic FFT size determination
302 if ( realfft && HasFlag(HalfSpectrum) )
303 nfft = 2*(src.size()-1); //assume even fft size
307 dst.derived().resize( nfft );
309 // check for nfft that does not fit the input data size
310 Index resize_input= ( realfft && HasFlag(HalfSpectrum) )
311 ? ( (nfft/2+1) - src.size() )
312 : ( nfft - src.size() );
314 if ( src.innerStride() != 1 || resize_input ) {
315 // if the vector is strided, then we need to copy it to a packed temporary
316 Matrix<src_type,1,Dynamic> tmp;
317 if ( resize_input ) {
318 size_t ncopy = (std::min)(src.size(),src.size() + resize_input);
319 tmp.setZero(src.size() + resize_input);
320 if ( realfft && HasFlag(HalfSpectrum) ) {
321 // pad at the Nyquist bin
322 tmp.head(ncopy) = src.head(ncopy);
323 tmp(ncopy-1) = real(tmp(ncopy-1)); // enforce real-only Nyquist bin
326 nhead = 1+ncopy/2-1; // range [0:pi)
327 ntail = ncopy/2-1; // range (-pi:0)
328 tmp.head(nhead) = src.head(nhead);
329 tmp.tail(ntail) = src.tail(ntail);
330 if (resize_input<0) { //shrinking -- create the Nyquist bin as the average of the two bins that fold into it
331 tmp(nhead) = ( src(nfft/2) + src( src.size() - nfft/2 ) )*real_type(.5);
332 }else{ // expanding -- split the old Nyquist bin into two halves
333 tmp(nhead) = src(nhead) * real_type(.5);
334 tmp(tmp.size()-nhead) = tmp(nhead);
340 inv( &dst[0],&tmp[0], nfft);
342 inv( &dst[0],&src[0], nfft);
346 template <typename _Output>
348 void inv( std::vector<_Output> & dst, const std::vector<Complex> & src,Index nfft=-1)
351 nfft = ( NumTraits<_Output>::IsComplex == 0 && HasFlag(HalfSpectrum) ) ? 2*(src.size()-1) : src.size();
353 inv( &dst[0],&src[0],nfft);
358 // TODO: multi-dimensional FFTs
360 void inv2(Complex * dst, const Complex * src, int n0,int n1)
362 m_impl.inv2(dst,src,n0,n1);
363 if ( HasFlag( Unscaled ) == false)
364 scale(dst,1./(n0*n1),n0*n1);
369 impl_type & impl() {return m_impl;}
372 template <typename T_Data>
374 void scale(T_Data * x,Scalar s,Index nx)
377 for (int k=0;k<nx;++k)
380 if ( ((ptrdiff_t)x) & 15 )
381 Matrix<T_Data, Dynamic, 1>::Map(x,nx) *= s;
383 Matrix<T_Data, Dynamic, 1>::MapAligned(x,nx) *= s;
384 //Matrix<T_Data, Dynamic, Dynamic>::Map(x,nx) * s;
389 void ReflectSpectrum(Complex * freq, Index nfft)
391 // create the implicit right-half spectrum (conjugate-mirror of the left-half)
392 Index nhbins=(nfft>>1)+1;
393 for (Index k=nhbins;k < nfft; ++k )
394 freq[k] = conj(freq[nfft-k]);
401template<typename T_SrcMat,typename T_FftIfc>
402template<typename T_DestMat> inline
403void fft_fwd_proxy<T_SrcMat,T_FftIfc>::evalTo(T_DestMat& dst) const
405 m_ifc.fwd( dst, m_src, m_nfft);
408template<typename T_SrcMat,typename T_FftIfc>
409template<typename T_DestMat> inline
410void fft_inv_proxy<T_SrcMat,T_FftIfc>::evalTo(T_DestMat& dst) const
412 m_ifc.inv( dst, m_src, m_nfft);
417#include "../../Eigen/src/Core/util/ReenableStupidWarnings.h"