mirror of
https://gitlab.com/libeigen/eigen.git
synced 2024-12-21 07:19:46 +08:00
297 lines
10 KiB
C++
297 lines
10 KiB
C++
// This file is part of Eigen, a lightweight C++ template library
|
|
// for linear algebra.
|
|
//
|
|
// Copyright (C) 2009 Mark Borgerding mark a borgerding net
|
|
//
|
|
// Eigen is free software; you can redistribute it and/or
|
|
// modify it under the terms of the GNU Lesser General Public
|
|
// License as published by the Free Software Foundation; either
|
|
// version 3 of the License, or (at your option) any later version.
|
|
//
|
|
// Alternatively, you can redistribute it and/or
|
|
// modify it under the terms of the GNU General Public License as
|
|
// published by the Free Software Foundation; either version 2 of
|
|
// the License, or (at your option) any later version.
|
|
//
|
|
// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
|
|
// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
|
// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
|
|
// GNU General Public License for more details.
|
|
//
|
|
// You should have received a copy of the GNU Lesser General Public
|
|
// License and a copy of the GNU General Public License along with
|
|
// Eigen. If not, see <http://www.gnu.org/licenses/>.
|
|
|
|
#ifndef EIGEN_FFT_H
|
|
#define EIGEN_FFT_H
|
|
|
|
#include <complex>
|
|
#include <vector>
|
|
#include <map>
|
|
#include <Eigen/Core>
|
|
|
|
|
|
/** \ingroup Unsupported_modules
|
|
* \defgroup FFT_Module Fast Fourier Transform module
|
|
*
|
|
* \code
|
|
* #include <unsupported/Eigen/FFT>
|
|
* \endcode
|
|
*
|
|
* This module provides Fast Fourier transformation, with a configurable backend
|
|
* implementation.
|
|
*
|
|
* The default implementation is based on kissfft. It is a small, free, and
|
|
* reasonably efficient default.
|
|
*
|
|
* There are currently two implementation backend:
|
|
*
|
|
* - fftw (http://www.fftw.org) : faster, GPL -- incompatible with Eigen in LGPL form, bigger code size.
|
|
* - MKL (http://en.wikipedia.org/wiki/Math_Kernel_Library) : fastest, commercial -- may be incompatible with Eigen in GPL form.
|
|
*
|
|
* \section FFTDesign Design
|
|
*
|
|
* The following design decisions were made concerning scaling and
|
|
* half-spectrum for real FFT.
|
|
*
|
|
* The intent is to facilitate generic programming and ease migrating code
|
|
* from Matlab/octave.
|
|
* We think the default behavior of Eigen/FFT should favor correctness and
|
|
* generality over speed. Of course, the caller should be able to "opt-out" from this
|
|
* behavior and get the speed increase if they want it.
|
|
*
|
|
* 1) %Scaling:
|
|
* Other libraries (FFTW,IMKL,KISSFFT) do not perform scaling, so there
|
|
* is a constant gain incurred after the forward&inverse transforms , so
|
|
* IFFT(FFT(x)) = Kx; this is done to avoid a vector-by-value multiply.
|
|
* The downside is that algorithms that worked correctly in Matlab/octave
|
|
* don't behave the same way once implemented in C++.
|
|
*
|
|
* How Eigen/FFT differs: invertible scaling is performed so IFFT( FFT(x) ) = x.
|
|
*
|
|
* 2) Real FFT half-spectrum
|
|
* Other libraries use only half the frequency spectrum (plus one extra
|
|
* sample for the Nyquist bin) for a real FFT, the other half is the
|
|
* conjugate-symmetric of the first half. This saves them a copy and some
|
|
* memory. The downside is the caller needs to have special logic for the
|
|
* number of bins in complex vs real.
|
|
*
|
|
* How Eigen/FFT differs: The full spectrum is returned from the forward
|
|
* transform. This facilitates generic template programming by obviating
|
|
* separate specializations for real vs complex. On the inverse
|
|
* transform, only half the spectrum is actually used if the output type is real.
|
|
*/
|
|
|
|
|
|
#ifdef EIGEN_FFTW_DEFAULT
|
|
// FFTW: faster, GPL -- incompatible with Eigen in LGPL form, bigger code size
|
|
# include <fftw3.h>
|
|
namespace Eigen {
|
|
# include "src/FFT/ei_fftw_impl.h"
|
|
//template <typename T> typedef struct ei_fftw_impl default_fft_impl; this does not work
|
|
template <typename T> struct default_fft_impl : public ei_fftw_impl<T> {};
|
|
}
|
|
#elif defined EIGEN_MKL_DEFAULT
|
|
// TODO
|
|
// intel Math Kernel Library: fastest, commercial -- may be incompatible with Eigen in GPL form
|
|
namespace Eigen {
|
|
# include "src/FFT/ei_imklfft_impl.h"
|
|
template <typename T> struct default_fft_impl : public ei_imklfft_impl {};
|
|
}
|
|
#else
|
|
// ei_kissfft_impl: small, free, reasonably efficient default, derived from kissfft
|
|
//
|
|
namespace Eigen {
|
|
# include "src/FFT/ei_kissfft_impl.h"
|
|
template <typename T>
|
|
struct default_fft_impl : public ei_kissfft_impl<T> {};
|
|
}
|
|
#endif
|
|
|
|
namespace Eigen {
|
|
|
|
template <typename _Scalar,
|
|
typename _Impl=default_fft_impl<_Scalar> >
|
|
class FFT
|
|
{
|
|
public:
|
|
typedef _Impl impl_type;
|
|
typedef typename impl_type::Scalar Scalar;
|
|
typedef typename impl_type::Complex Complex;
|
|
|
|
enum Flag {
|
|
Default=0, // goof proof
|
|
Unscaled=1,
|
|
HalfSpectrum=2,
|
|
// SomeOtherSpeedOptimization=4
|
|
Speedy=32767
|
|
};
|
|
|
|
FFT( const impl_type & impl=impl_type() , Flag flags=Default ) :m_impl(impl),m_flag(flags) { }
|
|
|
|
inline
|
|
bool HasFlag(Flag f) const { return (m_flag & (int)f) == f;}
|
|
|
|
inline
|
|
void SetFlag(Flag f) { m_flag |= (int)f;}
|
|
|
|
inline
|
|
void ClearFlag(Flag f) { m_flag &= (~(int)f);}
|
|
|
|
inline
|
|
void fwd( Complex * dst, const Scalar * src, int nfft)
|
|
{
|
|
m_impl.fwd(dst,src,nfft);
|
|
if ( HasFlag(HalfSpectrum) == false)
|
|
ReflectSpectrum(dst,nfft);
|
|
}
|
|
|
|
inline
|
|
void fwd( Complex * dst, const Complex * src, int nfft)
|
|
{
|
|
m_impl.fwd(dst,src,nfft);
|
|
}
|
|
|
|
/*
|
|
inline
|
|
void fwd2(Complex * dst, const Complex * src, int n0,int n1)
|
|
{
|
|
m_impl.fwd2(dst,src,n0,n1);
|
|
}
|
|
*/
|
|
|
|
template <typename _Input>
|
|
inline
|
|
void fwd( std::vector<Complex> & dst, const std::vector<_Input> & src)
|
|
{
|
|
if ( NumTraits<_Input>::IsComplex == 0 && HasFlag(HalfSpectrum) )
|
|
dst.resize( (src.size()>>1)+1); // half the bins + Nyquist bin
|
|
else
|
|
dst.resize(src.size());
|
|
fwd(&dst[0],&src[0],static_cast<int>(src.size()));
|
|
}
|
|
|
|
template<typename InputDerived, typename ComplexDerived>
|
|
inline
|
|
void fwd( MatrixBase<ComplexDerived> & dst, const MatrixBase<InputDerived> & src)
|
|
{
|
|
EIGEN_STATIC_ASSERT_VECTOR_ONLY(InputDerived)
|
|
EIGEN_STATIC_ASSERT_VECTOR_ONLY(ComplexDerived)
|
|
EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(ComplexDerived,InputDerived) // size at compile-time
|
|
EIGEN_STATIC_ASSERT((ei_is_same_type<typename ComplexDerived::Scalar, Complex>::ret),
|
|
YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
|
|
EIGEN_STATIC_ASSERT(int(InputDerived::Flags)&int(ComplexDerived::Flags)&DirectAccessBit,
|
|
THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_WITH_DIRECT_MEMORY_ACCESS_SUCH_AS_MAP_OR_PLAIN_MATRICES)
|
|
|
|
if ( NumTraits< typename InputDerived::Scalar >::IsComplex == 0 && HasFlag(HalfSpectrum) )
|
|
dst.derived().resize( (src.size()>>1)+1);
|
|
else
|
|
dst.derived().resize(src.size());
|
|
|
|
if (src.stride() != 1) {
|
|
Matrix<typename InputDerived::Scalar,1,Dynamic> tmp = src;
|
|
fwd( &dst[0],&tmp[0],src.size() );
|
|
}else{
|
|
fwd( &dst[0],&src[0],src.size() );
|
|
}
|
|
}
|
|
|
|
inline
|
|
void inv( Complex * dst, const Complex * src, int nfft)
|
|
{
|
|
m_impl.inv( dst,src,nfft );
|
|
if ( HasFlag( Unscaled ) == false)
|
|
scale(dst,Scalar(1./nfft),nfft); // scale the time series
|
|
}
|
|
|
|
inline
|
|
void inv( Scalar * dst, const Complex * src, int nfft)
|
|
{
|
|
m_impl.inv( dst,src,nfft );
|
|
if ( HasFlag( Unscaled ) == false)
|
|
scale(dst,Scalar(1./nfft),nfft); // scale the time series
|
|
}
|
|
|
|
template<typename OutputDerived, typename ComplexDerived>
|
|
inline
|
|
void inv( MatrixBase<OutputDerived> & dst, const MatrixBase<ComplexDerived> & src, int nfft=-1)
|
|
{
|
|
EIGEN_STATIC_ASSERT_VECTOR_ONLY(OutputDerived)
|
|
EIGEN_STATIC_ASSERT_VECTOR_ONLY(ComplexDerived)
|
|
EIGEN_STATIC_ASSERT_SAME_VECTOR_SIZE(ComplexDerived,OutputDerived) // size at compile-time
|
|
EIGEN_STATIC_ASSERT((ei_is_same_type<typename ComplexDerived::Scalar, Complex>::ret),
|
|
YOU_MIXED_DIFFERENT_NUMERIC_TYPES__YOU_NEED_TO_USE_THE_CAST_METHOD_OF_MATRIXBASE_TO_CAST_NUMERIC_TYPES_EXPLICITLY)
|
|
EIGEN_STATIC_ASSERT(int(OutputDerived::Flags)&int(ComplexDerived::Flags)&DirectAccessBit,
|
|
THIS_METHOD_IS_ONLY_FOR_EXPRESSIONS_WITH_DIRECT_MEMORY_ACCESS_SUCH_AS_MAP_OR_PLAIN_MATRICES)
|
|
|
|
if (nfft<1) {
|
|
nfft = ( NumTraits<typename OutputDerived::Scalar>::IsComplex == 0 && HasFlag(HalfSpectrum) ) ? 2*(src.size()-1) : src.size();
|
|
}
|
|
dst.derived().resize( nfft );
|
|
if (src.stride() != 1) {
|
|
// if the vector is strided, then we need to copy it to a packed temporary
|
|
Matrix<typename ComplexDerived::Scalar,1,Dynamic> tmp = src;
|
|
inv( &dst[0],&tmp[0], nfft);
|
|
}else{
|
|
inv( &dst[0],&src[0], nfft);
|
|
}
|
|
}
|
|
|
|
template <typename _Output>
|
|
inline
|
|
void inv( std::vector<_Output> & dst, const std::vector<Complex> & src,int nfft=-1)
|
|
{
|
|
if (nfft<1)
|
|
nfft = ( NumTraits<_Output>::IsComplex == 0 && HasFlag(HalfSpectrum) ) ? 2*(src.size()-1) : src.size();
|
|
dst.resize( nfft );
|
|
inv( &dst[0],&src[0],nfft);
|
|
}
|
|
|
|
|
|
/*
|
|
// TODO: multi-dimensional FFTs
|
|
inline
|
|
void inv2(Complex * dst, const Complex * src, int n0,int n1)
|
|
{
|
|
m_impl.inv2(dst,src,n0,n1);
|
|
if ( HasFlag( Unscaled ) == false)
|
|
scale(dst,1./(n0*n1),n0*n1);
|
|
}
|
|
*/
|
|
|
|
inline
|
|
impl_type & impl() {return m_impl;}
|
|
private:
|
|
|
|
template <typename T_Data>
|
|
inline
|
|
void scale(T_Data * x,Scalar s,int nx)
|
|
{
|
|
#if 1
|
|
for (int k=0;k<nx;++k)
|
|
*x++ *= s;
|
|
#else
|
|
if ( ((ptrdiff_t)x) & 15 )
|
|
Matrix<T_Data, Dynamic, 1>::Map(x,nx) *= s;
|
|
else
|
|
Matrix<T_Data, Dynamic, 1>::MapAligned(x,nx) *= s;
|
|
//Matrix<T_Data, Dynamic, Dynamic>::Map(x,nx) * s;
|
|
#endif
|
|
}
|
|
|
|
inline
|
|
void ReflectSpectrum(Complex * freq,int nfft)
|
|
{
|
|
// create the implicit right-half spectrum (conjugate-mirror of the left-half)
|
|
int nhbins=(nfft>>1)+1;
|
|
for (int k=nhbins;k < nfft; ++k )
|
|
freq[k] = conj(freq[nfft-k]);
|
|
}
|
|
|
|
impl_type m_impl;
|
|
int m_flag;
|
|
};
|
|
}
|
|
#endif
|
|
/* vim: set filetype=cpp et sw=2 ts=2 ai: */
|