xref: /aosp_15_r20/external/eigen/unsupported/Eigen/FFT (revision bf2c37156dfe67e5dfebd6d394bad8b2ab5804d4) 
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#include "../../Eigen/src/Core/util/DisableStupidWarnings.h"
72
73#ifdef EIGEN_FFTW_DEFAULT
74// FFTW: faster, GPL -- incompatible with Eigen in LGPL form, bigger code size
75#  include <fftw3.h>
76#  include "src/FFT/ei_fftw_impl.h"
77   namespace Eigen {
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> {};
80   }
81#elif defined EIGEN_MKL_DEFAULT
82// TODO
83// intel Math Kernel Library: fastest, commercial -- may be incompatible with Eigen in GPL form
84#  include "src/FFT/ei_imklfft_impl.h"
85   namespace Eigen {
86     template <typename T> struct default_fft_impl : public internal::imklfft_impl {};
87   }
88#else
89// internal::kissfft_impl:  small, free, reasonably efficient default, derived from kissfft
90//
91# include "src/FFT/ei_kissfft_impl.h"
92  namespace Eigen {
93     template <typename T>
94       struct default_fft_impl : public internal::kissfft_impl<T> {};
95  }
96#endif
97
98namespace Eigen {
99
100
101//
102template<typename T_SrcMat,typename T_FftIfc> struct fft_fwd_proxy;
103template<typename T_SrcMat,typename T_FftIfc> struct fft_inv_proxy;
104
105namespace internal {
106template<typename T_SrcMat,typename T_FftIfc>
107struct traits< fft_fwd_proxy<T_SrcMat,T_FftIfc> >
108{
109  typedef typename T_SrcMat::PlainObject ReturnType;
110};
111template<typename T_SrcMat,typename T_FftIfc>
112struct traits< fft_inv_proxy<T_SrcMat,T_FftIfc> >
113{
114  typedef typename T_SrcMat::PlainObject ReturnType;
115};
116}
117
118template<typename T_SrcMat,typename T_FftIfc>
119struct fft_fwd_proxy
120 : public ReturnByValue<fft_fwd_proxy<T_SrcMat,T_FftIfc> >
121{
122  typedef DenseIndex Index;
123
124  fft_fwd_proxy(const T_SrcMat& src,T_FftIfc & fft, Index nfft) : m_src(src),m_ifc(fft), m_nfft(nfft) {}
125
126  template<typename T_DestMat> void evalTo(T_DestMat& dst) const;
127
128  Index rows() const { return m_src.rows(); }
129  Index cols() const { return m_src.cols(); }
130protected:
131  const T_SrcMat & m_src;
132  T_FftIfc & m_ifc;
133  Index m_nfft;
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;
152};
153
154
155template <typename T_Scalar,
156         typename T_Impl=default_fft_impl<T_Scalar> >
157class FFT
158{
159  public:
160    typedef T_Impl impl_type;
161    typedef DenseIndex Index;
162    typedef typename impl_type::Scalar Scalar;
163    typedef typename impl_type::Complex Complex;
164
165    enum Flag {
166      Default=0, // goof proof
167      Unscaled=1,
168      HalfSpectrum=2,
169      // SomeOtherSpeedOptimization=4
170      Speedy=32767
171    };
172
173    FFT( const impl_type & impl=impl_type() , Flag flags=Default ) :m_impl(impl),m_flag(flags) { }
174
175    inline
176    bool HasFlag(Flag f) const { return (m_flag & (int)f) == f;}
177
178    inline
179    void SetFlag(Flag f) { m_flag |= (int)f;}
180
181    inline
182    void ClearFlag(Flag f) { m_flag &= (~(int)f);}
183
184    inline
185    void fwd( Complex * dst, const Scalar * src, Index nfft)
186    {
187        m_impl.fwd(dst,src,static_cast<int>(nfft));
188        if ( HasFlag(HalfSpectrum) == false)
189          ReflectSpectrum(dst,nfft);
190    }
191
192    inline
193    void fwd( Complex * dst, const Complex * src, Index nfft)
194    {
195        m_impl.fwd(dst,src,static_cast<int>(nfft));
196    }
197
198    /*
199    inline
200    void fwd2(Complex * dst, const Complex * src, int n0,int n1)
201    {
202      m_impl.fwd2(dst,src,n0,n1);
203    }
204    */
205
206    template <typename _Input>
207    inline
208    void fwd( std::vector<Complex> & dst, const std::vector<_Input> & src)
209    {
210      if ( NumTraits<_Input>::IsComplex == 0 && HasFlag(HalfSpectrum) )
211        dst.resize( (src.size()>>1)+1); // half the bins + Nyquist bin
212      else
213        dst.resize(src.size());
214      fwd(&dst[0],&src[0],src.size());
215    }
216
217    template<typename InputDerived, typename ComplexDerived>
218    inline
219    void fwd( MatrixBase<ComplexDerived> & dst, const MatrixBase<InputDerived> & src, Index nfft=-1)
220    {
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)
230
231      if (nfft<1)
232        nfft = src.size();
233
234      if ( NumTraits< src_type >::IsComplex == 0 && HasFlag(HalfSpectrum) )
235        dst.derived().resize( (nfft>>1)+1);
236      else
237        dst.derived().resize(nfft);
238
239      if ( src.innerStride() != 1 || src.size() < nfft ) {
240        Matrix<src_type,1,Dynamic> tmp;
241        if (src.size()<nfft) {
242          tmp.setZero(nfft);
243          tmp.block(0,0,src.size(),1 ) = src;
244        }else{
245          tmp = src;
246        }
247        fwd( &dst[0],&tmp[0],nfft );
248      }else{
249        fwd( &dst[0],&src[0],nfft );
250      }
251    }
252
253    template<typename InputDerived>
254    inline
255    fft_fwd_proxy< MatrixBase<InputDerived>, FFT<T_Scalar,T_Impl> >
256    fwd( const MatrixBase<InputDerived> & src, Index nfft=-1)
257    {
258      return fft_fwd_proxy< MatrixBase<InputDerived> ,FFT<T_Scalar,T_Impl> >( src, *this,nfft );
259    }
260
261    template<typename InputDerived>
262    inline
263    fft_inv_proxy< MatrixBase<InputDerived>, FFT<T_Scalar,T_Impl> >
264    inv( const MatrixBase<InputDerived> & src, Index nfft=-1)
265    {
266      return  fft_inv_proxy< MatrixBase<InputDerived> ,FFT<T_Scalar,T_Impl> >( src, *this,nfft );
267    }
268
269    inline
270    void inv( Complex * dst, const Complex * src, Index nfft)
271    {
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
275    }
276
277    inline
278    void inv( Scalar * dst, const Complex * src, Index nfft)
279    {
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
283    }
284
285    template<typename OutputDerived, typename ComplexDerived>
286    inline
287    void inv( MatrixBase<OutputDerived> & dst, const MatrixBase<ComplexDerived> & src, Index nfft=-1)
288    {
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)
300
301      if (nfft<1) { //automatic FFT size determination
302        if ( realfft && HasFlag(HalfSpectrum) )
303          nfft = 2*(src.size()-1); //assume even fft size
304        else
305          nfft = src.size();
306      }
307      dst.derived().resize( nfft );
308
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() );
313
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
324          }else{
325            size_t nhead,ntail;
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);
335            }
336          }
337        }else{
338          tmp = src;
339        }
340        inv( &dst[0],&tmp[0], nfft);
341      }else{
342        inv( &dst[0],&src[0], nfft);
343      }
344    }
345
346    template <typename _Output>
347    inline
348    void inv( std::vector<_Output> & dst, const std::vector<Complex> & src,Index nfft=-1)
349    {
350      if (nfft<1)
351        nfft = ( NumTraits<_Output>::IsComplex == 0 && HasFlag(HalfSpectrum) ) ? 2*(src.size()-1) : src.size();
352      dst.resize( nfft );
353      inv( &dst[0],&src[0],nfft);
354    }
355
356
357    /*
358    // TODO: multi-dimensional FFTs
359    inline
360    void inv2(Complex * dst, const Complex * src, int n0,int n1)
361    {
362      m_impl.inv2(dst,src,n0,n1);
363      if ( HasFlag( Unscaled ) == false)
364          scale(dst,1./(n0*n1),n0*n1);
365    }
366  */
367
368    inline
369    impl_type & impl() {return m_impl;}
370  private:
371
372    template <typename T_Data>
373    inline
374    void scale(T_Data * x,Scalar s,Index nx)
375    {
376#if 1
377      for (int k=0;k<nx;++k)
378        *x++ *= s;
379#else
380      if ( ((ptrdiff_t)x) & 15 )
381        Matrix<T_Data, Dynamic, 1>::Map(x,nx) *= s;
382      else
383        Matrix<T_Data, Dynamic, 1>::MapAligned(x,nx) *= s;
384         //Matrix<T_Data, Dynamic, Dynamic>::Map(x,nx) * s;
385#endif
386    }
387
388    inline
389    void ReflectSpectrum(Complex * freq, Index nfft)
390    {
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]);
395    }
396
397    impl_type m_impl;
398    int m_flag;
399};
400
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
404{
405    m_ifc.fwd( dst, m_src, m_nfft);
406}
407
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
411{
412    m_ifc.inv( dst, m_src, m_nfft);
413}
414
415}
416
417#include "../../Eigen/src/Core/util/ReenableStupidWarnings.h"
418
419#endif
420