mercredi 4 octobre 2023

PocketFFT++ usage: Forward/Backward transforms not giving back original data

I am trying to use PocketFFT++ for a 2D FFT.

My sample code simply takes an 8x8 vector of floats and calls the forward transform (PocketFFT++ function r2c) and then runs the inverse transform (PocketFFT++ function c2r). The demo code is given below:

#include <complex>
#include <cmath>
#include <vector>
#include <iostream>

#include "..\..\pocketfft_hdronly.h"

template<typename T>
std::ostream& operator<<(std::ostream& os, const std::vector<T>& vec)
{
    os << '[';
    for (auto val : vec)
        os << val << ',';
    os << ']';

    return os;
}

int main()
{
    std::cout << "PocketFFT++ test" << std::endl;

    using namespace std;
    using namespace pocketfft;

    /*
    shape_t shape_in;               // dimensions of the input shape
    stride_t stride_in;             // must have the size of each element. Must have size() equal to shape_in.size()
    stride_t stride_out;            // must have the size of each element. Must have size() equal to shape_in.size()
    shape_t axes;                   // 0 to shape.size()-1 inclusive
    bool forward;                   // FORWARD or BACKWARD
    float* data_in;                 // input data (reals)
    complex<float>* data_out;       // output data (FFT(input))
    float fct;                      // scaling factor

    r2c(const shape_t & shape_in,
        const stride_t & stride_in, const stride_t & stride_out, const shape_t & axes,
        bool forward, const T * data_in, complex<T> *data_out, T fct,
        size_t nthreads = 1)
    */

    shape_t shape_in{8,8};                                              // dimensions of the input shape
    stride_t stride_in{sizeof(float),sizeof(float)};                    // must have the size of each element. Must have size() equal to shape_in.size()
    stride_t stride_out{sizeof(complex<float>),sizeof(complex<float>)}; // must have the size of each element. Must have size() equal to shape_in.size()
    shape_t axes{0,1};                                                  // 0 to shape.size()-1 inclusive
    bool forward{ FORWARD };                                            // FORWARD or BACKWARD
    vector<float> data_in{
        1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f,
        2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 1.0f,
        3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 1.0f, 2.0f,
        4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 1.0f, 2.0f, 3.0f,
        5.0f, 6.0f, 7.0f, 8.0f, 1.0f, 2.0f, 3.0f, 4.0f,
        6.0f, 7.0f, 8.0f, 1.0f, 2.0f, 3.0f, 4.0f, 5.0f,
        7.0f, 8.0f, 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 7.0f,
        8.0f, 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 7.0f, 1.0f
    };                                                                  // input data (reals)
    vector<complex<float>> data_out(64);                                // output data (FFT(input))
    float fct{ 1.0f };                                                  // scaling factor

    std::cout << "data_in: " << data_in << std::endl;

    r2c(
        shape_in, 
        stride_in, 
        stride_out, 
        axes,
        forward, 
        data_in.data(), 
        data_out.data(), 
        fct
    );

    std::cout << "data_out: " << data_out << std::endl;

    /* inverse
    template<typename T> void c2r(const shape_t &shape_out,
      const stride_t &stride_in, const stride_t &stride_out, const shape_t &axes,
      bool forward, const complex<T> *data_in, T *data_out, T fct,
      size_t nthreads=1)
    */

    shape_t &inv_shape_out = shape_in;
    stride_t &inv_stride_in = stride_out;
    stride_t &inv_stride_out = stride_in;
    shape_t inv_axes = {0, 1};
    bool inv_forward = BACKWARD;
    vector<complex<float>>& inv_data_in = data_out;
    vector<float> inv_data_out(data_in.size());
    float inv_fct = 1.0f;

    std::cout << "inv_data_in: " << inv_data_in << std::endl;

    c2r(inv_shape_out,
        inv_stride_in,
        inv_stride_out,
        inv_axes,
        inv_forward,
        inv_data_in.data(),
        inv_data_out.data(),
        inv_fct);

    std::cout << "inv_data_out: " << inv_data_out << std::endl;
}

I was expecting the forward/backward transforms to give me back the original data (subject to minor rounding errors). However, that is not the case. Clearly, I am calling this wrong and I am not able to figure out my mistake even after reading the API documentation.

For the inverse transform, I am considering the input to be the output of the forward transform (i.e., vector<complex<float>>).

Here's the console output from the above demo program:

PocketFFT++ test
data_in: [1,2,3,4,5,6,7,8,2,3,4,5,6,7,8,1,3,4,5,6,7,8,1,2,4,5,6,7,8,1,2,3,5,6,7,8,1,2,3,4,6,7,8,1,2,3,4,5,7,8,1,2,3,4,5,7,8,1,2,3,4,5,7,1,]
data_out: [(316,0),(-25,9.65685),(-4,71.598),(127.735,160.049),(-235.598,539.127),(32.5685,-175.029),(-11.5979,1.65685),(135.764,45.3238),(134.108,105.421),(59.3827,71.7645),(30.3431,-110.912),(245.019,-516.617),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),]
inv_data_in: [(316,0),(-25,9.65685),(-4,71.598),(127.735,160.049),(-235.598,539.127),(32.5685,-175.029),(-11.5979,1.65685),(135.764,45.3238),(134.108,105.421),(59.3827,71.7645),(30.3431,-110.912),(245.019,-516.617),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),(0,0),]
inv_data_out: [2055.96,-2332.6,-1000.79,-1000.81,-54.1178,892.817,958.362,2704.16,-765.797,1036.13,-93.4214,385.087,-565.421,-836.476,4521.54,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,]

What mistake am I making in the calling of the functions?

Aucun commentaire:

Enregistrer un commentaire