Fourier Transform - Help Request

[countzero]

I want to write a c function which implements fourier transform in 2d in a quick manner. I wrote my own version by looking at the formula and it seems to work ok. However I was not very pleased with the speed so I wrote a fast fourier function from the book numerical recipes in c. However, when I check the results don't match. There is a -1 coefficient difference in the imaginary set with my inverse fourier and numerical recipes in c's fast fourier version.  I was wondering how can I go about and check which one is correct ? does anyone have an idea ?

numerical recipes in c
wikipedia fourier transform(where I implemented the formula from)
another page with fft c code(probably buggy though)

my humble code for those willing to debug : (which assumes all data is real)

Code: [Select]



double m2pi = -2.0*pi;
for (u=0; u<M; u++)
for (v=0; v<N; v++) {
temp = 0, itemp = 0;
for (x=0; x<M; x++)
for (y=0; y<N; y++) {
temp += A[x][y].re*cos(m2pi*u*x/M+m2pi*v*y/N);
itemp += A[x][y].re*sin(m2pi*u*x/M+m2pi*v*y/N);
}
B[u][v].re = temp;
B[u][v].im = itemp;
}


and the inverse

Code: [Select]


double m2pi = 2.0*pi;
for (u=0; u for (v=0; v temp = 0, itemp = 0;
for (x=0; x for (y=0; y constant = m2pi*u*x/M+m2pi*v*y/N;
cos_coef = cos(constant);
sin_coef = sin(constant);
temp += B[x][y].re*cos_coef - B[x][y].im*sin_coef;
itemp += B[x][y].im*cos_coef + B[x][y].re*sin_coef;
}
C[u][v].re = temp/(M*N);
C[u][v].im = itemp/(M*N);
}



One for Piru & Karlos I reckon.

[Oliver]

I don't have time to debug your code right now, but I suggest you check all your output with a known maths package like MATLAB.  I believe MATLAB has 2D FFT transforms available.  It's also very easy to use.

P.S. Is this for image analysis?

P.P.S. I'm pretty sure there's a lot of literature available on this topic.  You could have a look at some journals if general web pages are not sufficient.

[Cymric]

Numerical Recipes doesn't use zero-offset arrays, but against ALL GOOD PROGRAMMING PRACTICE has them start at 1. This has been an endless source of confusion and screws up efficient coding because you cannot use reverse loops (for (i=MAX; i--;)) easily.

What I'm trying to say is: did you take this into account when testing?

[countzero]

yes, it was a nuisiance but I believe I fixed it.

now, when I use the numerical recipes code, it just works both inverse and forward. I believe I'm doing something wrong in my code, cause I get the negative of the imaginary part. However I can still get the original sequence with the inverse transform (my code). my code vs num. rep. just doesn't work because of the difference of the imaginary parts.

no, it's not for image processing. It was for somekind of contest, which I already failed to submit in time.

[Karlos]

Last time I looked, jpeglib had a utility file full of fft functions...

[Cyberus]

Oliver wrote:
I don't have time to debug your code right now, but I suggest you check all your output with a known maths package like MATLAB.  I believe MATLAB has 2D FFT transforms available.  It's also very easy to use.

P.S. Is this for image analysis?

P.P.S. I'm pretty sure there's a lot of literature available on this topic.  You could have a look at some journals if general web pages are not sufficient.

When I was doing my Imaging course, we used Matlab to do FFTs. The great thing being you can easily get decent graphical output as well.

I'm guessing you could you also use Fortran?

[countzero]

ok, I just setup matlab to test results and it turns out that my code was correct after all and the num reps was wrong. probably I made a mistake when I was fiddling with the offsets. Now I'm trying to understand what actually num. rep. code does and debug it. I wonder if it would be easier if I set about to writing my own fft code.

[Cymric]

Because I couldn't sleep from the heat in these parts (I think it was 23 or 24 degrees C outside, I worked a bit on your code, and concluded it contained no errors. But you already saw that yourself too.

I looked up the code for the 2D-FFT, and it looks like there are quite a lot of assumptions being made which could all too easily be broken if you don't follow NR's coding standard exactly. You need to grab the PDF for subchapter 12.4 to get an idea. But I trust you already figured that out too. Good luck with the code study, in any case.

[Hyperspeed]

[Speelgoedmannetje]

Cymric wrote:
Numerical Recipes doesn't use zero-offset arrays, but against ALL GOOD PROGRAMMING PRACTICE has them start at 1. This has been an endless source of confusion and screws up efficient coding because you cannot use reverse loops (for (i=MAX; i--;)) easily.

How can I agree more. Zero is an offset, not a number.


o.t.
I'll check your code later.