mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-03-01 18:26:24 +08:00
Code Issues Packages Projects Releases Wiki Activity

added ei_fftw_impl

Browse Source
This commit is contained in:
Mark Borgerding 2009-05-30 17:55:47 -04:00
parent f13e000b45
commit 1fd6dfe428
3 changed files with 214 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
unsupported/Eigen/FFT
View File

@ -30,9 +30,8 @@
#define DEFAULT_FFT_IMPL ei_kissfft_impl
// FFTW: faster, GPL -- incompatible with Eigen in LGPL form, bigger code size
#ifdef FFTW_PATIENT // definition of FFTW_PATIENT indicates the caller has included fftw3.h, we can use FFTW routines
// TODO
// #include "src/FFT/ei_fftw_impl.h"
#ifdef FFTW_ESTIMATE // definition of FFTW_ESTIMATE indicates the caller has included fftw3.h, we can use FFTW routines
#include "src/FFT/ei_fftw_impl.h"
// #define DEFAULT_FFT_IMPL ei_fftw_impl
#endif

198
unsupported/Eigen/src/FFT/ei_fftw_impl.h Normal file
View File

@ -0,0 +1,198 @@
// 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/>.
namespace Eigen {
// 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>
T * ei_fftw_cast(const T* p)
{
return const_cast<T*>( p);
}
fftw_complex * ei_fftw_cast( const std::complex<double> * p)
{
return const_cast<fftw_complex*>( reinterpret_cast<const fftw_complex*>(p) );
}
fftwf_complex * ei_fftw_cast( const std::complex<float> * p)
{
return const_cast<fftwf_complex*>( reinterpret_cast<const fftwf_complex*>(p) );
}
fftwl_complex * ei_fftw_cast( const std::complex<long double> * p)
{
return const_cast<fftwl_complex*>( reinterpret_cast<const fftwl_complex*>(p) );
}
template <typename T>
struct ei_fftw_plan {};
template <>
struct ei_fftw_plan<float>
{
typedef float scalar_type;
typedef fftwf_complex complex_type;
fftwf_plan m_plan;
ei_fftw_plan() :m_plan(NULL) {}
~ei_fftw_plan() {if (m_plan) fftwf_destroy_plan(m_plan);}
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);
fftwf_execute_dft( m_plan, src,dst);
}
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);
fftwf_execute_dft( m_plan, src,dst);
}
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);
fftwf_execute_dft_r2c( m_plan,src,dst);
}
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);
fftwf_execute_dft_c2r( m_plan, src,dst);
}
};
template <>
struct ei_fftw_plan<double>
{
typedef double scalar_type;
typedef fftw_complex complex_type;
fftw_plan m_plan;
ei_fftw_plan() :m_plan(NULL) {}
~ei_fftw_plan() {if (m_plan) fftw_destroy_plan(m_plan);}
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_execute_dft( m_plan, src,dst);
}
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_execute_dft( m_plan, src,dst);
}
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_execute_dft_r2c( m_plan,src,dst);
}
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_execute_dft_c2r( m_plan, src,dst);
}
};
template <>
struct ei_fftw_plan<long double>
{
typedef long double scalar_type;
typedef fftwl_complex complex_type;
fftwl_plan m_plan;
ei_fftw_plan() :m_plan(NULL) {}
~ei_fftw_plan() {if (m_plan) fftwl_destroy_plan(m_plan);}
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);
fftwl_execute_dft( m_plan, src,dst);
}
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);
fftwl_execute_dft( m_plan, src,dst);
}
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);
fftwl_execute_dft_r2c( m_plan,src,dst);
}
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);
fftwl_execute_dft_c2r( m_plan, src,dst);
}
};
template <typename _Scalar>
struct ei_fftw_impl
{
typedef _Scalar Scalar;
typedef std::complex<Scalar> Complex;
void clear()
{
m_plans.clear();
}
void fwd( Complex * dst,const Complex *src,int nfft)
{
get_plan(nfft,false,dst,src).fwd(ei_fftw_cast(dst), ei_fftw_cast(src),nfft );
}
// real-to-complex forward FFT
void fwd( Complex * dst,const Scalar * src,int nfft)
{
get_plan(nfft,false,dst,src).fwd(ei_fftw_cast(dst), ei_fftw_cast(src) ,nfft);
int nhbins=(nfft>>1)+1;
for (int k=nhbins;k < nfft; ++k )
dst[k] = conj(dst[nfft-k]);
}
// inverse complex-to-complex
void inv(Complex * dst,const Complex *src,int nfft)
{
get_plan(nfft,true,dst,src).inv(ei_fftw_cast(dst), ei_fftw_cast(src),nfft );
// scaling
Scalar s = 1./nfft;
for (int k=0;k<nfft;++k)
dst[k] *= s;
}
// half-complex to scalar
void inv( Scalar * dst,const Complex * src,int nfft)
{
get_plan(nfft,true,dst,src).inv(ei_fftw_cast(dst), ei_fftw_cast(src),nfft );
Scalar s = 1./nfft;
for (int k=0;k<nfft;++k)
dst[k] *= s;
}
private:
typedef ei_fftw_plan<Scalar> PlanData;
typedef std::map<int,PlanData> PlanMap;
PlanMap m_plans;
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;
int key = (nfft<<3 ) | (inverse<<2) | (inplace<<1) | aligned;
return m_plans[key];
}
};
}

