| 24 | |
| 25 | template <class T> |
| 26 | void dotest(int nfft) |
| 27 | { |
| 28 | typedef kissfft<T> FFT; |
| 29 | typedef std::complex<T> cpx_type; |
| 30 | |
| 31 | cout << "type:" << typeid(T).name() << " nfft:" << nfft; |
| 32 | |
| 33 | FFT fft(nfft,false); |
| 34 | |
| 35 | vector<cpx_type> inbuf(nfft); |
| 36 | vector<cpx_type> outbuf(nfft); |
| 37 | for (int k=0;k<nfft;++k) |
| 38 | inbuf[k]= cpx_type( |
| 39 | (T)(rand()/(double)RAND_MAX - .5), |
| 40 | (T)(rand()/(double)RAND_MAX - .5) ); |
| 41 | fft.transform( &inbuf[0] , &outbuf[0] ); |
| 42 | |
| 43 | long double totalpower=0; |
| 44 | long double difpower=0; |
| 45 | |
| 46 | // Create long double constant for pi because M_PIl is not defined by |
| 47 | // all toolchains. |
| 48 | const long double pi = std::acosl(-1); |
| 49 | |
| 50 | for (int k0=0;k0<nfft;++k0) { |
| 51 | complex<long double> acc = 0; |
| 52 | long double phinc = 2*k0* pi / nfft; |
| 53 | for (int k1=0;k1<nfft;++k1) { |
| 54 | complex<long double> x(inbuf[k1].real(),inbuf[k1].imag()); |
| 55 | acc += x * exp( complex<long double>(0,-k1*phinc) ); |
| 56 | } |
| 57 | totalpower += norm(acc); |
| 58 | complex<long double> x(outbuf[k0].real(),outbuf[k0].imag()); |
| 59 | complex<long double> dif = acc - x; |
| 60 | difpower += norm(dif); |
| 61 | } |
| 62 | cout << " RMSE:" << sqrt(difpower/totalpower) << "\t"; |
| 63 | |
| 64 | double t0 = curtime(); |
| 65 | int nits=20e6/nfft; |
| 66 | for (int k=0;k<nits;++k) { |
| 67 | fft.transform( &inbuf[0] , &outbuf[0] ); |
| 68 | } |
| 69 | double t1 = curtime(); |
| 70 | cout << " MSPS:" << ( (nits*nfft)*1e-6/ (t1-t0) ) << endl; |
| 71 | } |
| 72 | |
| 73 | int main(int argc,char ** argv) |
| 74 | { |