Fftw.cpp

Go to the documentation of this file.

 
#include <cmath>  // for cos
#include <cstdio>
#include <map>     // for map
#include <string>  // for string
 
#include <boost/format.hpp>  // for format
 
#include <fftw3.h>
 
#include "ElementsKernel/MathConstants.h"   // for pi
#include "ElementsKernel/ProgramHeaders.h"  // for including all Program/related headers
#include "ElementsKernel/Unused.h"          // for ELEMENTS_UNUSED
 
using std::map;
using std::string;
 
constexpr std::size_t N = 32;
 
namespace Elements {
namespace Examples {
 
class Fftw : public Program {
 
public:
  ExitCode mainMethod(ELEMENTS_UNUSED map<string, VariableValue>& args) override {
 
    auto log = Logging::getLogger("FftwExample");
 
    fftw_complex in[N]; /* double [2] */
    fftw_complex out[N];
    fftw_complex in2[N];
    fftw_plan    p;
    fftw_plan    q;
 
    using std::cos;
 
    /* prepare a cosine wave */
    for (size_t i = 0; i < N; i++) {
      in[i][0] = cos(3.0 * 2.0 * Units::pi * static_cast<double>(i) / static_cast<double>(N));
      in[i][1] = 0;
    }
 
    /* forward Fourier transform, save the result in 'out' */
    p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(p);
    for (size_t i = 0; i < N; i++) {
      log.info() << boost::format("freq: %3d %+9.5f %+9.5f I") % i % out[i][0] % out[i][1];
    }
    fftw_destroy_plan(p);
 
    /* backward Fourier transform, save the result in 'in2' */
    printf("\nInverse transform:\n");
    q = fftw_plan_dft_1d(N, out, in2, FFTW_BACKWARD, FFTW_ESTIMATE);
    fftw_execute(q);
    /* normalize */
    for (size_t i = 0; i < N; i++) {
      in2[i][0] *= 1. / N;
      in2[i][1] *= 1. / N;
    }
    for (size_t i = 0; i < N; i++) {
      log.info() << boost::format("recover: %3d %+9.5f %+9.5f I vs. %+9.5f %+9.5f I") % i % in[i][0] % in[i][1] %
                        in2[i][0] % in2[i][1];
    }
    fftw_destroy_plan(q);
 
    fftw_cleanup();
 
    log.info() << "This is the end of the test";
 
    return ExitCode::OK;
  }
};
 
}  // namespace Examples
}  // namespace Elements
 
MAIN_FOR(Elements::Examples::Fftw)