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

added FFT inverse complex-to-scalar interface (not yet optimized)

Browse Source
This commit is contained in:
Mark Borgerding 2009-05-23 22:50:07 -04:00
parent 3047988172
commit 326ea77390
3 changed files with 81 additions and 44 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

73
bench/benchFFT.cpp
View File

@ -31,34 +31,67 @@
using namespace Eigen;
using namespace std;
#ifndef NFFT
#define NFFT 1024
#endif
template <typename T>
string nameof();
template <> string nameof<float>() {return "float";}
template <> string nameof<double>() {return "double";}
template <> string nameof<long double>() {return "long double";}
#ifndef TYPE
#define TYPE float
#endif
#ifndef NITS
#define NITS (10000000/NFFT)
#ifndef NFFT
#define NFFT 1024
#endif
#ifndef NDATA
#define NDATA 1000000
#endif
int main()
using namespace Eigen;
template <typename T>
void bench(int nfft)
{
vector<complex<TYPE> > inbuf(NFFT);
vector<complex<TYPE> > outbuf(NFFT);
Eigen::FFT<TYPE> fft;
typedef typename NumTraits<T>::Real Scalar;
typedef typename std::complex<Scalar> Complex;
int nits = NDATA/nfft;
vector<T> inbuf(nfft);
vector<Complex > outbuf(nfft);
FFT< Scalar > fft;
fft.fwd( outbuf , inbuf);
fft.fwd( outbuf , inbuf);
BenchTimer timer;
timer.reset();
for (int k=0;k<8;++k) {
timer.start();
for(int i = 0; i < NITS; i++)
fft.fwd( outbuf , inbuf);
timer.stop();
}
double mflops = 5.*NFFT*log2((double)NFFT) / (1e6 * timer.value() / (double)NITS );
cout << "NFFT=" << NFFT << " " << (double(1e-6*NFFT*NITS)/timer.value()) << " MS/s " << mflops << "MFLOPS\n";
BenchTimer timer;
timer.reset();
for (int k=0;k<8;++k) {
timer.start();
for(int i = 0; i < nits; i++)
fft.fwd( outbuf , inbuf);
timer.stop();
}
cout << nameof<Scalar>() << " ";
double mflops = 5.*nfft*log2((double)nfft) / (1e6 * timer.value() / (double)nits );
if ( NumTraits<T>::IsComplex ) {
cout << "complex";
}else{
cout << "real ";
mflops /= 2;
}
cout << " NFFT=" << nfft << " " << (double(1e-6*nfft*nits)/timer.value()) << " MS/s " << mflops << "MFLOPS\n";
}
int main(int argc,char ** argv)
{
bench<complex<float> >(NFFT);
bench<float>(NFFT);
bench<complex<double> >(NFFT);
bench<double>(NFFT);
bench<complex<long double> >(NFFT);
bench<long double>(NFFT);
return 0;
}

31
unsupported/Eigen/src/FFT/simple_fft_traits.h
View File

