我在尝试使用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
修改输入数组的原因。
...
...
该文档还告诉我们,标记FFTW_ESTIMATE
将保留输入。但是,最好的建议是一旦创建计划就初始化阵列。