1//  To use the simple FFT implementation
2//  g++ -o demofft -I.. -Wall -O3 FFT.cpp
3
4//  To use the FFTW implementation
5//  g++ -o demofft -I.. -DUSE_FFTW -Wall -O3 FFT.cpp -lfftw3 -lfftw3f -lfftw3l
6
7#ifdef USE_FFTW
8#include <fftw3.h>
9#endif
10
11#include <vector>
12#include <complex>
13#include <algorithm>
14#include <iterator>
15#include <iostream>
16#include <Eigen/Core>
17#include <unsupported/Eigen/FFT>
18
19using namespace std;
20using namespace Eigen;
21
22template <typename T>
23T mag2(T a)
24{
25    return a*a;
26}
27template <typename T>
28T mag2(std::complex<T> a)
29{
30    return norm(a);
31}
32
33template <typename T>
34T mag2(const std::vector<T> & vec)
35{
36    T out=0;
37    for (size_t k=0;k<vec.size();++k)
38        out += mag2(vec[k]);
39    return out;
40}
41
42template <typename T>
43T mag2(const std::vector<std::complex<T> > & vec)
44{
45    T out=0;
46    for (size_t k=0;k<vec.size();++k)
47        out += mag2(vec[k]);
48    return out;
49}
50
51template <typename T>
52vector<T> operator-(const vector<T> & a,const vector<T> & b )
53{
54    vector<T> c(a);
55    for (size_t k=0;k<b.size();++k)
56        c[k] -= b[k];
57    return c;
58}
59
60template <typename T>
61void RandomFill(std::vector<T> & vec)
62{
63    for (size_t k=0;k<vec.size();++k)
64        vec[k] = T( rand() )/T(RAND_MAX) - .5;
65}
66
67template <typename T>
68void RandomFill(std::vector<std::complex<T> > & vec)
69{
70    for (size_t k=0;k<vec.size();++k)
71        vec[k] = std::complex<T> ( T( rand() )/T(RAND_MAX) - .5, T( rand() )/T(RAND_MAX) - .5);
72}
73
74template <typename T_time,typename T_freq>
75void fwd_inv(size_t nfft)
76{
77    typedef typename NumTraits<T_freq>::Real Scalar;
78    vector<T_time> timebuf(nfft);
79    RandomFill(timebuf);
80
81    vector<T_freq> freqbuf;
82    static FFT<Scalar> fft;
83    fft.fwd(freqbuf,timebuf);
84
85    vector<T_time> timebuf2;
86    fft.inv(timebuf2,freqbuf);
87
88    long double rmse = mag2(timebuf - timebuf2) / mag2(timebuf);
89    cout << "roundtrip rmse: " << rmse << endl;
90}
91
92template <typename T_scalar>
93void two_demos(int nfft)
94{
95    cout << "     scalar ";
96    fwd_inv<T_scalar,std::complex<T_scalar> >(nfft);
97    cout << "    complex ";
98    fwd_inv<std::complex<T_scalar>,std::complex<T_scalar> >(nfft);
99}
100
101void demo_all_types(int nfft)
102{
103    cout << "nfft=" << nfft << endl;
104    cout << "   float" << endl;
105    two_demos<float>(nfft);
106    cout << "   double" << endl;
107    two_demos<double>(nfft);
108    cout << "   long double" << endl;
109    two_demos<long double>(nfft);
110}
111
112int main()
113{
114    demo_all_types( 2*3*4*5*7 );
115    demo_all_types( 2*9*16*25 );
116    demo_all_types( 1024 );
117    return 0;
118}
119