diff options
author | Mark Borgerding <mark@borgerding.net> | 2009-05-22 22:37:59 -0400 |
---|---|---|
committer | Mark Borgerding <mark@borgerding.net> | 2009-05-22 22:37:59 -0400 |
commit | 8b4afe3debb47bf15ea291a7f2d21d863d546536 (patch) | |
tree | d81260b937fa2e681015e76b561349192cacee24 /unsupported/test/FFT.cpp | |
parent | 68cad98bc935e53102a9432560085b81c5766743 (diff) |
added non-optimized real forward fft (no inverse yet)
Diffstat (limited to 'unsupported/test/FFT.cpp')
-rw-r--r-- | unsupported/test/FFT.cpp | 117 |
1 files changed, 80 insertions, 37 deletions
diff --git a/unsupported/test/FFT.cpp b/unsupported/test/FFT.cpp index 8347bb76b..ef03359e2 100644 --- a/unsupported/test/FFT.cpp +++ b/unsupported/test/FFT.cpp @@ -25,55 +25,98 @@ #include "main.h" #include <unsupported/Eigen/FFT.h> + using namespace std; +template < typename T> +complex<long double> promote(complex<T> x) { return complex<long double>(x.real(),x.imag()); } + +complex<long double> promote(float x) { return complex<long double>( x); } +complex<long double> promote(double x) { return complex<long double>( x); } +complex<long double> promote(long double x) { return complex<long double>( x); } + + + template <typename T1,typename T2> + long double fft_rmse( const vector<T1> & fftbuf,const vector<T2> & timebuf) + { + long double totalpower=0; + long double difpower=0; + for (size_t k0=0;k0<fftbuf.size();++k0) { + complex<long double> acc = 0; + long double phinc = -2.*k0* M_PIl / timebuf.size(); + for (size_t k1=0;k1<timebuf.size();++k1) { + acc += promote( timebuf[k1] ) * exp( complex<long double>(0,k1*phinc) ); + } + totalpower += norm(acc); + complex<long double> x = promote(fftbuf[k0]); + complex<long double> dif = acc - x; + difpower += norm(dif); + cerr << k0 << ":" << acc << " " << x << endl; + } + cerr << "rmse:" << sqrt(difpower/totalpower) << endl; + return sqrt(difpower/totalpower); + } + + template <typename T1,typename T2> + long double dif_rmse( const vector<T1> buf1,const vector<T2> buf2) + { + long double totalpower=0; + long double difpower=0; + size_t n = min( buf1.size(),buf2.size() ); + for (size_t k=0;k<n;++k) { + totalpower += (norm( buf1[k] ) + norm(buf2[k]) )/2.; + difpower += norm(buf1[k] - buf2[k]); + } + return sqrt(difpower/totalpower); + } + +template <class T> +void test_scalar(int nfft) +{ + typedef typename Eigen::FFT<T>::Complex Complex; + typedef typename Eigen::FFT<T>::Scalar Scalar; + + FFT<T> fft; + vector<Scalar> inbuf(nfft); + vector<Complex> outbuf; + 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 +} + template <class T> -void test_fft(int nfft) +void test_complex(int nfft) { typedef typename Eigen::FFT<T>::Complex Complex; FFT<T> fft; vector<Complex> inbuf(nfft); - vector<Complex> buf3(nfft); - vector<Complex> outbuf(nfft); + vector<Complex> outbuf; + vector<Complex> buf3; 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 + 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 + + fft.inv( buf3 , outbuf); + + VERIFY( dif_rmse(inbuf,buf3) < 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) )); + 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>(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) ); + CALL_SUBTEST( test_complex<float>(2*3*4*5) ); CALL_SUBTEST( test_complex<double>(2*3*4*5) ); CALL_SUBTEST( test_complex<long double>(2*3*4*5) ); + CALL_SUBTEST( test_complex<float>(2*3*4*5*7) ); CALL_SUBTEST( test_complex<double>(2*3*4*5*7) ); CALL_SUBTEST( test_complex<long double>(2*3*4*5*7) ); + + 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>(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) ); } |