impl.hpp

Go to the documentation of this file.

/// @file spectral/backends/fftw/impl.hpp
/// @brief FFTW3 backend wrappers.
/// Only included by src/spectral/fft.cpp.
#pragma once
 
#ifdef NUMERICS_HAS_FFTW
    #include "spectral/fft.hpp"
    #include <fftw3.h>
    #include <stdexcept>
 
namespace backends {
namespace fftw {
 
// std::complex<double> is layout-compatible with fftw_complex (double[2])
inline fftw_complex* fc(num::CVector& v) {
    return reinterpret_cast<fftw_complex*>(v.data());
}
inline const fftw_complex* fc(const num::CVector& v) {
    return reinterpret_cast<const fftw_complex*>(v.data());
}
 
inline void fft(const num::CVector& in, num::CVector& out) {
    fftw_plan p = fftw_plan_dft_1d(static_cast<int>(in.size()),
                                   const_cast<fftw_complex*>(fc(in)),
                                   fc(out),
                                   FFTW_FORWARD,
                                   FFTW_ESTIMATE);
    fftw_execute(p);
    fftw_destroy_plan(p);
}
 
inline void ifft(const num::CVector& in, num::CVector& out) {
    fftw_plan p = fftw_plan_dft_1d(static_cast<int>(in.size()),
                                   const_cast<fftw_complex*>(fc(in)),
                                   fc(out),
                                   FFTW_BACKWARD,
                                   FFTW_ESTIMATE);
    fftw_execute(p);
    fftw_destroy_plan(p);
}
 
inline void rfft(const num::Vector& in, num::CVector& out) {
    fftw_plan p = fftw_plan_dft_r2c_1d(static_cast<int>(in.size()),
                                       const_cast<double*>(in.data()),
                                       fc(out),
                                       FFTW_ESTIMATE);
    fftw_execute(p);
    fftw_destroy_plan(p);
}
 
inline void irfft(const num::CVector& in, int n, num::Vector& out) {
    fftw_plan p = fftw_plan_dft_c2r_1d(n,
                                       const_cast<fftw_complex*>(fc(in)),
                                       out.data(),
                                       FFTW_ESTIMATE);
    fftw_execute(p);
    fftw_destroy_plan(p);
}
 
struct FFTPlanImpl {
    fftw_plan plan;
    int       n;
 
    FFTPlanImpl(int n_, bool forward)
        : n(n_) {
        // Allocate dummy arrays -- FFTW_MEASURE overwrites them during planning
        num::CVector tmp_in(n_, num::cplx{0, 0}), tmp_out(n_, num::cplx{0, 0});
        plan = fftw_plan_dft_1d(n_,
                                fc(tmp_in),
                                fc(tmp_out),
                                forward ? FFTW_FORWARD : FFTW_BACKWARD,
                                FFTW_MEASURE);
    }
 
    ~FFTPlanImpl() {
        fftw_destroy_plan(plan);
    }
 
    void execute(const num::CVector& in, num::CVector& out) const {
        fftw_execute_dft(plan,
                         const_cast<fftw_complex*>(fc(in)),
                         const_cast<fftw_complex*>(fc(out)));
    }
};
 
} // namespace fftw
} // namespace backends
 
#endif // NUMERICS_HAS_FFTW