使用 Eigen 和 FFTW 进行二维傅里叶变换

Question

我正在尝试使用 FFTW 进行实值二维傅里叶变换。我的数据存储在动态大小的特征矩阵中。这是我写的包装类：

FFT2D.h：

#include <Eigen>

class FFT2D {
public:

    enum FFT_TYPE {FORWARD=0, REVERSE=1};
    FFT2D(EMatrix &input, EMatrix &output, FFT_TYPE type_ = FORWARD);
    ~FFT2D();

    void execute();

private:
    EMatrix& input;
    EMatrix& output;
    fftw_plan plan;
    FFT_TYPE type;
};

FFT2D.cpp：

#include "FFT2D.h"
#include <fftw3.h>
#include "Defs.h"


FFT2D::FFT2D(EMatrix &input_, EMatrix &output_, FFT_TYPE type_)
        : type(type_), input(input_), output(output_) {

    if (type == FORWARD)
        plan = fftw_plan_dft_2d((int) input.rows(), (int) input.cols(),
                (fftw_complex *) &input(0), (fftw_complex *) &output(0),
                FFTW_FORWARD, FFTW_ESTIMATE);
    else
        // placeholder for ifft-2d code, unwritten
}


FFT2D::~FFT2D() {
    fftw_destroy_plan(plan);
}

void FFT2D::execute() {
    fftw_execute(plan);  // seg-fault here
}

还有

EMatrix

的定义：

typedef Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> EMatrix;

问题是，我在

FFT2D::execute()

中遇到了段错误。我知道我在构造函数中设置了错误，并且我尝试了多种不同的方法，但我似乎无法在这方面取得任何进展。

我尝试过的事情包括：将

EMatrix

typedef 更改为

Eigen::ColMajor

，将

(fftw_complex *) input.data()

传递给

fftw_plan_dft_2d

，使用不同的 fftw 计划 (

fftw_plan_dft_r2c_2d

)。

我的 C++（显然）很生疏，但最终我需要的是对双精度实值 2D 本征矩阵进行 2D FT。预先感谢。

Answer 1

这里的主要问题是不存在“实值傅里叶变换”这样的东西。它只是虚部为零的傅里叶变换，但零仍然必须存在，正如您从

fftw_complex

定义中看到的：

typedef double fftw_complex[2];

这是有道理的，因为输出可以（并且可能）具有非零虚部。
_{输出将具有一些对称属性，即在一维变换的情况下它将是一个偶函数。}

因此，

(fftw_complex *) &input(0)

转换实际上不起作用 - FFTW 期望传递给它的

double

值是您传递给它的两倍。

解决方案是将矩阵原始数据与零交错，有多种方法可以做到这一点。几个例子：

• 您可以将整个矩阵复制到一个新数组中，然后将其传递给 FFTW，并在过程中添加零。
• 您可以在矩阵本身中保留零的空间 - 这样您就可以避免复制，但可能需要大量重构:)
• 我能想到的最好方法是使用

  std::complex<double>

  作为标量。这会在一定程度上损害您对“实值 FFT”的注意，但首先几乎不存在这样的事情。相反，您将能够保持所有实值操作不变，并且

  std::complex

  的布局将完美适合

  fftw_complex

  。

这里可能还有其他一些事情需要考虑，比如存储顺序（FFTW 按行主序对数组进行操作，因此本征矩阵应该遵守）以及对本征矩阵数据的线性访问的有效性（对我来说似乎没问题）。

Answer 2

如果您的输入是真实的，您可以使用：

fftw_plan fftw_plan_dft_r2c_2d(int n0, int n1, 双*输入，fftw_complex *输出， 未签名的标志）；

（这可以避免将实数转换为复数）。