25
unsupported/Eigen/src/FFT/ei_kissfft_impl.h
View File

@ -39,6 +39,7 @@ namespace Eigen {
std::vector<Complex> m_twiddles;
std::vector<int> m_stageRadix;
std::vector<int> m_stageRemainder;
std::vector<Complex> m_scratchBuf;
bool m_inverse;
void make_twiddles(int nfft,bool inverse)
@ -75,6 +76,8 @@ namespace Eigen {
n /= p;
m_stageRadix.push_back(p);
m_stageRemainder.push_back(n);
if ( p > 5 )
m_scratchBuf.resize(p); // scratchbuf will be needed in bfly_generic
}while(n>1);
}
@ -249,7 +252,7 @@ namespace Eigen {
Complex * twiddles = &m_twiddles[0];
Complex t;
int Norig = m_twiddles.size();
Complex * scratchbuf = (Complex*)alloca(p*sizeof(Complex) );
Complex * scratchbuf = &m_scratchBuf[0];
for ( u=0; u<m; ++u ) {
k=u;
@ -341,28 +344,28 @@ namespace Eigen {
void inv( Scalar * dst,const Complex * src,int nfft)
{
if (nfft&3) {
m_scratchBuf.resize(nfft);
inv(&m_scratchBuf[0],src,nfft);
m_tmpBuf.resize(nfft);
inv(&m_tmpBuf[0],src,nfft);
for (int k=0;k<nfft;++k)
dst[k] = m_scratchBuf[k].real();
dst[k] = m_tmpBuf[k].real();
}else{
// optimized version for multiple of 4
int ncfft = nfft>>1;
int ncfft2 = nfft>>2;
Complex * rtw = real_twiddles(ncfft2);
m_scratchBuf.resize(ncfft);
m_scratchBuf[0] = Complex( src[0].real() + src[ncfft].real(), src[0].real() - src[ncfft].real() );
m_tmpBuf.resize(ncfft);
m_tmpBuf[0] = Complex( src[0].real() + src[ncfft].real(), src[0].real() - src[ncfft].real() );
for (int k = 1; k <= ncfft / 2; ++k) {
Complex fk = src[k];
Complex fnkc = conj(src[ncfft-k]);
Complex fek = fk + fnkc;
Complex tmp = fk - fnkc;
Complex fok = tmp * conj(rtw[k-1]);
m_scratchBuf[k] = fek + fok;
m_scratchBuf[ncfft-k] = conj(fek - fok);
m_tmpBuf[k] = fek + fok;
m_tmpBuf[ncfft-k] = conj(fek - fok);
}
scale(&m_scratchBuf[0], ncfft, Scalar(1)/nfft );
get_plan(ncfft,true).work(0, reinterpret_cast<Complex*>(dst), &m_scratchBuf[0], 1,1);
scale(&m_tmpBuf[0], ncfft, Scalar(1)/nfft );
get_plan(ncfft,true).work(0, reinterpret_cast<Complex*>(dst), &m_tmpBuf[0], 1,1);
}
}
@ -372,7 +375,7 @@ namespace Eigen {
PlanMap m_plans;
std::map<int, std::vector<Complex> > m_realTwiddles;
std::vector<Complex> m_scratchBuf;
std::vector<Complex> m_tmpBuf;
int PlanKey(int nfft,bool isinverse) const { return (nfft<<1) | isinverse; }