Hi,

I am comparing four1() with some other implementations,
including my own. I found that the output from four1() is
twice the value of those from other implementations.

four1() is compared with following FFTs:
twofft(), realft(), Matlab fft(), and my own fft().

The input is a real vector of length 4, [0 1 2 3]. The vector is repackaged according to the input requirement of
each implementation. In the case of twofft(), I simply use two
identical vectors as input so that the output will also be the
same. Here is the result:

four1():

12
-4.0-4.0i
-4.0
-4.0+4.0i

twofft();
6.0 6.0
-2.0-2.0i -2.0-2.0i
-2.0 -2.0
-2.0+2.0i -2.0+2.0i

realft():
6.0
-2.0-2.0i
-2.0

Matlab fft() and my fft() gave the same results as twofft().

I can not understand the cause of value doubling from four1().
1. All others gave the same output. They have to be correct.
2. twofft() and realft() make use of four1(). Since these two
routines produce the same results as others, they must be
correct as well. Therefore, four1() has to be correct.

So, what is wrong here? Do I miss any thing? Please shed some
light on this. Thanks.

Steve

Hi, Steve!

Alas, I am not able to reproduce your results.

Here is my complete program (omiting only the guts of four1) :


#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr

void four1(float data[], unsigned long nn, int isign)
{
unsigned long n,mmax,m,j,istep,i;
//... rest of four1...

}
#undef SWAP

int main(int argc, char* argv[])
{

float data[8] = {0.,0.,1.,0.,2.,0.,3.,0.} ;

int i;
int isign=1;
unsigned long nn=4;

four1(data-1,nn,isign);
// note trick above to use the zero-based array "data"
for (i=0;i<8;i++) printf("%d %f\n",i,data[i]);

return 0;
}


This gives as output:


0 6.000000
1 0.000000
2 -2.000000
3 -2.000000
4 -2.000000
5 0.000000
6 -2.000000
7 2.000000


which is correct and (as you note) 1/2 of what you get. I think you've been bitten by a bug that is not in four1.

Cheers,
Bill P.

Hi Bill,

Thanks for your reply. It turns out that a mistake was made
in my test code, as you rightly pointed out in your response.

I simply assign the vector elements with the value of the loop
index. However, in test code for four1(), the index increments
by 2, whcih effectively created a real-valued vector twice the
magnitude of what I intended to. The loop index in test codes
for other implementation all increases by 1 at a time. That's why
their output was correct.

Mystery solved. Thank you.

Steve

Hi Bill,

I am having somewhat different problem; The spectrum seems to give -ve frequencies first and +ve later. I ran your example and get what you get. BUT, when I run it on MATLAB, this is what I get. I am reproducing exactly:


>> x = [0 1 2 3 ]'

x =

0
1
2
3

>> y = fft(x)

y =

6.0000
-2.0000 + 2.0000i
-2.0000
-2.0000 - 2.0000i

>>

Which is different from what you got. Can you please explain the difference. (FYI, MATLAB also claims to list +ve freq. before listing -ve)

Thanks

Here I am also stuck *badly* with four1(). I took care of the -1 shift in indexing and the isign = -1 vs isign = +1. I still get different (similar, but...) results compared to MatLab. Since, I know, my ML routine works (I do correlations of images with ML), I am using ML as a positive control. Also, ifft(fft_data) in ML gives the original data, while four1(fft_data,nn,1)/2*nn does not, which makes me "trust" ML more.

Examples:
//C++

float data[16] = {1,2,2,3,3,4,4,5,5,6,6,7,7,8,8,9};
float * fft_data = four1(data-1, 8, -1);

output:

36.000000
44.000000
-3.414214
1.414214
-8.000000
0.000000
-1.414214
-0.585786
-16.000000
-16.000000
-0.585786
-1.414214
0.000000
-8.000000
1.414214
-3.414214

MatLab:

a = [ 1+2*i 2+3*i 3+4*i 4+5*i 5+6*i 6+7*i 7+8*i 8+9*i];
b = fft(a);

ouptut:
36.0000 +44.0000i
-13.6569 + 5.6569i
-8.0000
-5.6569 - 2.3431i
-4.0000 - 4.0000i
-2.3431 - 5.6569i
0 - 8.0000i
5.6569 -13.6569i

From the descriptive (i.e. lame) standpoint, I only see that:
Zero frequencies are identical,
in ML, 1/(2*delta) is 4 times smaller
the others are 4 times bigger (abs values, I mean).

I got the exact same result with {1,2,2,3,3,4,4,5};

This is driving me nuts.

I may start looking into FFTW instead, although it's going to be more painful than copying a text from a .PDF and pasting it.

Andrey.

Hi Andrey,

I don't think this thread is being visited that frequently (at least by those people who can say something useful). If you need some C code for FFT, I can suggest this Application Note from Freescale Semiconductor (formerly Motorola).

http://www.freescale.com/files/32bit/doc/app_note/AN2114.pdf

It has fixed point FFT code that should perform better than four1()''s floating point performance.