fftw：为什么我的 2D DFT 输出与每行的 1D DFT 输出不一样？

Question

对于我的项目，我必须对大型 2D 输入矩阵进行 DFT，对其进行处理，然后使用 IDFT 将其转换回来，并将结果与​​输入矩阵进行比较。我的问题在于 2D DFT 步骤。我用一个小的简单数据集编写了一个测试，我在main() 中执行它。我将Eigen 库用于矩阵和向量。输出是这样的：

Using testM in a myTransform object:
1 2 3
4 5 6
7 8 9
calculateDFT took 33 Microseconds
DFT
          (6,0)
(-1.5,0.866025)
(-1.5,0.866025)
         (15,0)
(-1.5,0.866025)
(-1.5,0.866025)
         (24,0)
(-1.5,0.866025)
(-1.5,0.866025)
IDFT
1 2 3 4 5 6 7 8 9
Using testM in a myTransform2D object:
1 2 3
4 5 6
7 8 9
Default myTransform object created
DFT2D
          (45,0)  (-4.5,-2.59808)          (45,-0)
          (45,0) (-13.5,-7.79423)          (45,-0)
          (45,0)           (0,-0)          (45,-0)
IDFT
27.5 -0.5 -1.5 -8.5  8.5  7.5   -7   10    9

在下面的 sn-ps 中，this->N = this->nRows * this->nCols。测试 1 和测试 2 的结果应该是相同的，但它们明显不同。我一遍又一遍地阅读文档，但仍然找不到出错的原因。 fftw 进行行主要的多维变换，in 填充矩阵的每一行。 transfer_output 函数不对值本身做任何事情，只是将标准数组转换为 Eigen::Matrix。我哪里错了？任何帮助将不胜感激。我也试图在这里找到类似的帖子，但据我所知，没有一个有我的问题。

void test()
{
    RowVectorXf test(9);
    test << 1, 2, 3, 4, 5, 6, 7, 8, 9;
    // Prep matrix 
    Map<Matrix<float, 3, 3, RowMajor>> testM(test.data()); // convert test to rowmajor matrix

    // Test 1: feed the matrix to a myTransform object and take 1D DFTs and 1D IDFTs 
    std::cout << "Using testM in a myTransform object:\n" << testM << std::endl;
    myTransform testX1D(testM, 0);
    testX1D.vectorise();
    testX1D.calculateDFT();
    testX1D.calculateIDFT();
    std::cout << "DFT" << std::endl << testX1D.dft << std::endl;
    std::cout << "IDFT" << std::endl << testX1D.idft << std::endl; // works, too.

    .. Test 2: Feed the matrix to a myTransform2D object and take the 2D DFT and IDFT.
    std::cout << "Using testM in a myTransform2D object:\n" << testM << std::endl;
    myTransform2D testX(testM, 0); // 2D version 
    testX.vectorise(); // stored in testX.m which will hold the same as test but in a colmajor vector.
    testX.calculateDFT(); // where it goes wrong?
    std::cout << "DFT2D" << std::endl << testX.dft2D << std::endl;
    testX.calculateIDFT();
    std::cout << "IDFT" << std::endl << testX.idft << std::endl;
}

这就是我在每种情况下使用 fftw 库计算 DFT 的方式（fftwf 因为我使用单精度来节省内存，并且非测试数据的值在 -10000 到 10000 的数量级上，所以我认为这不是问题）。

void myTransform::calculateDFT()
/// Calculates discrete fourier transform of vectorised data `m`.
/** uses the FFTW library (https://fftw.org). The dft is stored in myTransform::dft*/
{
    //std::cout << m << std::endl;
    fftwf_complex* out;
    fftwf_plan p;
    out = (fftwf_complex*)fftw_malloc(sizeof(fftwf_complex) * this->nCols);
    float* in = new float[static_cast<const float&>(this->nCols)];
    p = fftwf_plan_dft_r2c_1d(this->nCols, in, out, FFTW_ESTIMATE);
    // calculate DFT for each trace and assign it to a segment of this->dft
    unsigned int factor = 0;
    auto check = std::chrono::high_resolution_clock::now();
    for (int k = 0; k < this->nRows; k++)
    {
        factor = k * this->nCols;
        //TODO: if possible, fix this slow in[i] = ... part. 
        for (int i = 0; i < this->nCols; i++)
        {
            in[i] = this->m[factor + i];
        }
        p = fftwf_plan_dft_r2c_1d(this->nCols, in, out, FFTW_ESTIMATE);
        fftwf_execute(p);
        this->transfer_output(out, k); // does nothing but add the output to a vector dft. 
    }
    delete [] in;
    fftwf_free(out);
    fftwf_destroy_plan(p);
}

