[136] | 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 | }
|
---|