Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 19 | /** |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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 | |
| 96 | namespace Eigen { |
| 97 | |
| 98 | |
| 99 | // |
| 100 | template<typename T_SrcMat,typename T_FftIfc> struct fft_fwd_proxy; |
| 101 | template<typename T_SrcMat,typename T_FftIfc> struct fft_inv_proxy; |
| 102 | |
| 103 | namespace internal { |
| 104 | template<typename T_SrcMat,typename T_FftIfc> |
| 105 | struct traits< fft_fwd_proxy<T_SrcMat,T_FftIfc> > |
| 106 | { |
| 107 | typedef typename T_SrcMat::PlainObject ReturnType; |
| 108 | }; |
| 109 | template<typename T_SrcMat,typename T_FftIfc> |
| 110 | struct traits< fft_inv_proxy<T_SrcMat,T_FftIfc> > |
| 111 | { |
| 112 | typedef typename T_SrcMat::PlainObject ReturnType; |
| 113 | }; |
| 114 | } |
| 115 | |
| 116 | template<typename T_SrcMat,typename T_FftIfc> |
| 117 | struct 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(); } |
| 128 | protected: |
| 129 | const T_SrcMat & m_src; |
| 130 | T_FftIfc & m_ifc; |
| 131 | Index m_nfft; |
| 132 | private: |
| 133 | fft_fwd_proxy& operator=(const fft_fwd_proxy&); |
| 134 | }; |
| 135 | |
| 136 | template<typename T_SrcMat,typename T_FftIfc> |
| 137 | struct 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(); } |
| 148 | protected: |
| 149 | const T_SrcMat & m_src; |
| 150 | T_FftIfc & m_ifc; |
| 151 | Index m_nfft; |
| 152 | private: |
| 153 | fft_inv_proxy& operator=(const fft_inv_proxy&); |
| 154 | }; |
| 155 | |
| 156 | |
| 157 | template <typename T_Scalar, |
| 158 | typename T_Impl=default_fft_impl<T_Scalar> > |
| 159 | class 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 | |
| 402 | template<typename T_SrcMat,typename T_FftIfc> |
| 403 | template<typename T_DestMat> inline |
| 404 | void fft_fwd_proxy<T_SrcMat,T_FftIfc>::evalTo(T_DestMat& dst) const |
| 405 | { |
| 406 | m_ifc.fwd( dst, m_src, m_nfft); |
| 407 | } |
| 408 | |
| 409 | template<typename T_SrcMat,typename T_FftIfc> |
| 410 | template<typename T_DestMat> inline |
| 411 | void 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: */ |