-- ifft2048.adb   inverse fft
with Ada.Numerics.Complex_Types; use Ada.Numerics.Complex_Types;
with Complex_Arrays; use Complex_Arrays;

procedure ifft2048(A : in out Complex_Vector) is   -- complex vector A(0..2047)
  pragma Suppress(All_Checks);
  E : Complex_Vector (0..10) := -- constants for inverse fft
   ((-1.0, 0.0), (0.0, -1.0),
   (0.707107, -0.707107), (0.92388, -0.382683),
   (0.980785, -0.19509), (0.995185, -0.0980171),
   (0.998795, -0.0490677), (0.999699, -0.0245412),
   (0.999925, -0.0122715), (0.999981, -0.00613588),
   (0.999995, -0.00306796));
  N : constant Integer := 2048;
  FN : constant Float := float(N);
  TMP : Complex;
  WX : Complex;
  W : Complex;
  I : Integer :=0;
  JS : Integer :=1;
  M : Integer;
  IX : Integer;
  ISP : Integer;
  MMAX : Integer :=1;
begin

  for IX in 1..N-1 loop  -- reorder data
    if JS > IX then
      TMP := A(JS-1);
      A(JS-1) := A(IX-1);
      A(IX-1) := TMP;
    end if;
    M := N / 2;
    while M < JS and M > 0 loop
      JS := JS - M;
      M := M / 2;
    end loop;
    JS := JS + M;
  end loop;
  while MMAX < N loop  -- compute transform
    ISP := MMAX + MMAX;
    WX := E(I);
    I := I+1;
    W := (1.0, 0.0);
    for M in 0..MMAX-1 loop
      IX := M;
      while IX+MMAX < N loop
        JS := IX + MMAX;
        TMP := W * A(JS);
        A(JS) := A(IX) - TMP;  -- basic butterfly
        A(IX) := A(IX) + TMP;
        IX := IX + ISP;
      end loop;
      W := W * WX;
    end loop;
    MMAX := ISP;
  end loop;
  -- only divide by N on inverse transform
  for I in 0..N-1 loop
    A(I) := A(I) / FN;
  end loop;
end ifft2048;