43.2 Radix‑2 FFT Algorithm (Cooley–Tukey)

Radix‑2 FFT requires N = 2^m. It splits the transform into even/odd halves recursively.

43.2.1 Bit‑Reversal Permutation

In‑place FFT needs input in bit‑reversed order:

void bitReverse(float[] x, int N) {
  int bits = Integer.numberOfTrailingZeros(N);
  for (int i = 0; i < N; i++) {
    int rev = Integer.reverse(i) >>> (32 - bits);
    if (rev > i) {
      float tre = x[2*i], tim = x[2*i+1];
      x[2*i] = x[2*rev]; x[2*i+1] = x[2*rev+1];
      x[2*rev] = tre; x[2*rev+1] = tim;
    }
  }
}

43.2.2 Iterative Radix‑2 FFT

void fftRadix2(float[] x, int N) {
  bitReverse(x, N);
  for (int s = 1; s <= Math.log(N)/Math.log(2); s++) {
    int m = 1 << s;
    float wm_re = (float)Math.cos(-2*Math.PI / m);
    float wm_im = (float)Math.sin(-2*Math.PI / m);
    for (int k = 0; k < N; k += m) {
      float w_re = 1f, w_im = 0f;
      for (int j = 0; j < m/2; j++) {
        int t = k + j + m/2;
        int u = k + j;
        float t_re = w_re * x[2*t]   - w_im * x[2*t+1];
        float t_im = w_re * x[2*t+1] + w_im * x[2*t];
        x[2*t]   = x[2*u]   - t_re;
        x[2*t+1] = x[2*u+1] - t_im;
        x[2*u]   += t_re;
        x[2*u+1] += t_im;
        // update w
        float new_w_re = w_re * wm_re - w_im * wm_im;
        float new_w_im = w_re * wm_im + w_im * wm_re;
        w_re = new_w_re; w_im = new_w_im;
      }
    }
  }
}

43.2.3 Inverse FFT

Swap sign in twiddle factors and scale by 1/N:

void ifft(float[] X, int N) {
  // conjugate
  for (int i = 0; i < N; i++) X[2*i+1] = -X[2*i+1];
  fftRadix2(X, N);
  // conjugate and scale
  float scale = 1f / N;
  for (int i = 0; i < N; i++) {
    X[2*i]   *= scale;
    X[2*i+1] *= -scale;
  }
}