fft.hpp

Go to the documentation of this file.

/// @file spectral/fft.hpp
/// @brief FFT interface with backend dispatch.
///
/// Each module defines its own backend enum for the choices relevant to it.
/// This enum covers spectral transforms; linalg uses num::Backend in core/policy.hpp.
///
/// Backends:
///   FFTBackend::seq   -- native Cooley-Tukey radix-2 DIT; always available.
///                        Input length must be a power of two.
///   FFTBackend::fftw  -- FFTW3 (AVX / NEON optimised); requires NUMERICS_HAS_FFTW.
///                        Falls back to FFTBackend::seq automatically if FFTW3 is absent.
///
/// Conventions (both backends):
///   Forward DFT:  X[k] = sum_{j=0}^{n-1} x[j] * exp(-2*pi*i*j*k/n)
///   Inverse DFT:  unnormalised (FFTW convention).  Divide by n to recover x.
///   rfft output:  n/2+1 complex bins  (Hermitian symmetry of real input).
#pragma once
 
#include "core/types.hpp"
#include "core/vector.hpp"
 
namespace num {
namespace spectral {
 
/// @brief Selects which backend handles an FFT operation.
enum class FFTBackend {
    seq,     ///< Native Cooley-Tukey radix-2 DIT -- always available (power-of-two only)
    simd,    ///< Handwritten AVX2 / NEON butterfly -- falls back to seq if unavailable
    stdsimd, ///< std::experimental::simd butterfly -- requires NUMERICS_HAS_STD_SIMD
    fftw,    ///< FFTW3 (mixed-radix, AVX / NEON optimised) -- requires NUMERICS_HAS_FFTW
};
 
// Convenience constants -- use these at call sites:
//   fft(in, out, num::spectral::fftw);
inline constexpr FFTBackend seq     = FFTBackend::seq;
inline constexpr FFTBackend fftw    = FFTBackend::fftw;
inline constexpr FFTBackend fft_simd    = FFTBackend::simd;
inline constexpr FFTBackend fft_stdsimd = FFTBackend::stdsimd;
 
/// True when FFTW3 was found at configure time.
inline constexpr bool has_fftw =
#ifdef NUMERICS_HAS_FFTW
    true;
#else
    false;
#endif
 
/// True when handwritten AVX2 or NEON backend is available.
inline constexpr bool has_fft_simd =
#if defined(NUMERICS_HAS_AVX2) || defined(NUMERICS_HAS_NEON)
    true;
#else
    false;
#endif
 
/// True when std::experimental::simd is available.
inline constexpr bool has_fft_stdsimd =
#ifdef NUMERICS_HAS_STD_SIMD
    true;
#else
    false;
#endif
 
/// Automatically selected at configure time: fftw > simd > seq.
inline constexpr FFTBackend default_fft_backend =
#ifdef NUMERICS_HAS_FFTW
    FFTBackend::fftw;
#elif defined(NUMERICS_HAS_AVX2) || defined(NUMERICS_HAS_NEON)
    FFTBackend::simd;
#else
    FFTBackend::seq;
#endif
 
// -- One-shot transforms --------------------------------------------------
 
/// @brief Forward complex DFT.  out must be pre-allocated to in.size().
void fft(const CVector& in, CVector& out,
         FFTBackend b = default_fft_backend);
 
/// @brief Inverse complex DFT (unnormalised: result = n * true_inverse).
void ifft(const CVector& in, CVector& out,
          FFTBackend b = default_fft_backend);
 
/// @brief Real-to-complex forward DFT.  out must be pre-allocated to n/2+1.
void rfft(const Vector& in, CVector& out,
          FFTBackend b = default_fft_backend);
 
/// @brief Complex-to-real inverse DFT (unnormalised).
/// @param n  Length of the real output (must satisfy in.size() == n/2+1).
void irfft(const CVector& in, int n, Vector& out,
           FFTBackend b = default_fft_backend);
 
// -- Reusable plan --------------------------------------------------------
 
/// @brief Precomputed plan for repeated complex transforms of a fixed size.
///
/// @code
/// num::spectral::FFTPlan plan(1024);          // forward, default backend
/// for (auto& frame : frames)
///     plan.execute(frame, spectrum);           // O(n log n), no allocation
/// @endcode
class FFTPlan {
public:
    /// @param n        Transform size (must be power of two for FFTBackend::seq)
    /// @param forward  true = forward DFT, false = inverse DFT
    /// @param b        Backend to use (default: default_fft_backend)
    explicit FFTPlan(int n, bool forward = true,
                     FFTBackend b = default_fft_backend);
    ~FFTPlan();
 
    FFTPlan(const FFTPlan&)            = delete;
    FFTPlan& operator=(const FFTPlan&) = delete;
 
    void execute(const CVector& in, CVector& out) const;
 
    int        size()    const { return n_; }
    FFTBackend backend() const { return backend_; }
 
private:
    int        n_;
    FFTBackend backend_;
    void*      impl_;   // backends::seq::FFTPlanImpl or backends::fftw::FFTPlanImpl
};
 
} // namespace spectral
} // namespace num