Errors with repeated FFTW calls

2
votes

I'm having a strange issue that I can't resolve. I made this as a simple example that demonstrates the problem. I have a sine wave defined between [0, 2*pi]. I take the Fourier transform using FFTW. Then I have a for loop where I repeatedly take the inverse Fourier transform. In each iteration, I take the average of my solution and print the results. I expect that the average stays the same with each iteration because there is no change to solution, y. However, when I pick N = 256 and other even values of N, I note that the average grows as if there are numerical errors. However, if I choose, say, N = 255 or N = 257, this is not the case and I get what is expect (avg = 0.0 for each iteration).

Code:

#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;
}

Output for 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

Any ideas?

2
transformfftfftw

2 Answers

1
votes

libfftw has the odious habit of modifying its inputs. Back up yhat if you want to do repeated inverse transforms.

OTOH, it's perverse, but why are you repeating the same operation if you don't expect it give different results? (Despite this being the case)

1
votes

As indicated in comments: "if you want to keep the input data unchanged, use the FFTW_PRESERVE_INPUT flag. Per http://www.fftw.org/doc/Planner-Flags.html"

For example:

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