最近调用FFTW来做reconstruction,遇到一个问题,先贴code:
2 fftw_complex in[N], out[N];
3 fftw_plan p;
4
5 p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
6
7 // 初始化in
8 for (int i = 0; i < N; i++)
9 {
10 in[i][0] = i; // 实数部分
11 in[i][1] = 0; // 虚数部分
12 }
13
14
15 // 执行fftw
16 fftw_execute(p);
17
18 fftw_complex out1[N];
19
20 fftw_plan p1;
21
22 p1 = fftw_plan_dft_1d(N, out, out1, FFTW_BACKWARD, FFTW_ESTIMATE);
23
24 fftw_execute(p1);
25
26 fftw_destroy_plan(p);
27 fftw_destroy_plan(p1);
28
29 double residual = 0;
30
31 for(int i = 0; i < N; i++)
32 {
33 out1[i][0] /= N;
34 out1[i][1] /= N;
35
36 residual += (out1[i][0] - in[i][0])*(out1[i][0] - in[i][0])
37 + (out1[i][1] - in[i][1])*(out1[i][1] - in[i][1]);
38 }
39
40
41 if (residual < 1e-6)
42 cout << "1D FFT passed successful!" << " --- " << "residual : " << residual << endl;
43 else
44 cout << "1D FFT failed!" << " --- " << "residual : " << residual << endl;
可是,却发现,如果如下调用
p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_MEASURE);
则结果会以一定的百分比运行成功,大概是运行十次,成功5次,成功几率50%.
而如果以
p = fftw_plan_dft_1d(N, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
则,结果总是运行成功。
查了一下,官方解释如下:
The flags
argument is usually either FFTW_MEASURE
or FFTW_ESTIMATE
. FFTW_MEASURE
instructs FFTW to run and measure the execution time of several FFTs in order to find the best way to compute the transform of size n
. This process takes some time (usually a few seconds), depending on your machine and on the size of the transform. FFTW_ESTIMATE
, on the contrary, does not run any computation and just builds a reasonable plan that is probably sub-optimal. In short, if your program performs many transforms of the same size and initialization time is not important, use FFTW_MEASURE
; otherwise use the estimate. The data in the in
/out
arrays is overwritten during FFTW_MEASURE
planning, so such planning should be done before the input is initialized by the user.
原因还需要进一步找。感觉还是有问题滴。