Home | History | Annotate | Download | only in examples 
      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>
     23 T mag2(T a)
     24 {
     25     return a*a;
     26 }
     27 template <typename T>
     28 T mag2(std::complex<T> a)
     29 {
     30     return norm(a);
     31 }
     32 
     33 template <typename T>
     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>
     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>
     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>
     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>
     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>
     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>
     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 
    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 
    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