mirror of
https://gitlab.com/libeigen/eigen.git
synced 2025-07-23 21:34:30 +08:00
262 lines
9.0 KiB
C++
262 lines
9.0 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
|
|
//
|
|
// This Source Code Form is subject to the terms of the Mozilla
|
|
// Public License v. 2.0. If a copy of the MPL was not distributed
|
|
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
|
|
|
|
namespace Eigen {
|
|
|
|
namespace internal {
|
|
|
|
// FFTW uses non-const arguments
|
|
// so we must use ugly const_cast calls for all the args it uses
|
|
//
|
|
// This should be safe as long as
|
|
// 1. we use FFTW_ESTIMATE for all our planning
|
|
// see the FFTW docs section 4.3.2 "Planner Flags"
|
|
// 2. fftw_complex is compatible with std::complex
|
|
// This assumes std::complex<T> layout is array of size 2 with real,imag
|
|
template <typename T>
|
|
inline
|
|
T * fftw_cast(const T* p)
|
|
{
|
|
return const_cast<T*>( p);
|
|
}
|
|
|
|
inline
|
|
fftw_complex * fftw_cast( const std::complex<double> * p)
|
|
{
|
|
return const_cast<fftw_complex*>( reinterpret_cast<const fftw_complex*>(p) );
|
|
}
|
|
|
|
inline
|
|
fftwf_complex * fftw_cast( const std::complex<float> * p)
|
|
{
|
|
return const_cast<fftwf_complex*>( reinterpret_cast<const fftwf_complex*>(p) );
|
|
}
|
|
|
|
inline
|
|
fftwl_complex * fftw_cast( const std::complex<long double> * p)
|
|
{
|
|
return const_cast<fftwl_complex*>( reinterpret_cast<const fftwl_complex*>(p) );
|
|
}
|
|
|
|
template <typename T>
|
|
struct fftw_plan {};
|
|
|
|
template <>
|
|
struct fftw_plan<float>
|
|
{
|
|
typedef float scalar_type;
|
|
typedef fftwf_complex complex_type;
|
|
fftwf_plan m_plan;
|
|
fftw_plan() :m_plan(NULL) {}
|
|
~fftw_plan() {if (m_plan) fftwf_destroy_plan(m_plan);}
|
|
|
|
inline
|
|
void fwd(complex_type * dst,complex_type * src,int nfft) {
|
|
if (m_plan==NULL) m_plan = fftwf_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
|
|
fftwf_execute_dft( m_plan, src,dst);
|
|
}
|
|
inline
|
|
void inv(complex_type * dst,complex_type * src,int nfft) {
|
|
if (m_plan==NULL) m_plan = fftwf_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
|
|
fftwf_execute_dft( m_plan, src,dst);
|
|
}
|
|
inline
|
|
void fwd(complex_type * dst,scalar_type * src,int nfft) {
|
|
if (m_plan==NULL) m_plan = fftwf_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
|
|
fftwf_execute_dft_r2c( m_plan,src,dst);
|
|
}
|
|
inline
|
|
void inv(scalar_type * dst,complex_type * src,int nfft) {
|
|
if (m_plan==NULL)
|
|
m_plan = fftwf_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
|
|
fftwf_execute_dft_c2r( m_plan, src,dst);
|
|
}
|
|
|
|
inline
|
|
void fwd2( complex_type * dst,complex_type * src,int n0,int n1) {
|
|
if (m_plan==NULL) m_plan = fftwf_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
|
|
fftwf_execute_dft( m_plan, src,dst);
|
|
}
|
|
inline
|
|
void inv2( complex_type * dst,complex_type * src,int n0,int n1) {
|
|
if (m_plan==NULL) m_plan = fftwf_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
|
|
fftwf_execute_dft( m_plan, src,dst);
|
|
}
|
|
|
|
};
|
|
template <>
|
|
struct fftw_plan<double>
|
|
{
|
|
typedef double scalar_type;
|
|
typedef fftw_complex complex_type;
|
|
::fftw_plan m_plan;
|
|
fftw_plan() :m_plan(NULL) {}
|
|
~fftw_plan() {if (m_plan) fftw_destroy_plan(m_plan);}
|
|
|
|
inline
|
|
void fwd(complex_type * dst,complex_type * src,int nfft) {
|
|
if (m_plan==NULL) m_plan = fftw_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
|
|
fftw_execute_dft( m_plan, src,dst);
|
|
}
|
|
inline
|
|
void inv(complex_type * dst,complex_type * src,int nfft) {
|
|
if (m_plan==NULL) m_plan = fftw_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
|
|
fftw_execute_dft( m_plan, src,dst);
|
|
}
|
|
inline
|
|
void fwd(complex_type * dst,scalar_type * src,int nfft) {
|
|
if (m_plan==NULL) m_plan = fftw_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
|
|
fftw_execute_dft_r2c( m_plan,src,dst);
|
|
}
|
|
inline
|
|
void inv(scalar_type * dst,complex_type * src,int nfft) {
|
|
if (m_plan==NULL)
|
|
m_plan = fftw_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
|
|
fftw_execute_dft_c2r( m_plan, src,dst);
|
|
}
|
|
inline
|
|
void fwd2( complex_type * dst,complex_type * src,int n0,int n1) {
|
|
if (m_plan==NULL) m_plan = fftw_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
|
|
fftw_execute_dft( m_plan, src,dst);
|
|
}
|
|
inline
|
|
void inv2( complex_type * dst,complex_type * src,int n0,int n1) {
|
|
if (m_plan==NULL) m_plan = fftw_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
|
|
fftw_execute_dft( m_plan, src,dst);
|
|
}
|
|
};
|
|
template <>
|
|
struct fftw_plan<long double>
|
|
{
|
|
typedef long double scalar_type;
|
|
typedef fftwl_complex complex_type;
|
|
fftwl_plan m_plan;
|
|
fftw_plan() :m_plan(NULL) {}
|
|
~fftw_plan() {if (m_plan) fftwl_destroy_plan(m_plan);}
|
|
|
|
inline
|
|
void fwd(complex_type * dst,complex_type * src,int nfft) {
|
|
if (m_plan==NULL) m_plan = fftwl_plan_dft_1d(nfft,src,dst, FFTW_FORWARD, FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
|
|
fftwl_execute_dft( m_plan, src,dst);
|
|
}
|
|
inline
|
|
void inv(complex_type * dst,complex_type * src,int nfft) {
|
|
if (m_plan==NULL) m_plan = fftwl_plan_dft_1d(nfft,src,dst, FFTW_BACKWARD , FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
|
|
fftwl_execute_dft( m_plan, src,dst);
|
|
}
|
|
inline
|
|
void fwd(complex_type * dst,scalar_type * src,int nfft) {
|
|
if (m_plan==NULL) m_plan = fftwl_plan_dft_r2c_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
|
|
fftwl_execute_dft_r2c( m_plan,src,dst);
|
|
}
|
|
inline
|
|
void inv(scalar_type * dst,complex_type * src,int nfft) {
|
|
if (m_plan==NULL)
|
|
m_plan = fftwl_plan_dft_c2r_1d(nfft,src,dst,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
|
|
fftwl_execute_dft_c2r( m_plan, src,dst);
|
|
}
|
|
inline
|
|
void fwd2( complex_type * dst,complex_type * src,int n0,int n1) {
|
|
if (m_plan==NULL) m_plan = fftwl_plan_dft_2d(n0,n1,src,dst,FFTW_FORWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
|
|
fftwl_execute_dft( m_plan, src,dst);
|
|
}
|
|
inline
|
|
void inv2( complex_type * dst,complex_type * src,int n0,int n1) {
|
|
if (m_plan==NULL) m_plan = fftwl_plan_dft_2d(n0,n1,src,dst,FFTW_BACKWARD,FFTW_ESTIMATE|FFTW_PRESERVE_INPUT);
|
|
fftwl_execute_dft( m_plan, src,dst);
|
|
}
|
|
};
|
|
|
|
template <typename _Scalar>
|
|
struct fftw_impl
|
|
{
|
|
typedef _Scalar Scalar;
|
|
typedef std::complex<Scalar> Complex;
|
|
|
|
inline
|
|
void clear()
|
|
{
|
|
m_plans.clear();
|
|
}
|
|
|
|
// complex-to-complex forward FFT
|
|
inline
|
|
void fwd( Complex * dst,const Complex *src,int nfft)
|
|
{
|
|
get_plan(nfft,false,dst,src).fwd(fftw_cast(dst), fftw_cast(src),nfft );
|
|
}
|
|
|
|
// real-to-complex forward FFT
|
|
inline
|
|
void fwd( Complex * dst,const Scalar * src,int nfft)
|
|
{
|
|
get_plan(nfft,false,dst,src).fwd(fftw_cast(dst), fftw_cast(src) ,nfft);
|
|
}
|
|
|
|
// 2-d complex-to-complex
|
|
inline
|
|
void fwd2(Complex * dst, const Complex * src, int n0,int n1)
|
|
{
|
|
get_plan(n0,n1,false,dst,src).fwd2(fftw_cast(dst), fftw_cast(src) ,n0,n1);
|
|
}
|
|
|
|
// inverse complex-to-complex
|
|
inline
|
|
void inv(Complex * dst,const Complex *src,int nfft)
|
|
{
|
|
get_plan(nfft,true,dst,src).inv(fftw_cast(dst), fftw_cast(src),nfft );
|
|
}
|
|
|
|
// half-complex to scalar
|
|
inline
|
|
void inv( Scalar * dst,const Complex * src,int nfft)
|
|
{
|
|
get_plan(nfft,true,dst,src).inv(fftw_cast(dst), fftw_cast(src),nfft );
|
|
}
|
|
|
|
// 2-d complex-to-complex
|
|
inline
|
|
void inv2(Complex * dst, const Complex * src, int n0,int n1)
|
|
{
|
|
get_plan(n0,n1,true,dst,src).inv2(fftw_cast(dst), fftw_cast(src) ,n0,n1);
|
|
}
|
|
|
|
|
|
protected:
|
|
typedef fftw_plan<Scalar> PlanData;
|
|
|
|
typedef Eigen::numext::int64_t int64_t;
|
|
|
|
typedef std::map<int64_t,PlanData> PlanMap;
|
|
|
|
PlanMap m_plans;
|
|
|
|
inline
|
|
PlanData & get_plan(int nfft,bool inverse,void * dst,const void * src)
|
|
{
|
|
bool inplace = (dst==src);
|
|
bool aligned = ( (reinterpret_cast<size_t>(src)&15) | (reinterpret_cast<size_t>(dst)&15) ) == 0;
|
|
int64_t key = ( (nfft<<3 ) | (inverse<<2) | (inplace<<1) | aligned ) << 1;
|
|
return m_plans[key];
|
|
}
|
|
|
|
inline
|
|
PlanData & get_plan(int n0,int n1,bool inverse,void * dst,const void * src)
|
|
{
|
|
bool inplace = (dst==src);
|
|
bool aligned = ( (reinterpret_cast<size_t>(src)&15) | (reinterpret_cast<size_t>(dst)&15) ) == 0;
|
|
int64_t key = ( ( (((int64_t)n0) << 30)|(n1<<3 ) | (inverse<<2) | (inplace<<1) | aligned ) << 1 ) + 1;
|
|
return m_plans[key];
|
|
}
|
|
};
|
|
|
|
} // end namespace internal
|
|
|
|
} // end namespace Eigen
|