20061102

The Impossible FFTW Bug

Ever had a computer bug that you spent weeks working on it? You might try every trick you know to solve it: debugging messages, profilers, Google, professors, coworkers, fellow students. You might try to explain the code to someone who doesn't even program just so you have to describe the logic in simple terms. Nothing seems to fix the problem.

I had this bug in my code for about two weeks. I had to write a program to compute the 2 dimensional FFT of an image. That's it. I'll just use FFTW, the most popular open source library for computing FFTs.

The Wrong Solution for computing FFTs with FFTW:

int j;
fftw_plan myplan;
fftw_complex in[x*y], out[x*y];

for (j = 0; j < x*y; j++) {
    in[j][0] = r[j];
    in[j][1] = i[j];
} // End For

myplan = (dir == 1)? fftw_plan_dft_2d(x, y, in, out, FFTW_FORWARD,  FFTW_MEASURE) :
                     fftw_plan_dft_2d(x, y, in, out, FFTW_BACKWARD, FFTW_MEASURE) ;

fftw_execute(myplan);

In these few lines of code lies a problem that I've spent weeks trying to find. There are no syntax errors, nor are there logic errors. Every time I ran this code, the result would be all zeros. If I ran it multiple times, I got the correct results for every execution after the first. I dropped the data into a prepared array, made a plan (as defined by the FFTW C library), then executed that library. Should work, right? Wrong.

The Correct Solution for computing FFTs with FFTW:

int j;
fftw_plan myplan;
fftw_complex in[x*y], out[x*y];

myplan = (dir == 1)? fftw_plan_dft_2d(x, y, in, out, FFTW_FORWARD,  FFTW_MEASURE) :
                     fftw_plan_dft_2d(x, y, in, out, FFTW_BACKWARD, FFTW_MEASURE) ;

for (j = 0; j < x*y; j++) {
    in[j][0] = r[j];
    in[j][1] = i[j];
} // End For

fftw_execute(myplan);

I found buried in the FFTW FAQ that I'm suppose to make a plan before giving it the data. I switched two sets of lines around and everything magically works for reasons I fail to understand. Let this be a lesson to anyone having to use FFTW.

2 comments:

Anonymous said...

That's because it measures the speed using different pieces of code. For that, it has to write to your reserved memory.

Yungui Xu said...

hi , thank you for your post here, and i have the same problem. with you infomation, i get i.

yungui.xu@gmail.com