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?
I found the errors. Posting this as an answer for anyone else who faces the same issue and comes here.
stride_in and stride_out have not been set correctly. The correct initialization for the demo code is:stride_t stride_in{sizeof(float),sizeof(float)*8}; // 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>)*8}; // must have the size of each element. Must have size() equal to shape_in.size()
r2c(shape_in, stride_in, stride_out, axes, forward, data_in.data(), data_out.data(), fct);
std::cout << "data_out: " << data_out << std::endl;
vector<complex<float>>& inv_data_in = data_out;
vector<float> inv_data_out(data_in.size());
float inv_fct = 1.0f/64.0f;
std::cout << "inv_data_in: " << inv_data_in << std::endl;
c2r(shape_in, stride_out, stride_in, axes, BACKWARD, data_out.data(), inv_data_out.data(), inv_fct);
The result is normalized by 1.0f/64.0f which is the reciprocal of the number of entries in the output matrix (8x8 = 64).
Now the input can be roundtripped correctly.
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: [(284,0),(-3.53553,-2.53553),(-1,-5),(3.53553,-4.53553),(6,0),(3.53553,4.53553),(-1,5),(-3.53553,2.53553),(-3.53553,-2.53553),(-33.4142,72.6691),(2.53553,-4.94975),(6,-1.41421),(4.94975,3.94975),(0,6),(-3.94975,3.53553),(-4.58579,-1.43051e-06),(-1,-5),(2.53553,-4.94975),(-26,30),(5.94975,3.53553),(1,7),(-4.53553,4.94975),(-6,0),(-3.94975,-3.53553),(3.53553,-4.53553),(6,-1.41421),(5.94975,3.53553),(-30.5858,20.669),(-4.94975,5.94975),(-7.41421,-2.98023e-07),(-4.53553,-4.94975),(9.53674e-07,-6),(6,0),(4.94975,3.94975),(1,7),(-4.94975,5.94975),(-40,0),(-4.94975,-5.94975),(1,-7),(4.94975,-3.94975),(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: [(284,0),(-3.53553,-2.53553),(-1,-5),(3.53553,-4.53553),(6,0),(3.53553,4.53553),(-1,5),(-3.53553,2.53553),(-3.53553,-2.53553),(-33.4142,72.6691),(2.53553,-4.94975),(6,-1.41421),(4.94975,3.94975),(0,6),(-3.94975,3.53553),(-4.58579,-1.43051e-06),(-1,-5),(2.53553,-4.94975),(-26,30),(5.94975,3.53553),(1,7),(-4.53553,4.94975),(-6,0),(-3.94975,-3.53553),(3.53553,-4.53553),(6,-1.41421),(5.94975,3.53553),(-30.5858,20.669),(-4.94975,5.94975),(-7.41421,-2.98023e-07),(-4.53553,-4.94975),(9.53674e-07,-6),(6,0),(4.94975,3.94975),(1,7),(-4.94975,5.94975),(-40,0),(-4.94975,-5.94975),(1,-7),(4.94975,-3.94975),(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: [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,]