annotate constant-q-cpp/src/ext/kissfft/test/testcpp.cc @ 372:af71cbdab621 tip

Update bqvec code
author Chris Cannam
date Tue, 19 Nov 2019 10:13:32 +0000
parents 5d0a2ebb4d17
children
rev   line source
Chris@366 1 #include "kissfft.hh"
Chris@366 2 #include <iostream>
Chris@366 3 #include <cstdlib>
Chris@366 4 #include <typeinfo>
Chris@366 5
Chris@366 6 #include <sys/time.h>
Chris@366 7 static inline
Chris@366 8 double curtime(void)
Chris@366 9 {
Chris@366 10 struct timeval tv;
Chris@366 11 gettimeofday(&tv, NULL);
Chris@366 12 return (double)tv.tv_sec + (double)tv.tv_usec*.000001;
Chris@366 13 }
Chris@366 14
Chris@366 15 using namespace std;
Chris@366 16
Chris@366 17 template <class T>
Chris@366 18 void dotest(int nfft)
Chris@366 19 {
Chris@366 20 typedef kissfft<T> FFT;
Chris@366 21 typedef std::complex<T> cpx_type;
Chris@366 22
Chris@366 23 cout << "type:" << typeid(T).name() << " nfft:" << nfft;
Chris@366 24
Chris@366 25 FFT fft(nfft,false);
Chris@366 26
Chris@366 27 vector<cpx_type> inbuf(nfft);
Chris@366 28 vector<cpx_type> outbuf(nfft);
Chris@366 29 for (int k=0;k<nfft;++k)
Chris@366 30 inbuf[k]= cpx_type(
Chris@366 31 (T)(rand()/(double)RAND_MAX - .5),
Chris@366 32 (T)(rand()/(double)RAND_MAX - .5) );
Chris@366 33 fft.transform( &inbuf[0] , &outbuf[0] );
Chris@366 34
Chris@366 35 long double totalpower=0;
Chris@366 36 long double difpower=0;
Chris@366 37 for (int k0=0;k0<nfft;++k0) {
Chris@366 38 complex<long double> acc = 0;
Chris@366 39 long double phinc = 2*k0* M_PIl / nfft;
Chris@366 40 for (int k1=0;k1<nfft;++k1) {
Chris@366 41 complex<long double> x(inbuf[k1].real(),inbuf[k1].imag());
Chris@366 42 acc += x * exp( complex<long double>(0,-k1*phinc) );
Chris@366 43 }
Chris@366 44 totalpower += norm(acc);
Chris@366 45 complex<long double> x(outbuf[k0].real(),outbuf[k0].imag());
Chris@366 46 complex<long double> dif = acc - x;
Chris@366 47 difpower += norm(dif);
Chris@366 48 }
Chris@366 49 cout << " RMSE:" << sqrt(difpower/totalpower) << "\t";
Chris@366 50
Chris@366 51 double t0 = curtime();
Chris@366 52 int nits=20e6/nfft;
Chris@366 53 for (int k=0;k<nits;++k) {
Chris@366 54 fft.transform( &inbuf[0] , &outbuf[0] );
Chris@366 55 }
Chris@366 56 double t1 = curtime();
Chris@366 57 cout << " MSPS:" << ( (nits*nfft)*1e-6/ (t1-t0) ) << endl;
Chris@366 58 }
Chris@366 59
Chris@366 60 int main(int argc,char ** argv)
Chris@366 61 {
Chris@366 62 if (argc>1) {
Chris@366 63 for (int k=1;k<argc;++k) {
Chris@366 64 int nfft = atoi(argv[k]);
Chris@366 65 dotest<float>(nfft); dotest<double>(nfft); dotest<long double>(nfft);
Chris@366 66 }
Chris@366 67 }else{
Chris@366 68 dotest<float>(32); dotest<double>(32); dotest<long double>(32);
Chris@366 69 dotest<float>(1024); dotest<double>(1024); dotest<long double>(1024);
Chris@366 70 dotest<float>(840); dotest<double>(840); dotest<long double>(840);
Chris@366 71 }
Chris@366 72 return 0;
Chris@366 73 }