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