JFFT.cc File Reference

Auxiliary program to run fast Fourier transformation (FFT). More...

#include <string>
#include <iostream>
#include <fstream>
#include <iomanip>
#include <vector>
#include <fftw3.h>
#include <cmath>
#include "JLang/gzstream.h"
#include "Jeep/JeepToolkit.hh"
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"

Go to the source code of this file.

Functions
int	main (int argc, char **argv)

Detailed Description

Auxiliary program to run fast Fourier transformation (FFT).

For more information, see FFTW.

Author

mdejong

Definition in file JFFT.cc.

Function Documentation

int main	(	int	argc,
		char **	argv
	)

Definition at line 23 of file JFFT.cc.

{
  using namespace std;
  using namespace JPP;

  string  inputFile;
  string  outputFile;
  double  binWidth;
  int     precision;
  int     debug;

  try { 

    JParser<> zap("Auxiliary program to run fast Fourier transformation (FFT).");
    
    zap['f'] = make_field(inputFile,  "input file (containing 1D array of values)");
    zap['o'] = make_field(outputFile, "output file (containing 2D array of values)")   = "fft.txt";
    zap['B'] = make_field(binWidth,   "bin width of input array (if zero, use index)") = 0.0;
    zap['F'] = make_field(precision,  "number of decimals if output format")           = 9;
    zap['d'] = make_field(debug)           = 2;
    
    zap(argc, argv);
  }
  catch(const exception& error) {
    FATAL(error.what() << endl);
  }

  // input
  
  vector<double> data;

  istream* in = open<istream>(inputFile);

  if (in == NULL) {
    FATAL("Invalid file " << inputFile << endl);
  }

  for (double x; *in >> x; ) {
    data.push_back(x);
  }

  close(in);

  const size_t N  = data.size(); 
 
  if (N == 0) {
    FATAL("Number of points: " << N << endl);
  }
  
  // FFT

  const size_t  NS     = N/2 + 1;
  double*       input  = (double*) malloc(sizeof(double) *  N);
  fftw_complex* output = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * NS);
  fftw_plan     plan   = fftw_plan_dft_r2c_1d(N, input, output, FFTW_MEASURE);

  for (size_t i = 0; i != N; ++i) {
    input[i] = data[i];
  }

  // execute plan

  fftw_execute(plan);
        
  // destroy plan

  fftw_destroy_plan(plan);

  // output

  const double xmin  = (binWidth != 0.0 ? 1.0 / (N * binWidth) : 0.0);
  const double xmax  = (binWidth != 0.0 ? 1.0 / (binWidth)     : 0.0);

  const int    width = precision + 7;

  ostream* out = open<ostream>(outputFile.c_str());

  for (size_t i = 0; i != NS; ++i) {

    const double real  = output[i][0];
    const double imag  = output[i][1];
    const double power = sqrt(real*real + imag*imag);

    if (binWidth != 0.0)
      *out << SCIENTIFIC(width,precision) << xmin + (double) i * (xmax - xmin) / (double) (N - 1);
    else
      *out << i;

    *out << ' ';
    *out << SCIENTIFIC(width,precision) << power << endl;
  }

  fftw_free(output);
  free(input);

  close(out);  
}