使用 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 以行优先顺序对数组进行操作，因此本征矩阵应该遵守）以及线性访问本征矩阵数据的有效性（对我来说似乎没问题）。