FFT_PFFFT.h

Go to the documentation of this file.

//
//  FFT_PFFFT.h
//  CQTDSP
//
//  Created by Juan Sierra on 4/8/25.
//

 
#pragma once
 
#include <cassert>
 
#include <Eigen/Core>
#include <pffft_double.h>
 
using namespace Eigen;
 
namespace jsa::cicuetea {
class DFTImpl
{
  public:
    DFTImpl(size_t fftSize) :
        fftSize(fftSize)
    {
        workData.resize(2 * fftSize);
        complexSetup = pffftd_new_setup(int(fftSize), PFFFT_COMPLEX);
        realSetup    = pffftd_new_setup(int(fftSize), PFFFT_REAL);
        assert(complexSetup && "Complex setup not initialized properly");
        assert(realSetup && "Real setup not initialized properly");
    }
 
    ~DFTImpl()
    {
        if (complexSetup) pffftd_destroy_setup(complexSetup);
        if (realSetup) pffftd_destroy_setup(realSetup);
    }
 
    void dft(const dcomplex* inPtr, dcomplex* outPtr)
    {
        double* inPtr_  = reinterpret_cast<double*>(const_cast<dcomplex*>(inPtr));
        double* outPtr_ = reinterpret_cast<double*>(outPtr);
        pffftd_transform_ordered(complexSetup, inPtr_, outPtr_, workData.data(), PFFFT_FORWARD);
    }
 
    void idft(const dcomplex* inPtr, dcomplex* outPtr)
    {
        double* inPtr_  = reinterpret_cast<double*>(const_cast<dcomplex*>(inPtr));
        double* outPtr_ = reinterpret_cast<double*>(outPtr);
        pffftd_transform_ordered(complexSetup, inPtr_, outPtr_, workData.data(), PFFFT_BACKWARD);
        Map<ArrayXcd>(outPtr, fftSize) *= (1.0 / fftSize);
    }
 
    void rdft(const double* inPtr, dcomplex* outPtr)
    {
        double* inPtr_  = reinterpret_cast<double*>(const_cast<double*>(inPtr));
        double* outPtr_ = reinterpret_cast<double*>(outPtr);
        pffftd_transform_ordered(realSetup, inPtr_, outPtr_, workData.data(), PFFFT_FORWARD);
    }
 
    void irdft(const dcomplex* inPtr, double* outPtr)
    {
        double* inPtr_  = reinterpret_cast<double*>(const_cast<dcomplex*>(inPtr));
        double* outPtr_ = reinterpret_cast<double*>(outPtr);
        pffftd_transform_ordered(realSetup, inPtr_, outPtr_, workData.data(), PFFFT_BACKWARD);
        Map<ArrayXd>(outPtr, fftSize) *= (1.0 / fftSize);
    }
 
    static std::string getName()
    {
        return "PFFFT";
    }
 
  private:
    const size_t   fftSize;
    Eigen::ArrayXd workData;
    PFFFTD_Setup*  complexSetup = nullptr;
    PFFFTD_Setup*  realSetup    = nullptr;
};
} // namespace jsa::cicuetea