fftw：为什么我的 2D DFT 输出与每行的 1D DFT 输出不同？
fftw: Why is my 2D DFT output not the same as the output from taking 1D DFT of each row?
对于我的项目，我必须对大型 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
在下面的片段中，this->N = this->nRows * this->nCols。测试一和测试二的结果应该是一样的，但是有明显的不同。我一遍又一遍地阅读文档，仍然找不到它出错的原因。 fftw 进行行主多维变换，矩阵的每行填充 in 。 transfer_output 函数不对值本身做任何事情，只是将标准数组转换为 Eigen::Matrix。我哪里错了？任何帮助将不胜感激。我也试图在这里找到类似的帖子，但 none 就我所能找到的而言遇到了我的问题。
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;
}
这就是我在每种情况下计算 DFT 的方法，使用 fftw 库（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. I allocate nRows * (nCols/2 + 1) single-precision floats as instructed here 上指定的 std::complex。对于 1D 情况，这是在 1D transfer_output 函数中完成的，其中 dft 填充 out[this->nCols - i] for 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 出错了。一维案例有效，所以我只展示二维案例。
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：删除图片，添加内容作为文本。
如果没有任何人都可以编译和测试的完整可重现示例，很难回答这个问题。所以我会给出执行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 的大小也是错误的，并且您将值复制到其中的方式也可能是错误的。
对于我的项目，我必须对大型 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
在下面的片段中，this->N = this->nRows * this->nCols。测试一和测试二的结果应该是一样的，但是有明显的不同。我一遍又一遍地阅读文档，仍然找不到它出错的原因。 fftw 进行行主多维变换，矩阵的每行填充 in 。 transfer_output 函数不对值本身做任何事情，只是将标准数组转换为 Eigen::Matrix。我哪里错了？任何帮助将不胜感激。我也试图在这里找到类似的帖子，但 none 就我所能找到的而言遇到了我的问题。
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;
}
这就是我在每种情况下计算 DFT 的方法，使用 fftw 库（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. I allocate nRows * (nCols/2 + 1) single-precision floats as instructed here 上指定的 std::complex。对于 1D 情况，这是在 1D transfer_output 函数中完成的，其中 dft 填充 out[this->nCols - i] for 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 出错了。一维案例有效，所以我只展示二维案例。
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：删除图片，添加内容作为文本。
如果没有任何人都可以编译和测试的完整可重现示例，很难回答这个问题。所以我会给出执行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 的大小也是错误的，并且您将值复制到其中的方式也可能是错误的。