blob: a62f6e0aa2caf3abce88147a7195f5a9ff715aeb (
plain) (
blame)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
|
/*
* Copyright (c) 2003-2010, Mark Borgerding. All rights reserved.
* This file is part of KISS FFT - https://github.com/mborgerding/kissfft
*
* SPDX-License-Identifier: BSD-3-Clause
* See COPYING file for more information.
*/
#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;
}
|
