test_DuneFFT.cxx
Go to the documentation of this file.
1 // test_DuneFFT.cxx
2 //
3 // David Adams
4 // April 2019
5 //
6 // Test DuneFFT.
7 
9 #include <string>
10 #include <iostream>
11 #include <fstream>
12 #include <sstream>
13 #include <vector>
14 #include <iomanip>
15 #include <cmath>
16 
17 #undef NDEBUG
18 #include <cassert>
19 
20 using std::string;
21 using std::cout;
22 using std::endl;
23 using std::ostringstream;
24 using std::ofstream;
25 using std::vector;
26 using std::setw;
27 using std::fixed;
28 
29 using Index = unsigned int;
31 using DFT = DuneFFT::DFT;
32 
33 //**********************************************************************
34 
35 int test_DuneFFT(Index ignorm, Index itnorm, int loglev, Index len) {
36  const string myname = "test_DuneFFT: ";
37 #ifdef NDEBUG
38  cout << myname << "NDEBUG must be off." << endl;
39  abort();
40 #endif
41  string line = "-----------------------------";
42 
43  RealDftNormalization norm(ignorm, itnorm);
44  cout << line << endl;
45  cout << myname << " Global norm: " << ignorm << endl;
46  cout << myname << " Term norm: " << itnorm << endl;
47 
48  cout << myname << line << endl;
49  cout << myname << "Create data." << endl;
50  vector<float> sams = { 3.0, 6.0, 16.1, 28.6, 30.2, 27.7, 16.3, 9.6, 4.2, -1.0,
51  -2.3, -4.2, -9.2, -18.6, -21.9, -29.0, -24.3, -14.2, -5.0, -3.0};
52  if ( len > 0 ) sams.resize(len, 0.0);
53  Index nsam = sams.size();
54  assert( sams.size() == nsam );
55  cout << myname << "Sample length: " << nsam << endl;
56  float samsum = 0.0;
57  for ( float sam : sams ) samsum += sam;
58  cout << myname << "Sample mean: " << samsum/nsam << endl;
59  Index namp = (nsam+2)/2;
60  Index npha = (nsam+1)/2;
61  cout << myname << " # samples: " << nsam << endl;
62  cout << myname << "Expected amp size: " << namp << endl;
63  cout << myname << "Expected pha size: " << npha << endl;
64 
65  cout << myname << line << endl;
66  cout << myname << "Create empty DFT." << endl;
67  DFT dft(norm);
68  assert( dft.size() == 0 );
69 
70  cout << myname << line << endl;
71  cout << myname << "Transform forward." << endl;
72  assert( DuneFFT::fftForward(sams, dft, loglev) == 0 );
73  assert( sams.size() == nsam );
74  cout << myname << "DFT size: " << dft.size() << endl;
75  assert( dft.size() == nsam );
76  assert( dft.nCompact() == namp );
77  assert( dft.nAmplitude() == namp );
78  assert( dft.nPhase() == npha );
79  cout << "Frequency components:" << endl;
80  for ( Index ifrq=0; ifrq<namp; ++ifrq ) {
81  cout << myname << setw(4) << ifrq << ": " << dft.amplitude(ifrq);
82  if ( ifrq < npha ) cout << " @ " << dft.phase(ifrq);
83  cout << endl;
84  }
85 
86  cout << myname << line << endl;
87  cout << myname << "Check power." << endl;
88  float pwr1 = 0.0;
89  for ( float sam : sams ) pwr1 += sam*sam;
90  float pwr2 = dft.power();
91  float pwr3 = 0.0;
92  float pwr4 = 0.0;
93  for ( Index ifrq=0; ifrq<nsam; ++ifrq ) {
94  float amp = dft.amplitude(ifrq);
95  float xre = dft.real(ifrq);
96  float xim = dft.imag(ifrq);
97  if ( !dft.normalization().isPower() || ifrq < dft.nCompact() ) {
98  pwr3 += amp*amp;
99  pwr4 += xre*xre + xim*xim;
100  }
101  }
102  double pfac = 1.0;
103  if ( dft.normalization().isStandard() ) pfac = 1.0/nsam;
104  if ( dft.normalization().isBin() ) pfac = nsam;
105  pwr3 *= pfac;
106  pwr4 *= pfac;
107  cout << myname << "Tick power: " << pwr1 << endl;
108  cout << myname << " DFT power: " << pwr2 << endl;
109  cout << myname << " Amp power: " << pwr3 << endl;
110  cout << myname << " ReI power: " << pwr4 << endl;
111  assert( fabs(pwr2 - pwr1) < 1.e-5*pwr1 );
112 
113  cout << myname << line << endl;
114  cout << myname << "Call inverse." << endl;
115  FloatVector sams2, xres2, xims2;
116  int rstat = DuneFFT::fftInverse(dft, sams2, loglev);
117  if ( rstat != 0 ) {
118  cout << myname << "FFT invert returned " << rstat << endl;
119  assert( false );
120  }
121 
122  assert( sams2.size() == nsam );
123  for ( Index isam=0; isam<nsam; ++isam ) {
124  cout << setw(4) << isam << ":" << setw(10) << fixed << sams[isam]
125  << " ?= " << setw(10) << fixed << sams2[isam] << endl;
126  assert( fabs(sams2[isam] - sams[isam]) < 1.e-4 );
127  }
128 
129  cout << myname << line << endl;
130  cout << myname << "Done." << endl;
131  return 0;
132 }
133 
134 //**********************************************************************
135 
136 int main(int argc, char* argv[]) {
137  Index ignorm = 1;
138  Index itnorm = 1;
139  int loglev = 0;
140  Index len = 0;
141  if ( argc > 1 ) {
142  string sarg(argv[1]);
143  if ( sarg == "-h" ) {
144  cout << "Usage: " << argv[0] << " [ignorm [itnorm [loglev [LEN]]]]" << endl;
145  return 0;
146  }
147  ignorm = std::stoi(sarg);
148  }
149  if ( argc > 2 ) {
150  string sarg(argv[2]);
151  itnorm = std::stoi(sarg);
152  }
153  if ( argc > 3 ) {
154  string sarg(argv[3]);
155  loglev = std::stoi(sarg);
156  }
157  if ( argc > 4 ) {
158  string sarg(argv[4]);
159  len = std::stoi(sarg);
160  }
161  return test_DuneFFT(ignorm, itnorm, loglev, len);
162 }
163 
164 //**********************************************************************
std::vector< float > FloatVector
Definition: DuneFFT.h:29
std::string string
Definition: nybbler.cc:12
struct vector vector
unsigned int Index
CompactRealDftData< float > DFT
Definition: DuneFFT.h:28
int test_DuneFFT(Index ignorm, Index itnorm, int loglev, Index len)
F amplitude(Index ifrq) const override
F imag(Index ifrq) const override
Index size() const
Definition: RealDftData.h:35
Index nCompact() const override
F phase(Index ifrq) const override
const double e
Index nAmplitude() const
int main(int argc, char *argv[])
static int fftForward(Index ntick, const float *psam, DFT &dft, Index logLevel=0)
Definition: DuneFFT.cxx:23
size_t size
Definition: lodepng.cpp:55
F power() const
Definition: RealDftData.h:95
auto norm(Vector const &v)
Return norm of the specified vector.
Q_EXPORT QTSManip setw(int w)
Definition: qtextstream.h:331
static int fftInverse(const DFT &dft, FloatVector &sams, Index logLevel=0)
Definition: DuneFFT.cxx:79
F real(Index ifrq) const override
void line(double t, double *p, double &x, double &y, double &z)
const Norm & normalization() const override
Dft::FloatVector FloatVector
QTextStream & endl(QTextStream &s)