FFT revision 7faaa9f3f0df9d23790277834d426c3d992ac3ba
1// This file is part of Eigen, a lightweight C++ template library 2// for linear algebra. 3// 4// Copyright (C) 2009 Mark Borgerding mark a borgerding net 5// 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/. 9 10#ifndef EIGEN_FFT_H 11#define EIGEN_FFT_H 12 13#include <complex> 14#include <vector> 15#include <map> 16#include <Eigen/Core> 17 18 19/** 20 * \defgroup FFT_Module Fast Fourier Transform module 21 * 22 * \code 23 * #include <unsupported/Eigen/FFT> 24 * \endcode 25 * 26 * This module provides Fast Fourier transformation, with a configurable backend 27 * implementation. 28 * 29 * The default implementation is based on kissfft. It is a small, free, and 30 * reasonably efficient default. 31 * 32 * There are currently two implementation backend: 33 * 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. 36 * 37 * \section FFTDesign Design 38 * 39 * The following design decisions were made concerning scaling and 40 * half-spectrum for real FFT. 41 * 42 * The intent is to facilitate generic programming and ease migrating code 43 * from Matlab/octave. 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. 47 * 48 * 1) %Scaling: 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++. 54 * 55 * How Eigen/FFT differs: invertible scaling is performed so IFFT( FFT(x) ) = x. 56 * 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. 63 * 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. 68 */ 69 70 71#ifdef EIGEN_FFTW_DEFAULT 72// FFTW: faster, GPL -- incompatible with Eigen in LGPL form, bigger code size 73# include <fftw3.h> 74# include "src/FFT/ei_fftw_impl.h" 75 namespace Eigen { 76 //template <typename T> typedef struct internal::fftw_impl default_fft_impl; this does not work 77 template <typename T> struct default_fft_impl : public internal::fftw_impl<T> {}; 78 } 79#elif defined EIGEN_MKL_DEFAULT 80// TODO 81// intel Math Kernel Library: fastest, commercial -- may be incompatible with Eigen in GPL form 82# include "src/FFT/ei_imklfft_impl.h" 83 namespace Eigen { 84 template <typename T> struct default_fft_impl : public internal::imklfft_impl {}; 85 } 86#else 87// internal::kissfft_impl: small, free, reasonably efficient default, derived from kissfft 88// 89# include "src/FFT/ei_kissfft_impl.h" 90 namespace Eigen { 91 template <typename T> 92 struct default_fft_impl : public internal::kissfft_impl<T> {}; 93 } 94#endif 95 96namespace Eigen { 97 98 99// 100template<typename T_SrcMat,typename T_FftIfc> struct fft_fwd_proxy; 101template<typename T_SrcMat,typename T_FftIfc> struct fft_inv_proxy; 102 103namespace internal { 104template<typename T_SrcMat,typename T_FftIfc> 105struct traits< fft_fwd_proxy<T_SrcMat,T_FftIfc> > 106{ 107 typedef typename T_SrcMat::PlainObject ReturnType; 108}; 109template<typename T_SrcMat,typename T_FftIfc> 110struct traits< fft_inv_proxy<T_SrcMat,T_FftIfc> > 111{ 112 typedef typename T_SrcMat::PlainObject ReturnType; 113}; 114} 115 116template<typename T_SrcMat,typename T_FftIfc> 117struct fft_fwd_proxy 118 : public ReturnByValue<fft_fwd_proxy<T_SrcMat,T_FftIfc> > 119{ 120 typedef DenseIndex Index; 121 122 fft_fwd_proxy(const T_SrcMat& src,T_FftIfc & fft, Index nfft) : m_src(src),m_ifc(fft), m_nfft(nfft) {} 123 124 template<typename T_DestMat> void evalTo(T_DestMat& dst) const; 125 126 Index rows() const { return m_src.rows(); } 127 Index cols() const { return m_src.cols(); } 128protected: 129 const T_SrcMat & m_src; 130 T_FftIfc & m_ifc; 131 Index m_nfft; 132private: 133 fft_fwd_proxy& operator=(const fft_fwd_proxy&); 134}; 135 136template<typename T_SrcMat,typename T_FftIfc> 137struct fft_inv_proxy 138 : public ReturnByValue<fft_inv_proxy<T_SrcMat,T_FftIfc> > 139{ 140 typedef DenseIndex Index; 141 142 fft_inv_proxy(const T_SrcMat& src,T_FftIfc & fft, Index nfft) : m_src(src),m_ifc(fft), m_nfft(nfft) {} 143 144 template<typename T_DestMat> void evalTo(T_DestMat& dst) const; 145 146 Index rows() const { return m_src.rows(); } 147 Index cols() const { return m_src.cols(); } 148protected: 149 const T_SrcMat & m_src; 150 T_FftIfc & m_ifc; 151 Index m_nfft; 152private: 153 fft_inv_proxy& operator=(const fft_inv_proxy&); 154}; 155 156 157template <typename T_Scalar, 158 typename T_Impl=default_fft_impl<T_Scalar> > 159class FFT 160{ 161 public: 162 typedef T_Impl impl_type; 163 typedef DenseIndex Index; 164 typedef typename impl_type::Scalar Scalar; 165 typedef typename impl_type::Complex Complex; 166 167 enum Flag { 168 Default=0, // goof proof 169 Unscaled=1, 170 HalfSpectrum=2, 171 // SomeOtherSpeedOptimization=4 172 Speedy=32767 173 }; 174 175 FFT( const impl_type & impl=impl_type() , Flag flags=Default ) :m_impl(impl),m_flag(flags) { } 176 177 inline 178 bool HasFlag(Flag f) const { return (m_flag & (int)f) == f;} 179 180 inline 181 void SetFlag(Flag f) { m_flag |= (int)f;} 182 183 inline 184 void ClearFlag(Flag f) { m_flag &= (~(int)f);} 185 186 inline 187 void fwd( Complex * dst, const Scalar * src, Index nfft) 188 { 189 m_impl.fwd(dst,src,static_cast<int>(nfft)); 190 if ( HasFlag(HalfSpectrum) == false) 191 ReflectSpectrum(dst,nfft); 192 } 193 194 inline 195 void fwd( Complex * dst, const Complex * src, Index nfft) 196 { 197 m_impl.fwd(dst,src,static_cast<int>(nfft)); 198 } 199 200 /* 201 inline 202 void fwd2(Complex * dst, const Complex * src, int n0,int n1) 203 { 204 m_impl.fwd2(dst,src,n0,n1); 205 } 206 */ 207 208 template <typename _Input> 209 inline 210 void fwd( std::vector<Complex> & dst, const std::vector<_Input> & src) 211 { 212 if ( NumTraits<_Input>::IsComplex == 0 && HasFlag(HalfSpectrum) ) 213 dst.resize( (src.size()>>1)+1); // half the bins + Nyquist bin 214 else 215 dst.resize(src.size()); 216 fwd(&dst[0],&src[0],src.size()); 217 } 218 219 template<typename InputDerived, typename ComplexDerived> 220 inline 221 void fwd( MatrixBase<ComplexDerived> & dst, const MatrixBase<InputDerived> & src, Index nfft=-1) 222 { 223 typedef typename ComplexDerived::Scalar dst_type; 224 typedef typename InputDerived::Scalar src_type; 225 EIGEN_STATIC_ASSERT_VECTOR_ONLY(InputDerived) 226 EIGEN_STATIC_ASSERT_VECTOR_ONLY(ComplexDerived) 227 EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(ComplexDerived,InputDerived) // size at compile-time 228 EIGEN_STATIC_ASSERT((internal::is_same<dst_type, Complex>::value), 229 YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) 230 EIGEN_STATIC_ASSERT(int(InputDerived::Flags)&int(ComplexDerived::Flags)&DirectAccessBit, 231 THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_WITH_DIRECT_MEMORY_ACCESS_SUCH_AS_MAP_OR_PLAIN_MATRICES) 232 233 if (nfft<1) 234 nfft = src.size(); 235 236 if ( NumTraits< src_type >::IsComplex == 0 && HasFlag(HalfSpectrum) ) 237 dst.derived().resize( (nfft>>1)+1); 238 else 239 dst.derived().resize(nfft); 240 241 if ( src.innerStride() != 1 || src.size() < nfft ) { 242 Matrix<src_type,1,Dynamic> tmp; 243 if (src.size()<nfft) { 244 tmp.setZero(nfft); 245 tmp.block(0,0,src.size(),1 ) = src; 246 }else{ 247 tmp = src; 248 } 249 fwd( &dst[0],&tmp[0],nfft ); 250 }else{ 251 fwd( &dst[0],&src[0],nfft ); 252 } 253 } 254 255 template<typename InputDerived> 256 inline 257 fft_fwd_proxy< MatrixBase<InputDerived>, FFT<T_Scalar,T_Impl> > 258 fwd( const MatrixBase<InputDerived> & src, Index nfft=-1) 259 { 260 return fft_fwd_proxy< MatrixBase<InputDerived> ,FFT<T_Scalar,T_Impl> >( src, *this,nfft ); 261 } 262 263 template<typename InputDerived> 264 inline 265 fft_inv_proxy< MatrixBase<InputDerived>, FFT<T_Scalar,T_Impl> > 266 inv( const MatrixBase<InputDerived> & src, Index nfft=-1) 267 { 268 return fft_inv_proxy< MatrixBase<InputDerived> ,FFT<T_Scalar,T_Impl> >( src, *this,nfft ); 269 } 270 271 inline 272 void inv( Complex * dst, const Complex * src, Index nfft) 273 { 274 m_impl.inv( dst,src,static_cast<int>(nfft) ); 275 if ( HasFlag( Unscaled ) == false) 276 scale(dst,Scalar(1./nfft),nfft); // scale the time series 277 } 278 279 inline 280 void inv( Scalar * dst, const Complex * src, Index nfft) 281 { 282 m_impl.inv( dst,src,static_cast<int>(nfft) ); 283 if ( HasFlag( Unscaled ) == false) 284 scale(dst,Scalar(1./nfft),nfft); // scale the time series 285 } 286 287 template<typename OutputDerived, typename ComplexDerived> 288 inline 289 void inv( MatrixBase<OutputDerived> & dst, const MatrixBase<ComplexDerived> & src, Index nfft=-1) 290 { 291 typedef typename ComplexDerived::Scalar src_type; 292 typedef typename OutputDerived::Scalar dst_type; 293 const bool realfft= (NumTraits<dst_type>::IsComplex == 0); 294 EIGEN_STATIC_ASSERT_VECTOR_ONLY(OutputDerived) 295 EIGEN_STATIC_ASSERT_VECTOR_ONLY(ComplexDerived) 296 EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(ComplexDerived,OutputDerived) // size at compile-time 297 EIGEN_STATIC_ASSERT((internal::is_same<src_type, Complex>::value), 298 YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY) 299 EIGEN_STATIC_ASSERT(int(OutputDerived::Flags)&int(ComplexDerived::Flags)&DirectAccessBit, 300 THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_WITH_DIRECT_MEMORY_ACCESS_SUCH_AS_MAP_OR_PLAIN_MATRICES) 301 302 if (nfft<1) { //automatic FFT size determination 303 if ( realfft && HasFlag(HalfSpectrum) ) 304 nfft = 2*(src.size()-1); //assume even fft size 305 else 306 nfft = src.size(); 307 } 308 dst.derived().resize( nfft ); 309 310 // check for nfft that does not fit the input data size 311 Index resize_input= ( realfft && HasFlag(HalfSpectrum) ) 312 ? ( (nfft/2+1) - src.size() ) 313 : ( nfft - src.size() ); 314 315 if ( src.innerStride() != 1 || resize_input ) { 316 // if the vector is strided, then we need to copy it to a packed temporary 317 Matrix<src_type,1,Dynamic> tmp; 318 if ( resize_input ) { 319 size_t ncopy = (std::min)(src.size(),src.size() + resize_input); 320 tmp.setZero(src.size() + resize_input); 321 if ( realfft && HasFlag(HalfSpectrum) ) { 322 // pad at the Nyquist bin 323 tmp.head(ncopy) = src.head(ncopy); 324 tmp(ncopy-1) = real(tmp(ncopy-1)); // enforce real-only Nyquist bin 325 }else{ 326 size_t nhead,ntail; 327 nhead = 1+ncopy/2-1; // range [0:pi) 328 ntail = ncopy/2-1; // range (-pi:0) 329 tmp.head(nhead) = src.head(nhead); 330 tmp.tail(ntail) = src.tail(ntail); 331 if (resize_input<0) { //shrinking -- create the Nyquist bin as the average of the two bins that fold into it 332 tmp(nhead) = ( src(nfft/2) + src( src.size() - nfft/2 ) )*src_type(.5); 333 }else{ // expanding -- split the old Nyquist bin into two halves 334 tmp(nhead) = src(nhead) * src_type(.5); 335 tmp(tmp.size()-nhead) = tmp(nhead); 336 } 337 } 338 }else{ 339 tmp = src; 340 } 341 inv( &dst[0],&tmp[0], nfft); 342 }else{ 343 inv( &dst[0],&src[0], nfft); 344 } 345 } 346 347 template <typename _Output> 348 inline 349 void inv( std::vector<_Output> & dst, const std::vector<Complex> & src,Index nfft=-1) 350 { 351 if (nfft<1) 352 nfft = ( NumTraits<_Output>::IsComplex == 0 && HasFlag(HalfSpectrum) ) ? 2*(src.size()-1) : src.size(); 353 dst.resize( nfft ); 354 inv( &dst[0],&src[0],nfft); 355 } 356 357 358 /* 359 // TODO: multi-dimensional FFTs 360 inline 361 void inv2(Complex * dst, const Complex * src, int n0,int n1) 362 { 363 m_impl.inv2(dst,src,n0,n1); 364 if ( HasFlag( Unscaled ) == false) 365 scale(dst,1./(n0*n1),n0*n1); 366 } 367 */ 368 369 inline 370 impl_type & impl() {return m_impl;} 371 private: 372 373 template <typename T_Data> 374 inline 375 void scale(T_Data * x,Scalar s,Index nx) 376 { 377#if 1 378 for (int k=0;k<nx;++k) 379 *x++ *= s; 380#else 381 if ( ((ptrdiff_t)x) & 15 ) 382 Matrix<T_Data, Dynamic, 1>::Map(x,nx) *= s; 383 else 384 Matrix<T_Data, Dynamic, 1>::MapAligned(x,nx) *= s; 385 //Matrix<T_Data, Dynamic, Dynamic>::Map(x,nx) * s; 386#endif 387 } 388 389 inline 390 void ReflectSpectrum(Complex * freq, Index nfft) 391 { 392 // create the implicit right-half spectrum (conjugate-mirror of the left-half) 393 Index nhbins=(nfft>>1)+1; 394 for (Index k=nhbins;k < nfft; ++k ) 395 freq[k] = conj(freq[nfft-k]); 396 } 397 398 impl_type m_impl; 399 int m_flag; 400}; 401 402template<typename T_SrcMat,typename T_FftIfc> 403template<typename T_DestMat> inline 404void fft_fwd_proxy<T_SrcMat,T_FftIfc>::evalTo(T_DestMat& dst) const 405{ 406 m_ifc.fwd( dst, m_src, m_nfft); 407} 408 409template<typename T_SrcMat,typename T_FftIfc> 410template<typename T_DestMat> inline 411void fft_inv_proxy<T_SrcMat,T_FftIfc>::evalTo(T_DestMat& dst) const 412{ 413 m_ifc.inv( dst, m_src, m_nfft); 414} 415 416} 417#endif 418/* vim: set filetype=cpp et sw=2 ts=2 ai: */ 419