mirror/eigen
2
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-01-30 17:40:05 +08:00
Code Issues Packages Projects Releases Wiki Activity

initial pass of FFT module -- includes complex 1-d case only

Browse Source
This commit is contained in:
Mark Borgerding 2009-05-19 00:21:04 -04:00
parent 72f66d717d
commit 92ca9fc032
4 changed files with 466 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

84
unsupported/Eigen/FFT.h Normal file
View File

@ -0,0 +1,84 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra. Eigen itself is part of the KDE project.
//
// 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
// simple_fft_traits: small, free, reasonably efficient default, derived from kissfft
#include "src/FFT/simple_fft_traits.h"
#define DEFAULT_FFT_TRAITS simple_fft_traits
// FFTW: faster, GPL-not LGPL, 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/fftw_traits.h"
// #define DEFAULT_FFT_TRAITS fftw_traits
#endif
// intel Math Kernel Library: fastest, commerical
#ifdef _MKL_DFTI_H_ // mkl_dfti.h has been included, we can use MKL FFT routines
// TODO
// #include "src/FFT/imkl_traits.h"
// #define DEFAULT_FFT_TRAITS imkl_traits
#endif
namespace Eigen {
template <typename _Scalar,
typename _Traits=DEFAULT_FFT_TRAITS<_Scalar>
>
class FFT
{
public:
typedef _Traits traits_type;
typedef typename traits_type::Scalar Scalar;
typedef typename traits_type::Complex Complex;
FFT(const traits_type & traits=traits_type() ) :m_traits(traits) { }
void fwd( Complex * dst, const Complex * src, int nfft)
{
m_traits.prepare(nfft,false,dst,src);
m_traits.exec(dst,src);
m_traits.postprocess(dst);
}
void inv( Complex * dst, const Complex * src, int nfft)
{
m_traits.prepare(nfft,true,dst,src);
m_traits.exec(dst,src);
m_traits.postprocess(dst);
}
// TODO: fwd,inv for Scalar
// TODO: multi-dimensional FFTs
// TODO: handle Eigen MatrixBase
traits_type & traits() {return m_traits;}
private:
traits_type m_traits;
};
#undef DEFAULT_FFT_TRAITS
}
#endif

296
unsupported/Eigen/src/FFT/simple_fft_traits.h Normal file
View File