@ -54,28 +54,14 @@ namespace Eigen {
work(0, dst, src, 1,1);
}else{
int ncfft = nfft>>1;
int ncfft2 = nfft>>2;
// use optimized mode for even real
fwd( dst, reinterpret_cast<const Complex*> (src),ncfft);
make_real_twiddles(nfft);
std::cerr << "dst[0] = " << dst[0] << "\n";
Complex dc = dst[0].real() + dst[0].imag();
Complex nyquist = dst[0].real() - dst[0].imag();
int k;
#if 0
using namespace std;
cerr << "desired:\n";
for ( k=1;k <= (ncfft>>1) ; ++k ) {
Complex x = exp( Complex(0,-3.14159265358979323846264338327 * ((double) (k) / ncfft + .5) ) );
cerr << k << " " << x << "angle(x):" << arg(x) << "\n";
}
dc=0;
cerr << "twiddles:\n";
for (k=0;k<ncfft;++k) {
Complex x = m_twiddles[k];
cerr << k << " " << x << "angle(-x):" << arg(-x) << "\n";
}
#endif
for ( k=1;k <= (ncfft>>1) ; ++k ) {
for ( k=1;k <= ncfft2 ; ++k ) {
Complex fpk = dst[k];
Complex fpnk = conj(dst[ncfft-k]);
Complex f1k = fpk + fpnk;
@ -87,7 +73,6 @@ namespace Eigen {
dst[ncfft-k] = conj(f1k -tw)*Scalar(.5);
}
// place conjugate-symmetric half at the end for completeness
// TODO: make this configurable ( opt-out )
for ( k=1;k < ncfft ; ++k )
@ -98,6 +83,16 @@ namespace Eigen {
}
}
// half-complex to scalar
void inv( Scalar * dst,const Complex * src,int nfft)
{
// TODO add optimized version for even numbers
std::vector<Complex> tmp(nfft);
inv(&tmp[0],src,nfft);
for (int k=0;k<nfft;++k)
dst[k] = tmp[k].real();
}
void inv(Complex * dst,const Complex *src,int nfft)
{
prepare(nfft,true);
@ -156,7 +151,7 @@ namespace Eigen {
default: p += 2; break;
}
if (p*p>n)
p=n;// no more factors
p=n;// impossible to have a factor > sqrt(n)
}
n /= p;
m_stageRadix.push_back(p);

21
unsupported/test/FFT.cpp
View File

@ -28,6 +28,10 @@
using namespace std;
float norm(float x) {return x*x;}
double norm(double x) {return x*x;}
long double norm(long double x) {return x*x;}
template < typename T>
complex<long double> promote(complex<T> x) { return complex<long double>(x.real(),x.imag()); }
@ -83,7 +87,11 @@ void test_scalar(int nfft)
for (int k=0;k<nfft;++k)
inbuf[k]= (T)(rand()/(double)RAND_MAX - .5);
fft.fwd( outbuf,inbuf);
VERIFY( fft_rmse(outbuf,inbuf) < 1e-5 );// gross check
VERIFY( fft_rmse(outbuf,inbuf) < test_precision<T>() );// gross check
vector<Scalar> buf3;
fft.inv( buf3 , outbuf);
VERIFY( dif_rmse(inbuf,buf3) < test_precision<T>() );// gross check
}
template <class T>
@ -100,18 +108,18 @@ void test_complex(int nfft)
inbuf[k]= Complex( (T)(rand()/(double)RAND_MAX - .5), (T)(rand()/(double)RAND_MAX - .5) );
fft.fwd( outbuf , inbuf);
VERIFY( fft_rmse(outbuf,inbuf) < 1e-5 );// gross check
VERIFY( fft_rmse(outbuf,inbuf) < test_precision<T>() );// gross check
fft.inv( buf3 , outbuf);
VERIFY( dif_rmse(inbuf,buf3) < 1e-5 );// gross check
VERIFY( dif_rmse(inbuf,buf3) < test_precision<T>() );// gross check
}
void test_FFT()
{
#if 0
#if 1
CALL_SUBTEST( test_complex<float>(32) ); CALL_SUBTEST( test_complex<double>(32) ); CALL_SUBTEST( test_complex<long double>(32) );
CALL_SUBTEST( test_complex<float>(1024) ); CALL_SUBTEST( test_complex<double>(1024) ); CALL_SUBTEST( test_complex<long double>(1024) );
CALL_SUBTEST( test_complex<float>(256) ); CALL_SUBTEST( test_complex<double>(256) ); CALL_SUBTEST( test_complex<long double>(256) );
CALL_SUBTEST( test_complex<float>(3*8) ); CALL_SUBTEST( test_complex<double>(3*8) ); CALL_SUBTEST( test_complex<long double>(3*8) );
CALL_SUBTEST( test_complex<float>(5*32) ); CALL_SUBTEST( test_complex<double>(5*32) ); CALL_SUBTEST( test_complex<long double>(5*32) );
CALL_SUBTEST( test_complex<float>(2*3*4) ); CALL_SUBTEST( test_complex<double>(2*3*4) ); CALL_SUBTEST( test_complex<long double>(2*3*4) );
@ -120,8 +128,9 @@ void test_FFT()
#endif
#if 1
CALL_SUBTEST( test_scalar<float>(45) ); CALL_SUBTEST( test_scalar<double>(45) ); CALL_SUBTEST( test_scalar<long double>(45) );
CALL_SUBTEST( test_scalar<float>(32) ); CALL_SUBTEST( test_scalar<double>(32) ); CALL_SUBTEST( test_scalar<long double>(32) );
CALL_SUBTEST( test_scalar<float>(1024) ); CALL_SUBTEST( test_scalar<double>(1024) ); CALL_SUBTEST( test_scalar<long double>(1024) );
CALL_SUBTEST( test_scalar<float>(256) ); CALL_SUBTEST( test_scalar<double>(256) ); CALL_SUBTEST( test_scalar<long double>(256) );
CALL_SUBTEST( test_scalar<float>(2*3*4*5*7) ); CALL_SUBTEST( test_scalar<double>(2*3*4*5*7) ); CALL_SUBTEST( test_scalar<long double>(2*3*4*5*7) );
#endif
}