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