@ -0,0 +1,296 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra. Eigen itself is part of the KDE project.
//
// 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/>.
#include <complex>
#include <vector>
namespace Eigen {
template <typename _Scalar>
struct simple_fft_traits
{
typedef _Scalar Scalar;
typedef std::complex<Scalar> Complex;
simple_fft_traits() : m_nfft(0) {}
void prepare(int nfft,bool inverse,Complex * dst,const Complex *src)
{
if (m_nfft == nfft) {
// reuse the twiddles, conjugate if necessary
if (inverse != m_inverse)
for (int i=0;i<nfft;++i)
m_twiddles[i] = conj( m_twiddles[i] );
m_inverse = inverse;
return;
}
m_nfft = nfft;
m_inverse = inverse;
m_twiddles.resize(nfft);
Scalar phinc = (inverse?2:-2)* acos( (Scalar) -1) / nfft;
for (int i=0;i<nfft;++i)
m_twiddles[i] = exp( Complex(0,i*phinc) );
m_stageRadix.resize(0);
m_stageRemainder.resize(0);
//factorize
//start factoring out 4's, then 2's, then 3,5,7,9,...
int n= nfft;
int p=4;
do {
while (n % p) {
switch (p) {
case 4: p = 2; break;
case 2: p = 3; break;
default: p += 2; break;
}
if (p*p>n)
p=n;// no more factors
}
n /= p;
m_stageRadix.push_back(p);
m_stageRemainder.push_back(n);
}while(n>1);
}
void exec(Complex * dst, const Complex * src)
{
work(0, dst, src, 1,1);
}
void postprocess(Complex *dst)
{
if (m_inverse) {
Scalar scale = 1./m_nfft;
for (int k=0;k<m_nfft;++k)
dst[k] *= scale;
}
}
private:
void work( int stage,Complex * Fout, const Complex * f, size_t fstride,size_t in_stride)
{
int p = m_stageRadix[stage];
int m = m_stageRemainder[stage];
Complex * Fout_beg = Fout;
Complex * Fout_end = Fout + p*m;
if (m==1) {
do{
*Fout = *f;
f += fstride*in_stride;
}while(++Fout != Fout_end );
}else{
do{
// recursive call:
// DFT of size m*p performed by doing
// p instances of smaller DFTs of size m,
// each one takes a decimated version of the input
work(stage+1, Fout , f, fstride*p,in_stride);
f += fstride*in_stride;
}while( (Fout += m) != Fout_end );
}
Fout=Fout_beg;
// recombine the p smaller DFTs
switch (p) {
case 2: bfly2(Fout,fstride,m); break;
case 3: bfly3(Fout,fstride,m); break;
case 4: bfly4(Fout,fstride,m); break;
case 5: bfly5(Fout,fstride,m); break;
default: bfly_generic(Fout,fstride,m,p); break;
}
}
void bfly2( Complex * Fout, const size_t fstride, int m)
{
for (int k=0;k<m;++k) {
Complex t = Fout[m+k] * m_twiddles[k*fstride];
Fout[m+k] = Fout[k] - t;
Fout[k] += t;
}
}
void bfly4( Complex * Fout, const size_t fstride, const size_t m)
{
Complex scratch[7];
int negative_if_inverse = m_inverse * -2 +1;
for (size_t k=0;k<m;++k) {
scratch[0] = Fout[k+m] * m_twiddles[k*fstride];
scratch[1] = Fout[k+2*m] * m_twiddles[k*fstride*2];
scratch[2] = Fout[k+3*m] * m_twiddles[k*fstride*3];
scratch[5] = Fout[k] - scratch[1];
Fout[k] += scratch[1];
scratch[3] = scratch[0] + scratch[2];
scratch[4] = scratch[0] - scratch[2];
scratch[4] = Complex( scratch[4].imag()*negative_if_inverse , -scratch[4].real()* negative_if_inverse );
Fout[k+2*m] = Fout[k] - scratch[3];
Fout[k] += scratch[3];
Fout[k+m] = scratch[5] + scratch[4];
Fout[k+3*m] = scratch[5] - scratch[4];
}
}
void bfly3( Complex * Fout, const size_t fstride, const size_t m)
{
size_t k=m;
const size_t m2 = 2*m;
Complex *tw1,*tw2;
Complex scratch[5];
Complex epi3;
epi3 = m_twiddles[fstride*m];
tw1=tw2=&m_twiddles[0];
do{
scratch[1]=Fout[m] * *tw1;
scratch[2]=Fout[m2] * *tw2;
scratch[3]=scratch[1]+scratch[2];
scratch[0]=scratch[1]-scratch[2];
tw1 += fstride;
tw2 += fstride*2;
Fout[m] = Complex( Fout->real() - .5*scratch[3].real() , Fout->imag() - .5*scratch[3].imag() );
scratch[0] *= epi3.imag();
*Fout += scratch[3];
Fout[m2] = Complex( Fout[m].real() + scratch[0].imag() , Fout[m].imag() - scratch[0].real() );
Fout[m] += Complex( -scratch[0].imag(),scratch[0].real() );
++Fout;
}while(--k);
}
void bfly5( Complex * Fout, const size_t fstride, const size_t m)
{
Complex *Fout0,*Fout1,*Fout2,*Fout3,*Fout4;
size_t u;
Complex scratch[13];
Complex * twiddles = &m_twiddles[0];
Complex *tw;
Complex ya,yb;
ya = twiddles[fstride*m];
yb = twiddles[fstride*2*m];
Fout0=Fout;
Fout1=Fout0+m;
Fout2=Fout0+2*m;
Fout3=Fout0+3*m;
Fout4=Fout0+4*m;
tw=twiddles;
for ( u=0; u<m; ++u ) {
scratch[0] = *Fout0;
scratch[1] = *Fout1 * tw[u*fstride];
scratch[2] = *Fout2 * tw[2*u*fstride];
scratch[3] = *Fout3 * tw[3*u*fstride];
scratch[4] = *Fout4 * tw[4*u*fstride];
scratch[7] = scratch[1] + scratch[4];
scratch[10] = scratch[1] - scratch[4];
scratch[8] = scratch[2] + scratch[3];
scratch[9] = scratch[2] - scratch[3];
*Fout0 += scratch[7];
*Fout0 += scratch[8];
scratch[5] = scratch[0] + Complex(
(scratch[7].real()*ya.real() ) + (scratch[8].real() *yb.real() ),
(scratch[7].imag()*ya.real()) + (scratch[8].imag()*yb.real())
);
scratch[6] = Complex(
(scratch[10].imag()*ya.imag()) + (scratch[9].imag()*yb.imag()),
-(scratch[10].real()*ya.imag()) - (scratch[9].real()*yb.imag())
);
*Fout1 = scratch[5] - scratch[6];
*Fout4 = scratch[5] + scratch[6];
scratch[11] = scratch[0] +
Complex(
(scratch[7].real()*yb.real()) + (scratch[8].real()*ya.real()),
(scratch[7].imag()*yb.real()) + (scratch[8].imag()*ya.real())
);
scratch[12] = Complex(
-(scratch[10].imag()*yb.imag()) + (scratch[9].imag()*ya.imag()),
(scratch[10].real()*yb.imag()) - (scratch[9].real()*ya.imag())
);
*Fout2=scratch[11]+scratch[12];
*Fout3=scratch[11]-scratch[12];
++Fout0;++Fout1;++Fout2;++Fout3;++Fout4;
}
}
/* perform the butterfly for one stage of a mixed radix FFT */
void bfly_generic(
Complex * Fout,
const size_t fstride,
int m,
int p
)
{
int u,k,q1,q;
Complex * twiddles = &m_twiddles[0];
Complex t;
int Norig = m_nfft;
Complex * scratchbuf = (Complex*)alloca(p*sizeof(Complex) );
for ( u=0; u<m; ++u ) {
k=u;
for ( q1=0 ; q1<p ; ++q1 ) {
scratchbuf[q1] = Fout[ k ];
k += m;
}
k=u;
for ( q1=0 ; q1<p ; ++q1 ) {
int twidx=0;
Fout[ k ] = scratchbuf[0];
for (q=1;q<p;++q ) {
twidx += fstride * k;
if (twidx>=Norig) twidx-=Norig;
t=scratchbuf[q] * twiddles[twidx];
Fout[ k ] += t;
}
k += m;
}
}
}
int m_nfft;
bool m_inverse;
std::vector<Complex> m_twiddles;
std::vector<int> m_stageRadix;
std::vector<int> m_stageRemainder;
};
}

