使用 FFTW_MEASURE 时 FFT 输出为空白,但使用 FFTW_ESTIMATE 时可以正常工作

FFT output is blank when using FFTW_MEASURE, but works fine with FFTW_ESTIMATE

我在尝试使用 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

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

documentation of the planner flags of FFTW 详细说明了正在发生的事情。我很确定您已经猜到 FFTW_ESTIMATE 保留输入数组并 FTTW_MEASURE 修改它的原因。

Important: the planner overwrites the input array during planning unless a saved plan (see Wisdom) is available for that problem, so you should initialize your input data after creating the plan.*** The only exceptions to this are the FFTW_ESTIMATE and FFTW_WISDOM_ONLY flags, as mentioned below.

...

  • FFTW_ESTIMATE specifies that, instead of actual measurements of different algorithms, a simple heuristic is used to pick a (probably sub-optimal) plan quickly. With this flag, the input/output arrays are not overwritten during planning.
  • FFTW_MEASURE tells FFTW to find an optimized plan by actually computing several FFTs and measuring their execution time. Depending on your machine, this can take some time (often a few seconds). FFTW_MEASURE is the default planning option.

...

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