Fast Fourier transform

next - skip - up - start

   FFT(X, Y, F, G);                         // F=X and G=Y are OK
where X, Y, F, G are column vectors. X and Y are the real and imaginary input vectors; F and G are the real and imaginary output vectors. The lengths of X and Y must be equal and should be the product of numbers less than about 10 for fast execution.

The formula is 

          n-1
   h[k] = SUM  z[j] exp (-2 pi i jk/n)
          j=0
where z[j] is stored complex and stored in X(j+1) and Y(j+1). Likewise h[k] is complex and stored in F(k+1) and G(k+1). The fast Fourier algorithm takes order n log(n) operations (for good values of n) rather than n**2 that straight evaluation would take.

I use the method of Carl de Boor (1980), Siam J Sci Stat Comput, pp 173-8. The sines and cosines are calculated explicitly. This gives better accuracy, at an expense of being a little slower than is otherwise possible.

Related functions 

   FFTI(F, G, X, Y);                        // X=F and Y=G are OK
   RealFFT(X, F, G);
   RealFFTI(F, G, X);
FFTI is the inverse trasform for FFT. RealFFT is for the case when the input vector is real, that is Y = 0. I assume the length of X, denoted by n, is even. The program sets the lengths of F and G to n/2 + 1. RealFFTI is the inverse of RealFFT.