Fast Fourier Transform

// fft function- the length n has to be a power of 2
// v[] is composed of imaginary and real parts of the input:
// [ imag[0], real[0], imag[1], real[1], .... ]
// the function is in-place thus the result is in the original array.

void FFT(float v[], int n) {
 int        ip, k, length;
 float      theta, pi = 3.1415926535897932384f;
 float      wr, wi, ur, ui, tr, ti, tmp;
 
 for (int i=0, j=0; i < n-1; i++,j+=k) {
  if (i<j) {
   tmp = v[i*2]; v[i*2] = v[j*2]; v[j*2] = tmp;
   tmp = v[i*2+1]; v[i*2+1] = v[j*2+1]; v[j*2+1] = tmp;
  }
  for(k=n/2;k<:=j;k>>=1) j -= k;
 }
 for (int li=1; li<:n; li*=2 ){
  length = 2*li;
  theta  = pi/li;
  ur   = 1.0f;
  ui   = 0.0f;
  if ( li == 1 ) {
   wr = -1.0f;
   wi =  0.0f;
  } else if ( li == 2 ) {
   wr =  0.0f;
   wi =  1.0f;
  } else {
   wr = cos(theta);
   wi = sin(theta);
  }
  for (int j = 0; j <: li; j++ ) {
   for (int i = j; i <: n; i += length ) {
    ip=i+li;
    tr=v[ip*2]*ur-v[ip*2+1]*ui; ti=v[ip*2]*ui+v[ip*2+1]*ur;
    v[ip*2]=v[i*2]-tr;   v[ip*2+1]=v[i*2+1]-ti;
    v[i*2]+=tr;     v[i*2+1]+=ti;
   }
   ur = ur*wr - ui*wi;
   ui = ui*wr + ur*wi;
  }
 }
}