重复 FFTW 调用的错误答案

【问题标题】：Errors with repeated FFTW calls重复 FFTW 调用的错误
【发布时间】：2015-04-21 01:06:05
【问题描述】：

我遇到了一个无法解决的奇怪问题。我把它作为一个简单的例子来演示这个问题。我在 [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

有什么想法吗？

【问题讨论】：

标签： transform fft fftw

【解决方案1】：

libfftw 有修改其输入的可恶习惯。如果您想进行重复的逆变换，请备份 yhat。

OTOH，这很反常，但是如果您不期望它会产生不同的结果，为什么还要重复相同的操作呢？（尽管如此）

【讨论】：

好吧，在实际代码中，结果应该会改变（每次迭代都会修改 yhat），但我遇到了同样的问题，偶数 N 的解决方案会增长，而奇数 N 会给出正确的结果。在这里，我想让示例尽可能简单。但是，我猜 yhat 每次迭代都会改变真正的问题，这意味着我不能使用单个备份。
如果要保持输入数据不变，请使用FFTW_PRESERVE_INPUT 标志。每fftw.org/doc/Planner-Flags.html
我修改了c2r 计划以包含FFTW_ESTIMATE | FFTW_PRESERVE_INPUT，它似乎在这里和我的实际问题中都成功了。谢谢！
感谢 doco，我找不到。 C2r 默认是破坏性的。 FFTW 现在在我的书中不那么可憎了。
你能把这个写成@symmetric 的答案并接受吗？

【解决方案2】：

如 cmets 所示：“如果要保持输入数据不变，请使用FFTW_PRESERVE_INPUT 标志。每http://www.fftw.org/doc/Planner-Flags.html”

例如：

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

【讨论】：