PocketFFT++ 用法:前向/后向变换不返回原始数据

问题描述 投票:0回答:1

我正在尝试使用 PocketFFT++ 进行 2D FFT。

我的示例代码仅采用 8x8 浮点数向量并调用正向变换(PocketFFT++ 函数

r2c
),然后运行逆变换(PocketFFT++ 函数
c2r
)。演示代码如下:

#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;
}

我期望前向/后向变换能够返回原始数据(存在较小的舍入误差)。然而,事实并非如此。显然,我认为这是错误的,即使在阅读了 API 文档后我也无法找出我的错误。

对于逆变换,我将输入视为正向变换的输出(即

vector<complex<float>>
)。

这是上述演示程序的控制台输出:

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,]

我在调用函数时犯了什么错误?

c++ c++11 fft
1个回答
0
投票

我发现了错误。将此发布作为对其他面临同样问题并来到这里的人的答案。

  1. stride_in 和 stride_out 设置不正确。演示代码的正确初始化是:
    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()

  1. r2c 和 c2r 参数如下:
    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);

结果通过

1.0f/64.0f
进行标准化,它是输出矩阵 (
8x8 = 64
) 中条目数的倒数。

现在输入可以正确往返。

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,]
© www.soinside.com 2019 - 2024. All rights reserved.