重复 FFTW 调用的错误

Question

我遇到了一个无法解决的奇怪问题。我把它作为一个简单的例子来演示这个问题。我在 [0, 2*pi] 之间定义了一个正弦波。我使用 FFTW 进行傅里叶变换。然后我有一个 for 循环，我在其中重复进行逆傅立叶变换。在每次迭代中，我取我的解决方案的平均值并打印结果。我希望每次迭代的平均值保持不变，因为解决方案 y 没有变化。但是，当我选择 N = 256 和 N 的其他偶数时，我注意到平均值会增长，就好像存在数值错误一样。但是，如果我选择 N = 255 或 N = 257，则情况并非如此，我得到了预期的结果（每次迭代平均 = 0.0）。

代码：

#include <stdio.h>
#include <stdlib.h>
#include <fftw3.h>
#include <math.h>

int main(void)
{
  int N = 256;
  double dx = 2.0 * M_PI / (double)N, dt = 1.0e-3;
  double *x, *y;

  x = (double *) malloc (sizeof (double) * N); 
  y = (double *) malloc (sizeof (double) * N);

  // initial conditions
  for (int i = 0; i < N; i++) {
    x[i] = (double)i * dx;
    y[i] = sin(x[i]);
  }

  fftw_complex yhat[N/2 + 1];
  fftw_plan fftwplan, fftwplan2;    

  // forward plan
  fftwplan = fftw_plan_dft_r2c_1d(N, y, yhat, FFTW_ESTIMATE);
  fftw_execute(fftwplan);

  // set N/2th mode to zero if N is even
  if (N % 2 < 1.0e-13) {
    yhat[N/2][0] = 0.0;
    yhat[N/2][1] = 0.0;
  }

  // backward plan
  fftwplan2 = fftw_plan_dft_c2r_1d(N, yhat, y, FFTW_ESTIMATE);

  for (int i = 0; i < 50; i++) {
    // yhat to y  
    fftw_execute(fftwplan2);

    // rescale
    for (int j = 0; j < N; j++) {
      y[j] = y[j] / (double)N;
    }

    double avg = 0.0;
    for (int j = 0; j < N; j++) {
      avg += y[j];
    }
    printf("%.15f\n", avg/N);
  }

  fftw_destroy_plan(fftwplan);
  fftw_destroy_plan(fftwplan2);
  void fftw_cleanup(void);
  free(x);
  free(y);
  return 0;
}

N = 256 的输出：

0.000000000000000
0.000000000000000
0.000000000000000
-0.000000000000000
0.000000000000000
0.000000000000022
-0.000000000000007
-0.000000000000039
0.000000000000161
-0.000000000000314
0.000000000000369
0.000000000004775
-0.000000000007390
-0.000000000079126
-0.000000000009457
-0.000000000462023
0.000000000900855
-0.000000000196451
0.000000000931323
-0.000000009895302
0.000000039348379
0.000000133179128
0.000000260770321
-0.000003233551979
0.000008285045624
-0.000016331672668
0.000067450106144
-0.000166893005371
0.001059055328369
-0.002521514892578
0.005493164062500
-0.029907226562500
0.093383789062500
-0.339111328125000
1.208251953125000
-3.937500000000000
13.654296875000000
-43.812500000000000
161.109375000000000
-479.250000000000000
1785.500000000000000
-5369.000000000000000
19376.000000000000000
-66372.000000000000000
221104.000000000000000
-753792.000000000000000
2387712.000000000000000
-8603776.000000000000000
29706240.000000000000000
-96833536.000000000000000

有什么想法吗？

Answer 1

libfftw 有修改输入的恶习。如果要重复逆变换，请备份yhat。

OTOH，这是有悖常理的，但如果您不希望它产生不同的结果，为什么要重复相同的操作？（尽管如此）

Answer 2

如评论中所示："if you want to keep the input data unchanged, use the FFTW_PRESERVE_INPUT flag. Per http://www.fftw.org/doc/Planner-Flags.html"

例如：

// backward plan
fftwplan2 = fftw_plan_dft_c2r_1d(N, yhat, y, FFTW_ESTIMATE | FFTW_PRESERVE_INPUT);

重复 FFTW 调用的错误

Errors with repeated FFTW calls

fft

transform

fftw