c@409: #include "kissfft.hh"
c@409: #include <iostream>
c@409: #include <cstdlib>
c@409: #include <typeinfo>
c@409: 
c@409: #include <sys/time.h>
c@409: static inline
c@409: double curtime(void)
c@409: {
c@409:     struct timeval tv;
c@409:     gettimeofday(&tv, NULL);
c@409:     return (double)tv.tv_sec + (double)tv.tv_usec*.000001;
c@409: }
c@409: 
c@409: using namespace std;
c@409: 
c@409: template <class T>
c@409: void dotest(int nfft)
c@409: {
c@409:     typedef kissfft<T> FFT;
c@409:     typedef std::complex<T> cpx_type;
c@409: 
c@409:     cout << "type:" << typeid(T).name() << " nfft:" << nfft;
c@409: 
c@409:     FFT fft(nfft,false);
c@409: 
c@409:     vector<cpx_type> inbuf(nfft);
c@409:     vector<cpx_type> outbuf(nfft);
c@409:     for (int k=0;k<nfft;++k)
c@409:         inbuf[k]= cpx_type( 
c@409:                 (T)(rand()/(double)RAND_MAX - .5),
c@409:                 (T)(rand()/(double)RAND_MAX - .5) );
c@409:     fft.transform( &inbuf[0] , &outbuf[0] );
c@409: 
c@409:     long double totalpower=0;
c@409:     long double difpower=0;
c@409:     for (int k0=0;k0<nfft;++k0) {
c@409:         complex<long double> acc = 0;
c@409:         long double phinc = 2*k0* M_PIl / nfft;
c@409:         for (int k1=0;k1<nfft;++k1) {
c@409:             complex<long double> x(inbuf[k1].real(),inbuf[k1].imag()); 
c@409:             acc += x * exp( complex<long double>(0,-k1*phinc) );
c@409:         }
c@409:         totalpower += norm(acc);
c@409:         complex<long double> x(outbuf[k0].real(),outbuf[k0].imag()); 
c@409:         complex<long double> dif = acc - x;
c@409:         difpower += norm(dif);
c@409:     }
c@409:     cout << " RMSE:" << sqrt(difpower/totalpower) << "\t";
c@409: 
c@409:     double t0 = curtime();
c@409:     int nits=20e6/nfft;
c@409:     for (int k=0;k<nits;++k) {
c@409:         fft.transform( &inbuf[0] , &outbuf[0] );
c@409:     }
c@409:     double t1 = curtime();
c@409:     cout << " MSPS:" << ( (nits*nfft)*1e-6/ (t1-t0) ) << endl;
c@409: }
c@409: 
c@409: int main(int argc,char ** argv)
c@409: {
c@409:     if (argc>1) {
c@409:         for (int k=1;k<argc;++k) {
c@409:             int nfft = atoi(argv[k]);
c@409:             dotest<float>(nfft); dotest<double>(nfft); dotest<long double>(nfft);
c@409:         }
c@409:     }else{
c@409:         dotest<float>(32); dotest<double>(32); dotest<long double>(32);
c@409:         dotest<float>(1024); dotest<double>(1024); dotest<long double>(1024);
c@409:         dotest<float>(840); dotest<double>(840); dotest<long double>(840);
c@409:     }
c@409:     return 0;
c@409: }