LAPACK/ScaLAPACK Development

I have recently begun working with Fourier Transforms and have reached a point that confuses me.

(Even though it not a LAPACK question, hopefully, some members might be able to help me.)

I have input the same data to two programs from the NETLIB/FFTPACK library: EZFFTF and CFFTF

The N data inputs are all real, evenly-spaced.

EZFFTF outputs the Fourier coefficients in the form A_k and B_k, where k goes from 1 to N/2
CFFTF outputs the Fourier coefficients in complex form: g_k = gr_k + igi_k, where k goes from 1 to N -1. gr_k and gi_k are the k-th real and imaginary components, respectively. (No subscripts possible in this forum?!)

The--one--reference that I have states that the g_k values can be found from the A_k and B_k values according to the following formula: g_k = (A_k - i B_k)/2, where i indicates the imaginary component.

However, the numbers I am getting are not working out according to this formula.
Either I have used the program(s) incorrectly, or the formula is not correct.

Can anybody here confirm that the formula above is, indeed, the correct way to reconcile output from CFFTF with output from EZFFTF?
Any chance somebody here is running EZFFTF, and could run my numbers, to confirm my results?
I have attached the test data, in case anybody is interested. 30 real data points representing sin(E) values for one period of elliptical motion.

For a description of storage of real and imaginary parts for the FFT of real data, see

http://www.nag.co.uk/numeric/fl/nagdoc_fl24/html/C06/c06intro.html#background12

In EZFFTF, the length, L, of A and B is:

Code: Select all

 
   N odd, (N-1)/2
   N even, N/2


The A(k) and B(k) for k=1,...L give you G(k) = A(K) + i B(K)

The remaining G(K) are implied using the properties

A(N-K+2) = A(K) for K = 2 ,...

and

B(N-k+2) = -B(K) for K = 2,...

You will see that g(n-k+2) is the complex conjugate of g(k) for k = 2,...

Lawrence

A brief follow-up:

The difference in output was caused by a difference in scaling.
EZFFTF scales the A and B coefficients by 1/N.
CFFTF scales the coefficients by 1/(square root(N)).

And the definition of CFFTF did not include a negative sign in the exponent, so the sign of the gi_k was opposite to the sign of B_k.

So ...

Multiplying the A and B coefficients by the square root of N makes the numbers identical.

Multiplying the B coefficients by -1 could be done, but the results are still valid anyhow. Whether the first half of the outputs are used, or the complex conjugates are used--but not both--the results are correct.

Whew!! I have now finished translating EZFFTF from FORTRAN to C++ and JavaScript. A live, working, version is posted here:

Discrete Fourier Transform

If anybody would like to run it through a few tests, to confirm it works, I would welcome any feedback.