c++ - fftw:为什么我的 2D DFT 输出与每行的 1D DFT 输出不同？

Question

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

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

在下面的代码片段中，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 。我按照指示分配 nRows * (nCols/2 + 1) 单精度 float here 。对于一维情况，这是在一维 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 转换回时域，同样是一维和二维。我不确定 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 的大小也是错误的，并且您将值复制到其中的方式可能是也错了。