1
unsupported/test/CMakeLists.txt
View File

@ -19,3 +19,4 @@ ei_add_test(autodiff)
ei_add_test(BVH)
ei_add_test(matrixExponential)
ei_add_test(alignedvector3)
ei_add_test(FFT)

85
unsupported/test/FFT.cpp Normal file
View File

@ -0,0 +1,85 @@
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra. Eigen itself is part of the KDE project.
//
// 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/>.
#include "main.h"
#include <unsupported/Eigen/FFT.h>
//#include <iostream>
//#include <cstdlib>
//#include <typeinfo>
using namespace std;
template <class T>
void test_fft(int nfft)
{
typedef typename Eigen::FFT<T>::Complex Complex;
//cout << "type:" << typeid(T).name() << " nfft:" << nfft;
FFT<T> fft;
vector<Complex> inbuf(nfft);
vector<Complex> buf3(nfft);
vector<Complex> outbuf(nfft);
for (int k=0;k<nfft;++k)
inbuf[k]= Complex(
(T)(rand()/(double)RAND_MAX - .5),
(T)(rand()/(double)RAND_MAX - .5) );
fft.fwd( &outbuf[0] , &inbuf[0] ,nfft);
fft.inv( &buf3[0] , &outbuf[0] ,nfft);
long double totalpower=0;
long double difpower=0;
for (int k0=0;k0<nfft;++k0) {
complex<long double> acc = 0;
long double phinc = 2*k0* M_PIl / nfft;
for (int k1=0;k1<nfft;++k1) {
complex<long double> x(inbuf[k1].real(),inbuf[k1].imag());
acc += x * exp( complex<long double>(0,-k1*phinc) );
}
totalpower += norm(acc);
complex<long double> x(outbuf[k0].real(),outbuf[k0].imag());
complex<long double> dif = acc - x;
difpower += norm(dif);
}
long double rmse = sqrt(difpower/totalpower);
VERIFY( rmse < 1e-5 );// gross check
totalpower=0;
difpower=0;
for (int k=0;k<nfft;++k) {
totalpower += norm( inbuf[k] );
difpower += norm(inbuf[k] - buf3[k]);
}
rmse = sqrt(difpower/totalpower);
VERIFY( rmse < 1e-5 );// gross check
}
void test_FFT()
{
CALL_SUBTEST(( test_fft<float>(32) )); CALL_SUBTEST(( test_fft<double>(32) )); CALL_SUBTEST(( test_fft<long double>(32) ));
CALL_SUBTEST(( test_fft<float>(1024) )); CALL_SUBTEST(( test_fft<double>(1024) )); CALL_SUBTEST(( test_fft<long double>(1024) ));
CALL_SUBTEST(( test_fft<float>(2*3*4*5*7) )); CALL_SUBTEST(( test_fft<double>(2*3*4*5*7) )); CALL_SUBTEST(( test_fft<long double>(2*3*4*5*7) ));
}