FFT.cpp
Go to the documentation of this file.
00001 //  To use the simple FFT implementation
00002 //  g++ -o demofft -I.. -Wall -O3 FFT.cpp 
00003 
00004 //  To use the FFTW implementation
00005 //  g++ -o demofft -I.. -DUSE_FFTW -Wall -O3 FFT.cpp -lfftw3 -lfftw3f -lfftw3l
00006 
00007 #ifdef USE_FFTW
00008 #include <fftw3.h>
00009 #endif
00010 
00011 #include <vector>
00012 #include <complex>
00013 #include <algorithm>
00014 #include <iterator>
00015 #include <iostream>
00016 #include <Eigen/Core>
00017 #include <unsupported/Eigen/FFT>
00018 
00019 using namespace std;
00020 using namespace Eigen;
00021 
00022 template <typename T>
00023 T mag2(T a)
00024 {
00025     return a*a;
00026 }
00027 template <typename T>
00028 T mag2(std::complex<T> a)
00029 {
00030     return norm(a);
00031 }
00032 
00033 template <typename T>
00034 T mag2(const std::vector<T> & vec)
00035 {
00036     T out=0;
00037     for (size_t k=0;k<vec.size();++k)
00038         out += mag2(vec[k]);
00039     return out;
00040 }
00041 
00042 template <typename T>
00043 T mag2(const std::vector<std::complex<T> > & vec)
00044 {
00045     T out=0;
00046     for (size_t k=0;k<vec.size();++k)
00047         out += mag2(vec[k]);
00048     return out;
00049 }
00050 
00051 template <typename T>
00052 vector<T> operator-(const vector<T> & a,const vector<T> & b )
00053 {
00054     vector<T> c(a);
00055     for (size_t k=0;k<b.size();++k) 
00056         c[k] -= b[k];
00057     return c;
00058 }
00059 
00060 template <typename T>
00061 void RandomFill(std::vector<T> & vec)
00062 {
00063     for (size_t k=0;k<vec.size();++k)
00064         vec[k] = T( rand() )/T(RAND_MAX) - .5;
00065 }
00066 
00067 template <typename T>
00068 void RandomFill(std::vector<std::complex<T> > & vec)
00069 {
00070     for (size_t k=0;k<vec.size();++k)
00071         vec[k] = std::complex<T> ( T( rand() )/T(RAND_MAX) - .5, T( rand() )/T(RAND_MAX) - .5);
00072 }
00073 
00074 template <typename T_time,typename T_freq>
00075 void fwd_inv(size_t nfft)
00076 {
00077     typedef typename NumTraits<T_freq>::Real Scalar;
00078     vector<T_time> timebuf(nfft);
00079     RandomFill(timebuf);
00080 
00081     vector<T_freq> freqbuf;
00082     static FFT<Scalar> fft;
00083     fft.fwd(freqbuf,timebuf);
00084 
00085     vector<T_time> timebuf2;
00086     fft.inv(timebuf2,freqbuf);
00087 
00088     long double rmse = mag2(timebuf - timebuf2) / mag2(timebuf);
00089     cout << "roundtrip rmse: " << rmse << endl;
00090 }
00091 
00092 template <typename T_scalar>
00093 void two_demos(int nfft)
00094 {
00095     cout << "     scalar ";
00096     fwd_inv<T_scalar,std::complex<T_scalar> >(nfft);
00097     cout << "    complex ";
00098     fwd_inv<std::complex<T_scalar>,std::complex<T_scalar> >(nfft);
00099 }
00100 
00101 void demo_all_types(int nfft)
00102 {
00103     cout << "nfft=" << nfft << endl;
00104     cout << "   float" << endl;
00105     two_demos<float>(nfft);
00106     cout << "   double" << endl;
00107     two_demos<double>(nfft);
00108     cout << "   long double" << endl;
00109     two_demos<long double>(nfft);
00110 }
00111 
00112 int main()
00113 {
00114     demo_all_types( 2*3*4*5*7 );
00115     demo_all_types( 2*9*16*25 );
00116     demo_all_types( 1024 );
00117     return 0;
00118 }


re_vision
Author(s): Dorian Galvez-Lopez
autogenerated on Sun Jan 5 2014 11:31:08