使用PocketFFT++执行2D FFT正逆变换无法还原原始数据的问题

核心错误点 #

1. 步长（stride）参数设置错误 #

PocketFFT的stride参数表示每个维度上相邻元素的字节偏移量，而非单个元素的大小。你的代码中：

• 实输入数组的stride_in设置为{sizeof(float), sizeof(float)}，第二个维度（行间）的步长错误，正确值应为一行的总字节数：8 * sizeof(float)（8个float元素的大小）。
• 复数输出数组的stride_out设置为{sizeof(complex<float>), sizeof(complex<float>)}，同样第二个维度步长错误，正确值应为(8/2 +1) * sizeof(complex<float>)（r2c输出的每行元素数是N/2+1）。

2. r2c输出数组的尺寸错误 #

实转复FFT（r2c）利用实信号FFT的共轭对称性，仅存储正频率部分，因此输出尺寸应为shape[0] * (shape[1]/2 +1)。对于8×8的实输入，正确的输出复数数组大小是8*(8/2+1)=40，而非你设置的64。多余的24个默认初始化的0元素会完全破坏后续逆变换的输入数据。

3. 逆变换缩放因子缺失 #

PocketFFT不自动进行归一化，正逆变换的缩放因子需要手动配置：

• 正变换（r2c）使用fct=1.0f是正确的。
• 逆变换（c2r）需要设置fct=1.0f/(8*8)=1/64.0f，否则逆变换结果会是原始数据的总元素数倍（64倍）。

4. c2r调用时的输入形状与步长不匹配 #

逆变换c2r的输入是r2c的输出，其形状为{8, 5}，因此需要对应设置输入步长，不能直接复用错误的stride_out参数。

修正后的示例代码 #

#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 << '[';
    size_t cnt = 0;
    for (auto val : vec) {
        os << val;
        if (++cnt != vec.size()) os << ',';
    }
    os << ']';
    return os;
}

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

    using namespace std;
    using namespace pocketfft;

    // 正向r2c变换配置
    shape_t shape_in{8,8};
    // 实输入步长：行内元素偏移sizeof(float)，行间偏移8*sizeof(float)
    stride_t stride_in{sizeof(float), 8 * sizeof(float)};
    // 复数输出形状：8行，每行8/2+1=5个元素
    shape_t shape_out_r2c{8, 8/2 +1};
    // 复数输出步长：行内偏移sizeof(complex<float>)，行间偏移5*sizeof(complex<float>)
    stride_t stride_out{sizeof(complex<float>), (8/2 +1) * sizeof(complex<float>)};
    shape_t axes{0,1};
    bool forward{ FORWARD };
    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
    };
    // 正确尺寸的复数输出数组
    vector<complex<float>> data_out(8*(8/2+1));
    float fct{1.0f};

    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 (r2c result): " << data_out << std::endl;

    // 逆向c2r变换配置
    shape_t inv_shape_out = shape_in; // 逆变换输出形状和原输入一致
    stride_t inv_stride_in = stride_out; // 逆变换输入步长和r2c输出步长一致
    stride_t inv_stride_out = stride_in; // 逆变换输出步长和原输入步长一致
    shape_t inv_axes = {0,1};
    bool inv_forward = BACKWARD;
    vector<float> inv_data_out(data_in.size());
    // 逆变换缩放因子：1/(总元素数)
    float inv_fct = 1.0f/(8*8);

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

    std::cout << "inv_data_out (c2r result): " << inv_data_out << std::endl;

    // 验证误差
    cout << "\nError check:" << endl;
    for (size_t i=0; i<data_in.size(); ++i) {
        float err = abs(inv_data_out[i] - data_in[i]);
        cout << "Index " << i << ": err = " << err << endl;
    }
}

修正后效果 #

修正后的代码执行后，逆变换输出inv_data_out会与原始data_in几乎完全一致，仅存在浮点运算带来的微小舍入误差（通常在1e-6量级以内）。

内容的提问来源于stack exchange，提问作者vvg