lasp/src/lasp/dsp/lasp_fft.cpp

// lasp_fft.c
//
// Author: J.A. de Jong - ASCEE
//
// Description:
// FFt implementation
//
//////////////////////////////////////////////////////////////////////
#include <memory>
#include <cassert>
#define DEBUGTRACE_ENABLED
#include "lasp_fft.h"
#include "debugtrace.hpp"
#include "lasp_config.h"
using rte = std::runtime_error;

#if LASP_FFT_BACKEND == Armadillo
#include "fftpack.h"
class Fft_impl {
  public:
    us nfft;
    Fft_impl(const us nfft) : nfft(nfft) {
      throw runtime_error(
          "This code does not output correct results, as it computes the full "
          "FFT. It needs to be reworked to single-sided spectra.");

      DEBUGTRACE_ENTER;
    }

    vd ifft(const vc &f) { return arma::ifft(f); }
    vc fft(const vd &time) { return arma::fft(time); }
};
#elif LASP_FFT_BACKEND == FFTW
#include <fftw3.h>

/**
 * @brief FFTW implementation
 */
class Fft_impl {
  public:
    us nfft;
    fftw_plan forward_plan = nullptr;
    fftw_plan reverse_plan = nullptr;
    fftw_complex *frequencyDomain = nullptr;
    d *timeDomain = nullptr;

    Fft_impl(const us nfft) : nfft(nfft) {

      timeDomain = (d *)fftw_malloc(sizeof(d) * nfft);

      frequencyDomain =
        (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * (nfft / 2 + 1));

      forward_plan =
        fftw_plan_dft_r2c_1d(nfft, timeDomain, frequencyDomain, FFTW_MEASURE);

      reverse_plan =
        fftw_plan_dft_c2r_1d(nfft, frequencyDomain, timeDomain, FFTW_MEASURE);
      if(!forward_plan || !reverse_plan || !timeDomain || !frequencyDomain) {
        throw rte("Error allocating FFT");
      }

    };

    ~Fft_impl() {
      fftw_destroy_plan(forward_plan);
      fftw_destroy_plan(reverse_plan);
      fftw_free(frequencyDomain);
      fftw_free(timeDomain);
    }

    vc fft(const vd &time) {
      if(time.n_rows != nfft) {
        throw rte("Invalid size of input vector, should be equal to nfft"); 
      }

      vc res(nfft / 2 + 1);
      memcpy(timeDomain, time.memptr(), sizeof(d) * nfft);

      fftw_execute(forward_plan);

      memcpy(reinterpret_cast<void *>(res.memptr()), frequencyDomain,
          sizeof(c) * (nfft / 2 + 1));

      return res;
    }
    vd ifft(const vc &f) {
      DEBUGTRACE_ENTER;
      if(f.n_rows != nfft/2+1) {
        throw rte("Invalid size of input vector, should be equal to nfft/2+1"); 
      }
      memcpy(frequencyDomain,
          reinterpret_cast<void *>(const_cast<c *>(f.memptr())),
          (nfft / 2 + 1) * sizeof(c));

      fftw_execute(reverse_plan);

      vd res(nfft);
      memcpy(res.memptr(), timeDomain, nfft * sizeof(d));

      /* Scale by dividing by nfft. Checked with numpy implementation
       * that this indeed needs to be done for FFTW. */
      res *= 1.0 / nfft;

      return res;
    }
};
#else
#error                                                                         \
  "Cannot compile lasp_ffc.c, no FFT backend specified. Should either be FFTPack, or FFTW"
#endif

Fft::Fft(const us nfft) {
  DEBUGTRACE_ENTER;
  if (nfft == 0) {
    throw rte("Invalid nfft: 0");
  }
  if (nfft >= LASP_MAX_NFFT) {
    throw rte("Invalid nfft, should be smaller than: " +
        std::to_string(LASP_MAX_NFFT));
  }
  _impl = std::make_unique<Fft_impl>(nfft);
}
Fft::~Fft() { }

us Fft::nfft() const { assert(_impl); return _impl->nfft; }

vc Fft::fft(const vd &timedata) { assert(_impl); return _impl->fft(timedata); }

cmat Fft::fft(const dmat &freqdata) {
  DEBUGTRACE_ENTER;
  assert(_impl); 
  cmat res(_impl->nfft/2+1, freqdata.n_cols);
#pragma omp parallel for
  for (us colno = 0; colno < freqdata.n_cols; colno++) {
    res.col(colno) = _impl->fft(freqdata.col(colno));
  }
  return res;
}

vd Fft::ifft(const vc &freqdata) {
  DEBUGTRACE_ENTER;
  assert(_impl); 
  return _impl->ifft(freqdata);
}
dmat Fft::ifft(const cmat &freqdata) {
  dmat res(_impl->nfft, freqdata.n_cols);
#pragma omp parallel for
  for (us colno = 0; colno < freqdata.n_cols; colno++) {
    res.col(colno) = _impl->ifft(freqdata.col(colno));
  }
  return res;
}

void load_fft_wisdom(const char *wisdom) {
#if LASP_FFT_BACKEND == Armadillo
#elif LASP_FFT_BACKEND == FFTW
  if (wisdom) {
    int rv = fftw_import_wisdom_from_string(wisdom);
    if (rv != 1) {
      fprintf(stderr, "Error loading FFTW wisdom");
    }
  }
#endif
}

char *store_fft_wisdom() {
#if LASP_FFT_BACKEND == Armadillo
  return NULL;
#elif LASP_FFT_BACKEND == FFTW
  return fftw_export_wisdom_to_string();
#endif
}