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