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)23T mag2(T a) 24 { 25 return a*a; 26 } 27 template <typename T> mag2(std::complex<T> a)28T mag2(std::complex<T> a) 29 { 30 return norm(a); 31 } 32 33 template <typename T> mag2(const std::vector<T> & vec)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 42 template <typename T> mag2(const std::vector<std::complex<T>> & vec)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 51 template <typename T> operator -(const vector<T> & a,const vector<T> & b)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 60 template <typename T> RandomFill(std::vector<T> & vec)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 67 template <typename T> RandomFill(std::vector<std::complex<T>> & vec)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 74 template <typename T_time,typename T_freq> fwd_inv(size_t nfft)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 92 template <typename T_scalar> two_demos(int nfft)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 demo_all_types(int nfft)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 main()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