doc/examples/FFT.cpp
Go to the documentation of this file.
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>
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) - T(.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) - T(.5), T( rand() )/T(RAND_MAX) - T(.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  T_time 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 }
void two_demos(int nfft)
SCALAR Scalar
Definition: bench_gemm.cpp:33
void fwd_inv(size_t nfft)
Scalar * b
Definition: benchVecAdd.cpp:17
Scalar Scalar * c
Definition: benchVecAdd.cpp:17
void demo_all_types(int nfft)
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
Definition: Half.h:150
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:150
vector< T > operator-(const vector< T > &a, const vector< T > &b)
int main()
Array33i a
void RandomFill(std::vector< T > &vec)
Eigen::Triplet< double > T
A small structure to hold a non zero as a triplet (i,j,value).
Definition: SparseUtil.h:154
T mag2(T a)


gtsam
Author(s):
autogenerated on Sat May 8 2021 02:42:04