Jpp  19.1.0-rc.1
the software that should make you happy
JFFT.cc
Go to the documentation of this file.
1 #include <string>
2 #include <iostream>
3 #include <fstream>
4 #include <iomanip>
5 #include <vector>
6 #include <fftw3.h>
7 #include <cmath>
8 
9 #include "JLang/gzstream.h"
10 
11 #include "Jeep/JeepToolkit.hh"
12 #include "Jeep/JParser.hh"
13 #include "Jeep/JMessage.hh"
14 
15 
16 /**
17  * \file
18  * Auxiliary program to run fast Fourier transformation (FFT).\n
19  * For more information, see <a href="https://www.fftw.org/">FFTW</a>.
20  *
21  * \author mdejong
22  */
23 int main(int argc, char **argv)
24 {
25  using namespace std;
26  using namespace JPP;
27 
28  string inputFile;
29  string outputFile;
30  double binWidth;
31  int precision;
32  int debug;
33 
34  try {
35 
36  JParser<> zap("Auxiliary program to run fast Fourier transformation (FFT).");
37 
38  zap['f'] = make_field(inputFile, "input file (containing 1D array of values)");
39  zap['o'] = make_field(outputFile, "output file (containing 2D array of values)") = "fft.txt";
40  zap['B'] = make_field(binWidth, "bin width of input array (if zero, use index)") = 0.0;
41  zap['F'] = make_field(precision, "number of decimals if output format") = 9;
42  zap['d'] = make_field(debug) = 2;
43 
44  zap(argc, argv);
45  }
46  catch(const exception& error) {
47  FATAL(error.what() << endl);
48  }
49 
50  // input
51 
53 
54  istream* in = open<istream>(inputFile);
55 
56  if (in == NULL) {
57  FATAL("Invalid file " << inputFile << endl);
58  }
59 
60  for (double x; *in >> x; ) {
61  data.push_back(x);
62  }
63 
64  close(in);
65 
66  const size_t N = data.size();
67 
68  if (N == 0) {
69  FATAL("Number of points: " << N << endl);
70  }
71 
72  // FFT
73 
74  const size_t NS = N/2 + 1;
75  double* input = (double*) malloc(sizeof(double) * N);
76  fftw_complex* output = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * NS);
77  fftw_plan plan = fftw_plan_dft_r2c_1d(N, input, output, FFTW_MEASURE);
78 
79  for (size_t i = 0; i != N; ++i) {
80  input[i] = data[i];
81  }
82 
83  // execute plan
84 
85  fftw_execute(plan);
86 
87  // destroy plan
88 
89  fftw_destroy_plan(plan);
90 
91  // output
92 
93  const double xmin = (binWidth != 0.0 ? 1.0 / (N * binWidth) : 0.0);
94  const double xmax = (binWidth != 0.0 ? 1.0 / (binWidth) : 0.0);
95 
96  const int width = precision + 7;
97 
98  ostream* out = open<ostream>(outputFile.c_str());
99 
100  for (size_t i = 0; i != NS; ++i) {
101 
102  const double real = output[i][0];
103  const double imag = output[i][1];
104  const double power = sqrt(real*real + imag*imag);
105 
106  if (binWidth != 0.0)
107  *out << SCIENTIFIC(width,precision) << xmin + (double) i * (xmax - xmin) / (double) (N - 1);
108  else
109  *out << i;
110 
111  *out << ' ';
112  *out << SCIENTIFIC(width,precision) << power << endl;
113  }
114 
115  fftw_free(output);
116  free(input);
117 
118  close(out);
119 }
120 
string outputFile
int main(int argc, char **argv)
Definition: JFFT.cc:23
General purpose messaging.
#define FATAL(A)
Definition: JMessage.hh:67
int debug
debug level
Definition: JSirene.cc:69
Utility class to parse command line options.
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition: JParser.hh:2158
Auxiliary methods for handling file names, type names and environment.
Utility class to parse command line options.
Definition: JParser.hh:1714
const double xmax
Definition: JQuadrature.cc:24
const double xmin
Definition: JQuadrature.cc:23
void close(std::istream *pf)
Close file.
Definition: JeepToolkit.hh:391
This name space includes all other name spaces (except KM3NETDAQ, KM3NET and ANTARES).
Definition: JSTDTypes.hh:14
Auxiliary data structure for floating point format specification.
Definition: JManip.hh:488