• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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 
19 using namespace std;
20 using namespace Eigen;
21 
22 template <typename T>
mag2(T a)23 T mag2(T a)
24 {
25     return a*a;
26 }
27 template <typename T>
mag2(std::complex<T> a)28 T mag2(std::complex<T> a)
29 {
30     return norm(a);
31 }
32 
33 template <typename T>
mag2(const std::vector<T> & vec)34 T 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 
42 template <typename T>
mag2(const std::vector<std::complex<T>> & vec)43 T 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 
51 template <typename T>
operator -(const vector<T> & a,const vector<T> & b)52 vector<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 
60 template <typename T>
RandomFill(std::vector<T> & vec)61 void 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 
67 template <typename T>
RandomFill(std::vector<std::complex<T>> & vec)68 void 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 
74 template <typename T_time,typename T_freq>
fwd_inv(size_t nfft)75 void 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 
92 template <typename T_scalar>
two_demos(int nfft)93 void 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 
demo_all_types(int nfft)101 void 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 
main()112 int 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