对于 2D DFT 案例：这里我使用fftw.org 中指定的 std::complex。我按照here 的指示分配nRows * (nCols/2 + 1) 单精度浮点数。对于一维情况，这是在一维transfer_output 函数中完成的，其中dft 填充有out[this->nCols - i] 用于i > this->nCols/2

void myTransform2D::calculateDFT()
/// Should calculate the DFT in 2D with fftwf_plan_dft_r2c_2d(n0, n1, *in, *out, flags).
{
    std::complex<float>* out;
    fftwf_plan p;
    out = (std::complex<float>*)fftwf_malloc(sizeof(std::complex<float>) * this->nRows * (this->nCols/2+1)); // Hermitian symmetry for r2c transforms
    float* in = new float[this->N];
    in = (float*)fftwf_malloc(sizeof(float) * this->N);
    p = fftwf_plan_dft_r2c_2d(this->nRows, this->nCols, in, reinterpret_cast<fftwf_complex*>(out), FFTW_ESTIMATE);
    // Fill input array
    for (int i = 0; i < this->nRows; i++)
    {
        int factor = i * this->nCols;
        for (int j = 0; j < this->nCols; j++)
        {
            in[factor + j] = this->m[factor + j];
        }
    }
    fftwf_execute(p);
    transfer_output(out);
    fftwf_free(in);
    fftwf_free(out);
    fftwf_destroy_plan(p);
}

我使用 IDFT 转换回时域，同样是 1D 和 2D。我不确定 2D 版本是否有效，因为 DFT 出错了。 1D 案例有效，所以我只展示 2D 案例。

void myTransform2D::calculateIDFT()
/// Calculates inverse fourier transform of `this->dft2D`.
/** Also uses the FFTW library. Results might not be perfect as we use floats
    instead of doubles because of large data sizes. */
{
    float* out = new float[this->N];
    std::complex<float>* in; 
    fftwf_plan pr;
    in = (std::complex<float>*)fftwf_malloc(sizeof(std::complex<float>) * this->N);
    out = (float*)fftwf_malloc(sizeof(float) * this->N);
    pr = fftwf_plan_dft_c2r_2d(this->nRows, this->nCols, reinterpret_cast<fftwf_complex*>(in), out, FFTW_ESTIMATE);

    for (int i = 0; i < this->nRows; i++)
    {
        for (int j = 0; j < this->nCols; j++)
        {
            in[i * this->nCols + j] = this->dft2D(i, j);
        }
    }

    fftwf_execute(pr);

    for (int i = 0; i < this->N; i++)
    {
        this->idft[i] = out[i] / this->N; // fftw does unnormalized inverse transforms.
    }

    fftwf_free(out);
    fftwf_free(in);
    fftwf_destroy_plan(pr);
}

编辑：按照建议删除了一些代码 EDIT2：删除图像，添加内容为文本。

Answer 1

如果没有任何人都可以编译和测试的完整可复制示例，很难回答这个问题。所以我会给出执行2D前向和后向变换的代码和参考结果。

转发。

const auto fft_size = n * (n / 2 + 1);
auto out = (std::complex<float>*)fftwf_malloc(sizeof(std::complex<float>) * fft_size);
auto p = fftwf_plan_dft_r2c_2d(n, n, in, (fftwf_complex*)(out), FFTW_ESTIMATE);
fftwf_execute(p);

for (std::size_t i = 0; i < fft_size; ++i)
    std::cout << *(out + i) << ' ';

对于行优先顺序的矩阵

1 2 3
4 5 6
7 8 9

正确的输出是：

(45,0) (-4.5,2.59808) (-13.5,7.79423) (0,0) (-13.5,-7.79423) (0,0)

向后。

auto in2 = (float*)fftwf_malloc(sizeof(float) * n * n);
auto p2 = fftwf_plan_dft_c2r_2d(n, n, (fftwf_complex*)(out), in2, FFTW_ESTIMATE);
fftwf_execute(p2);

for (std::size_t row = 0; row < n; ++row) {
    for (std::size_t col = 0; col < n; ++col)
        std::cout << *(in2 + col + row * n) / (n * n) << ' ';
    std::cout << std::endl;
}

这会输出原始矩阵。

---

请注意，前向变换 (fft_size) 的输出大小是 n * (n / 2 + 1) 而不是 n * n。在您的输出中，我看到 9 个复杂条目而不是 6 个。函数 calculateIDFT() 中 in 的大小也是错误的，并且您将值复制到其中的方式也可能是错误的。