comparison ext/kissfft/test/testcpp.cc @ 184:76ec2365b250

Bring in kissfft into this repo (formerly a subrepo, but the remote is not responding)

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