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 | }
|
---|