FFTW: Inverse of forward fft not equal to original function

Question

I'm trying to use FFTW to compute fast summations, and I've run into an issue:

int numFreq3 = numFreq*numFreq*numFreq;

FFTW_complex* dummy_sq_fft = (FFTW_complex*)FFTW_malloc( sizeof(FFTW_complex)*numFreq3 );

FFTW_complex* dummy_sq = (FFTW_complex*)FFTW_malloc( sizeof(FFTW_complex)*numFreq3 );

FFTW_complex* orig = (FFTW_complex*)FFTW_malloc( sizeof(FFTW_complex)*numFreq3 );

FFTW_plan dummyPlan = FFTW_plan_dft_3d( numFreq, numFreq, numFreq,
                  orig, dummy_sq_fft,
                  FFTW_FORWARD, FFTW_MEASURE );

FFTW_plan dummyInvPlan = FFTW_plan_dft_3d( numFreq, numFreq, numFreq,
                      dummy_sq_fft, dummy_sq,
                      FFTW_BACKWARD, FFTW_MEASURE );

for(int i= 0; i < numFreq3; i++) {
  orig[ i ][ 0 ] = sparseProfile02[ 0 ][ i ][ 0 ];
  //img. part == 0
  orig[ i ][ 1 ]  = sparseProfile02[ 0 ][ i ] [ 1 ];
}

FFTW_execute(dummyPlan);
FFTW_execute(dummyInvPlan);

int count = 0; 
for(int i=0; i<numFreq3; i++) {
  double factor = dummy_sq[ i ][ 0 ]/sparseProfile02[ 0 ][ i ][ 0 ];

  if(factor < 0) {
    count++;
  }
}

std::cout<<"Count "<<count<<"\n";

FFTW_free(dummy_sq_fft);
FFTW_free(dummy_sq);
FFTW_destroy_plan(dummyPlan);
FFTW_destroy_plan(dummyInvPlan);

(Here sparseProfile02[0] is of type FFTW_complex*, and contains only positive real data.)

Since we have dummy_sq = IFFT(FFT(sparseProfile02[ 0 ])), we must have dummy_sq = n^3*sparseProfile02. But this is true only some of the time; in fact, values on the dummy_sq grid are negative whenever corresponding values on the sparseProfile02 grid are zero (but not vice-versa). Does anyone know why this is happening?

Answer 1

The FFTs (forward and inverse) have rounding error, and I think this is what's biting you. In general, you shouldn't expect a zero to stay exactly zero through your process (although it could be zero for trivial test cases). In your test loop, is

fabs(dummy_sq[i][0] - numFreq*numFreq*numFreq*sparseProfile02[0][i][0])

large relative to the magnitude of your data?

As a really simple (pathological) example with just a 1D FFT of size 2 with real values:

ifft(fft([1e20, 1.0])) != [2e20, 2.0]

The 1.0 is lost in the 1e20 with a double precision FFT.

It's also possible you're getting some NaNs for factor when you divide by a zero sample in sparseProfile02.