使用 FFTW_MEASURE 时 FFT 输出为空，但使用 FFTW_ESTIMATE 时工作正常

Question

我在尝试使用 fftw3 时遇到以下问题。由于某种原因，每当我使用 FFTW_MEASURE 而不是 FFTW_ESTIMATE 进行 FFT 时，我都会得到空白输出。最终我试图实现 fft 卷积，因此下面的示例包括 FFT 和逆 FFT。

显然我错过了一些东西......有人能教育我吗？谢谢你！

我使用的是 Linux (OpenSUSE Leap 42.1)，使用包管理器中提供的 fftw3 版本。

最小工作示例：

#include <iostream>
#include <iomanip>
#include <cmath>
#include <fftw3.h>
using namespace std;


int main(int argc, char ** argv)
{
    int width = 10;    
    int height = 8;

    cout.setf(ios::fixed|ios::showpoint);
    cout << setprecision(2);

    double * inp = (double *) fftw_malloc(sizeof(double) * width * height);
    fftw_complex * cplx = (fftw_complex *)fftw_malloc(sizeof(fftw_complex) * height * (width/2 + 1));

    for(int i = 0; i < width * height; i++) inp[i] = sin(i);

    fftw_plan fft = fftw_plan_dft_r2c_2d(height, width, inp, cplx,  FFTW_MEASURE );
    fftw_plan ifft = fftw_plan_dft_c2r_2d(height, width, cplx, inp,  FFTW_MEASURE );

    fftw_execute(fft);

    for(int j = 0; j < height; j++)
    {
        for(int i = 0; i < (width/2 + 1); i++)
        {
            cout << cplx[i+width*j][0] << " ";
        }
        cout << endl;
    }
    cout << endl << endl;

    fftw_execute(ifft);

    for(int j = 0; j < height; j++)
    {
        for(int i = 0; i < width; i++)
        {
            cout << inp[i+width*j] << " ";
        }
        cout << endl;
    }

    fftw_destroy_plan(fft);
    fftw_destroy_plan(ifft);
    fftw_free(cplx);
    fftw_free(inp);

    return 0;
}

只需在 FFTW_ESTIMATE 和 FFTW_MEASURE 之间进行更改即可。

编译为：

g++ *.cpp -lm -lfftw3 --std=c++11

使用 FFTW_ESTIMATE 输出（第一个块是 FT 的实部，第二个块是逆 FT 之后）：

1.51 2.24 -1.52 -0.05 0.15 0.19 
0.23 0.15 1.77 1.19 0.54 0.41 
1.97 -0.15 -1.32 -2.51 -1.20 -3.38 
4.34 15.21 -24.82 -7.44 -4.16 -2.51 
-0.43 -0.06 1.55 2.93 -2.81 -0.42 
0.00 0.00 0.00 -nan 0.00 0.00 
0.00 0.00 0.00 0.00 0.00 -nan 
0.00 0.00 0.00 0.00 0.00 0.00 


0.00 67.32 72.74 11.29 -60.54 -76.71 -22.35 52.56 79.15 32.97 
-43.52 -80.00 -42.93 33.61 79.25 52.02 -23.03 -76.91 -60.08 11.99 
73.04 66.93 -0.71 -67.70 -72.45 -10.59 61.00 76.51 21.67 -53.09 
-79.04 -32.32 44.11 79.99 42.33 -34.25 -79.34 -51.48 23.71 77.10 
59.61 -12.69 -73.32 -66.54 1.42 68.07 72.14 9.89 -61.46 -76.30 
-20.99 53.62 78.93 31.67 -44.70 -79.98 -41.72 34.89 79.43 50.94 
-24.38 -77.29 -59.13 13.39 73.60 66.15 -2.12 -68.44 -71.83 -9.18 
61.91 76.08 20.31 -54.14 -78.81 -31.02 45.29 79.96 41.12 -35.53

使用 FFTW_MEASURE 输出（第一个块是 FT 的实部，第二个块是逆 FT 之后）：

0.00 0.00 0.00 0.00 0.00 0.00 
0.00 0.00 0.00 0.00 0.00 0.00 
0.00 0.00 0.00 0.00 0.00 0.00 
0.00 0.00 0.00 0.00 0.00 0.00 
0.00 0.00 0.00 0.00 0.00 0.00 
0.00 0.00 0.00 -nan 0.00 0.00 
0.00 0.00 0.00 0.00 0.00 -nan 
0.00 0.00 0.00 0.00 0.00 0.00 


0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

Answer 1

@Paul_R 的评论。足以解决问题。输入数组可以按fftw_plan_dft_r2c_2d()调用进行修改。因此，输入数组必须在创建 fftw 计划后初始化。

FFTW 的规划器标志的文档详细说明了正在发生的情况。FFTW_ESTIMATE我很确定您已经猜到了保留输入数组并修改它的原因FTTW_MEASURE。

重要提示：规划器会在规划期间覆盖输入数组，除非已保存的计划（请参阅 Wisdom）可用于该问题，因此您应该在创建计划后初始化输入数据。*** 唯一的例外是 FFTW_ESTIMATE 和 FFTW_WISDOM_ONLY 标志，如下所述。

...

• FFTW_ESTIMATE指定，不是使用不同算法的实际测量，而是使用简单的启发式来快速选择（可能是次优的）计划。使用此标志，在规划期间不会覆盖输入/输出数组。
• FFTW_MEASURE告诉 FFTW 通过实际计算多个 FFT 并测量其执行时间来找到优化计划。根据您的机器，这可能需要一些时间（通常是几秒钟）。FFTW_MEASURE 是默认计划选项。

...

文档还告诉我们该标志FFTW_ESTIMATE将保留输入。然而，最好的建议是在创建计划后初始化数组。