如何将 FFTW 与 Eigen 库一起使用？

Question

我正在尝试学习如何将 FFTW 库与 Eigen 一起使用。我不想使用 Eigen 不受支持的模块，因为我最终想将 FFTW 的智慧功能合并到我的代码中。但是，我正在努力实现一个非常基本的示例。到目前为止，这是我的代码：

void fft(Eigen::Ref<Eigen::VectorXcd> inVec, int N) {

    fftw_complex *in, *out;
    fftw_plan p;

    in = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    p = fftw_plan_dft_1d(N, in, in, FFTW_FORWARD, FFTW_ESTIMATE);

    in = reinterpret_cast<fftw_complex*>(inVec.data());

    fftw_execute(p);

    fftw_destroy_plan(p);

    // reassign input back to inVec here

    fftw_free(in); fftw_free(out);
}

这基本上是从 fftw 的 documentation 中的第 2.1 章复制而来的。文档上好像说它的基本接口不能对不同的数据集进行操作，但是你也必须在创建计划后初始化数据。我不明白这一点。为什么不简单地用新数据覆盖您第一次初始化的数据并再次执行计划？

此外，我尝试将特征向量转换为 fftw_complex，但我怀疑我在这里犯了一个错误。在 Eigen 不支持的 FFTW 模块中，还有一个 const_cast 调用。为什么会在那里，这是必要的，因为我的基础数据在这里不是 const 吗？

最后，如果 in 是指向 fftw_complex 数据的指针，我怎样才能将它重新分配给我的 inVec 然后释放 in？

Answer 1

I don't understand this point. Why not simply overwrite the data you initialize the first time with new data and execute the plan again?

是的，这正是 fftw 希望您做的。 in = reinterpret_cast<fftw_complex*>(inVec.data()); 行只是设置了一个指针。它不会复制数组。你需要把内容memcpy过来，意思是memcpy(in, invec.data(), N * sizeof(fftw_complex));

您想要的（在 FFTW 文档中有些隐藏）是 "new-array execute function"，它允许您为调用指定不同的数组。但是，请注意文档中的这一行：

The alignment issue is especially critical, because if you don’t use fftw_malloc then you may have little control over the alignment of arrays in memory. For example, neither the C++ new function nor the Fortran allocate statement provide strong enough guarantees about data alignment. If you don’t use fftw_malloc, therefore, you probably have to use FFTW_UNALIGNED (which disables most SIMD support). If possible, it is probably better for you to simply create multiple plans (creating a new plan is quick once one exists for a given size), or better yet re-use the same array for your transforms.

当您使用 Eigen 时，这可能不是问题，因为 Eigen 也会对齐其数组。 GCC 还使用 16 字节对齐，因此即使在调用 malloc 时也可能没问题，但这取决于您的平台。但是，您的函数接口不保证对齐，因为 Eigen::Ref 可能是更大数组的一部分。但是你可以在运行时检查。像这样：

unsigned flags = FFTW_ESTIMATE;
if(reinterpret_cast<size_t>(invec.data()) % 16)
    flags |=  FFTW_UNALIGNED;
p = fftw_plan_dft_1d(..., flags);

有关矢量化丢失的警告可能已过时。今天，正确对齐已不那么重要（尤其是在任何支持 AVX 的东西上），我怀疑（但尚未验证）FFTW 也会对未对齐的内存访问进行矢量化。

旁注：在 p = fftw_plan_dft_1d(N, in, in, FFTW_FORWARD, FFTW_ESTIMATE); 行中，第三个参数应该是 out，而不是 in；除非你想要一个 in-place 转换，在这种情况下你不需要分配 out 数组。

编辑：对齐检查可能会中断。我刚刚检查了源代码，FFTW 可能会根据编译标志定义不同的对齐方式。我没有找到一种好方法来弄清楚它使用的对齐方式。对于 AVX 可能是 32 字节，对于 AVX-512 可能是 64 字节。

另一方面，提供用于检查此类内容的 fftw_alignment_of 函数执行简单的模 16，就像我上面的代码一样。所以我不知道。破损应该非常明显。它只会崩溃，不会导致无效结果。

如何将 FFTW 与 Eigen 库一起使用？

How can I use FFTW with the Eigen library?

c++

eigen

fftw