-- fft16.adb
with Ada.Numerics.Complex_Types; use Ada.Numerics.Complex_Types;
with Complex_Arrays; use Complex_Arrays;

procedure fft16(Z : in out Complex_Vector) is --  (0..15)
  pragma Suppress(All_Checks);
  W: Complex_Vector(0..15);  -- scratch vector, used many times
  E: Complex_Vector(0..8) :=   -- constants for FFT algorithm
   ((1.0, 0.0), (0.92388, 0.382683),
   (0.707107, 0.707107), (0.382683, 0.92388),
   (0.0, 1.0), (-0.382683, 0.92388),
   (-0.707107, 0.707107), (-0.92388, 0.382683),
   (-1.0, 0.0));
  I, J, K, L, M : Integer;
begin

   M := 8;
   L := 1;
   loop
      K := 0;
      J := L;
      I := 0;
      loop
         loop
            W(I+K) := Z(I) + Z(M+I);
            W(I+J) := E(K) * (Z(I) - Z(M+I));
            I := I+1;
            if I >= J then exit; end if;
         end loop;
         K := J;
         J := K+L;
         if J > M then exit; end if;
      end loop;
      L := L+L;
                -- work back other way without copying
      K := 0;
      J := L;
      I := 0;
      loop
        loop
          Z(I+K) := W(I) + W(M+I);
          Z(I+J) := E(K) * (W(I) - W(M+I));
          I := I+1;
          if I >= J then exit; end if;
        end loop;
        K := J;
        J := K+L;
        if J > M then exit; end if;
      end loop;
      L := L+L;
      if L > M then exit; end if;  -- result is in Z
   end loop;
end fft16;