C++ FFTW3的fftw_plan_dft_r2c_2d()和Eigen::MatrixXd的傅里叶正变换_eigen fftw

项目场景：

C++中，对含有较多数据的二维实数矩阵Eigen::MatrixXd进行傅里叶变换时，ftw_plan_dft_r2c_2d()会比fftw_plan_dft_2d()快一些，但在使用ftw_plan_dft_r2c_2d()时，发现要输出与输入同等大小的矩阵时会缺失一些值（共轭值）。

---

问题描述

C++中，对含有较多数据的二维实数矩阵Eigen::MatrixXd进行傅里叶变换的时候，ftw_plan_dft_r2c_2d()如果输出与输入同等大小的矩阵时会缺失一些值（共轭值）。

1、C++使用fftw_plan_dft_2d()对实数矩阵进行傅里叶变换：先将MatrixXd赋值给MatrixXcd，然后再使用fftw_plan_dft_2d()。其输出结果和Matlab的结果相同。

#include <iostream>
#include <cmath>
#include "Eigen/Dense"
#include <opencv2/core/eigen.hpp>
#include "fftw3.h"
 
using namespace std;
using namespace Eigen;
 
int main() {
    MatrixXd a(5, 6);      //奇数行
    a << 11, 18, 20, 16, 10, 14,
          4,  1,  8, 29, 19, 21,
         17, 27, 15,  9, 25,  5,
         22,  7,  2, 24, 30, 23,
          3,  6, 28, 13, 12, 26;
    MatrixXcd A = a;
    MatrixXcd FTa(a.rows(), a.cols());
    fftw_plan P;
    P = fftw_plan_dft_2d(a.cols(), a.rows(), (fftw_complex*)A.data(), (fftw_complex*)FTa.data(), FFTW_FORWARD, FFTW_ESTIMATE);
    fftw_execute(P);
    cout << FTa << endl;
    fftw_destroy_plan(P);
    return 0;
};

2、Matlab的fft2输出结果：

3、但是，在C++使用fftw_plan_dft_r2c_2d()对实数矩阵进行傅里叶变换的时候，其输出结果其实只有Matlab/fftw_plan_dft_2d()结果的前 int(实矩阵行数 / 2) + 1 行的数据。

#include <iostream>
#include <cmath>
#include "Eigen/Dense"
#include <opencv2/core/eigen.hpp>
#include "fftw3.h"
 
using namespace std;
using namespace Eigen;
 
int main() {
    MatrixXd a(5, 6);         //奇数行实矩阵
    a << 11, 18, 20, 16, 10, 14,
          4,  1,  8, 29, 19, 21,
         17, 27, 15,  9, 25,  5,
         22,  7,  2, 24, 30, 23,
          3,  6, 28, 13, 12, 26;
    MatrixXcd FTa(a.rows(), a.cols());
    fftw_plan P;
    P = fftw_plan_dft_r2c_2d(a.cols(), a.rows(), a.data(), (fftw_complex*)FTa.data(), FFTW_ESTIMATE);
    fftw_execute(P);
    cout << FTa << endl;     //输出只有matlab/fftw_plan_dft_2d()输出的前3行数据
    return 0;
};

上图C++输出的红框数据对应下图Matlab输出的红框数据。

4、补充：偶数行实矩阵的情况

int main() {
	MatrixXd aa(4, 6);       //偶数行实矩阵
	aa << 11, 18, 20, 16, 10, 14,
		   4,  1,  8, 29, 19, 21,
		  17, 27, 15,  9, 25,  5,
		  22,  7,  2, 24, 30, 23;
	MatrixXcd FTaa(aa.rows(), aa.cols());
	fftw_plan PP;
	PP = fftw_plan_dft_r2c_2d(aa.cols(), aa.rows(), aa.data(), (fftw_complex*)FTaa.data(), FFTW_ESTIMATE);
	fftw_execute(PP);
	cout << FTaa << endl;
	return 0;
};

---

解决方案：

实际上，fftw_plan_dft_r2c_2d()输出缺少的 ceil(double(实矩阵行数) / 2) - 1 行的数据对应的共轭部分已经全部包含在它的输出里。

以Matlab的输出展示，黄框为fftw_plan_dft_r2c_2d()输出的缺失行，红框为缺失行对应的共轭区域。 因此，只需要对红框的数据取共轭并排列得当获得黄框数据，就能补全fftw_plan_dft_r2c_2d()的输出。

偶数行实矩阵的傅里叶变换的共轭区域与上述奇数行实矩阵的傅里叶变换的共轭区域略有不同，如图所示：

 观察可以发现，互为共轭的数据分为两部分，如蓝框所示：第1列部分与其余列部分。

 第一部分的元素在列方向上共轭对称；第二部分中红框的元素取共轭，并使红框区域翻转180°即为黄框区域。

#include <iostream>
#include <cmath>
#include "Eigen/Dense"
#include <opencv2/core/eigen.hpp>
#include "fftw3.h"
 
using namespace std;
using namespace Eigen;
 
int main(){
    MatrixXd a(5, 6);     //奇数行实矩阵
    a << 11, 18, 20, 16, 10, 14,
          4,  1,  8, 29, 19, 21,
         17, 27, 15,  9, 25,  5,
         22,  7,  2, 24, 30, 23,
          3,  6, 28, 13, 12, 26;
    /*
    MatrixXd a(4, 6);     //偶数行实矩阵
    a << 11, 18, 20, 16, 10, 14,
          4,  1,  8, 29, 19, 21,
         17, 27, 15,  9, 25,  5,
         22,  7,  2, 24, 30, 23; 
    */
 
    MatrixXcd temp(a.rows()/2 + 1, a.cols());
    fftw_plan P;
    P = fftw_plan_dft_r2c_2d(a.cols(), a.rows(), a.data(), (fftw_complex*)temp.data(), FFTW_ESTIMATE);
    fftw_execute(P);
 
    int ConjRows = ceil(double(a.rows()) / 2) - 1;
    MatrixXcd FTa(a.rows(), a.cols());
    FTa.topRows(temp.rows()) = temp;
    FTa.bottomLeftCorner(ConjRows, 1) = ((temp.block(1, 0, ConjRows, 1)).conjugate()).colwise().reverse();
    FTa.bottomRightCorner(ConjRows, a.cols() - 1) = (((temp.block(1, 1, ConjRows, a.cols() - 1)).conjugate()).rowwise().reverse()).colwise().reverse();
    
    cout << FTa << endl;
    fftw_destroy_plan(P);
}