-- Generic_Complex_Fast_Fourier_Transform

with ARRAY_EXCEPTIONS;

generic -- package YOUR_NAME is new GENERIC_COMPLEX_FAST_FOURIER_TRANSFORM
        --     ( REAL, COMPLEX, COMPLEX_VECTOR, COMPLEX_MATRIX);
        --        do the instantiation where the instantiation of
        --        GENERIC_COMPLEX_TYPES is visible

  type REAL is digits <> ;
  type COMPLEX is private;
  type COMPLEX_VECTOR is array(INTEGER range <>) of COMPLEX;
  type COMPLEX_MATRIX is array(INTEGER range <>,INTEGER range <>) of COMPLEX;
  with function COMPOSE_FROM_CARTESIAN(L,R:REAL) return COMPLEX is <>;
  with function RE(L:COMPLEX) return REAL is <>;
  with function IM(L:COMPLEX) return REAL is <>;
  with function "+"(L,R:COMPLEX) return COMPLEX is <>;
  with function "-"(L,R:COMPLEX) return COMPLEX is <>;
  with function "*"(L,R:COMPLEX) return COMPLEX is <>;
  with function "/"(L,R:COMPLEX) return COMPLEX is <>;
  with function "*"(L:COMPLEX;R:REAL) return COMPLEX is <>;
  with function "/"(L:COMPLEX;R:REAL) return COMPLEX is <>;

package GENERIC_COMPLEX_FAST_FOURIER_TRANSFORM is

  procedure FFT ( A : in out COMPLEX_VECTOR ) ;
  function  FFT ( A : COMPLEX_VECTOR ) return COMPLEX_VECTOR ;

  --      PURPOSE : COMPUTE DISCRETE FAST FOURIER TRANSFORM OF
  --                THE POINTS IN THE COMPLEX VECTOR
  --
  --      INPUT : THE COMPLEX VECTOR A ( e.g. TIME DOMAIN )
  --
  --      OUTPUT : THE TRANSFORM OF VECTOR A ( e.g. FREQUENCY DOMAIN )
  --               IS RETURNED IN PLACE FOR THE PROCEDURE,
  --               IS RETURNED AS THE VALUE OF THE FUNCTION
  --               SCALED BY DIVIDING BY NUMBER OF POINTS
  --
  --      METHOD : N LOG N FFT ALGORITHM

  procedure INVERSE_FFT ( A : in out COMPLEX_VECTOR ) ;
  function  INVERSE_FFT ( A : COMPLEX_VECTOR ) return COMPLEX_VECTOR ;

  --      PURPOSE : COMPUTE INVERSE DISCRETE FAST FOURIER TRANSFORM OF
  --                THE POINTS IN THE COMPLEX VECTOR
  --
  --      INPUT : THE COMPLEX VECTOR A ( e.g. FREQUENCY DOMAIN )
  --
  --      OUTPUT : THE INVERSE TRANSFORM OF VECTOR A ( e.g. TIME DOMAIN )
  --               IS RETURNED IN PLACE FOR THE PROCEDURE,
  --               IS RETURNED AS THE VALUE OF THE FUNCTION
  --               NO SCALING
  --
  --      METHOD : N LOG N FFT ALGORITHM

  procedure FEJER ( A : in out COMPLEX_VECTOR ) ;
  function  FEJER ( A : COMPLEX_VECTOR ) return COMPLEX_VECTOR ;

  --      PURPOSE : COMPUTE FEJER MODIFIED FOURIER TRANSFORM OF
  --                THE POINTS IN THE COMPLEX VECTOR
  --
  --      INPUT : THE COMPLEX VECTOR A ( e.g. TIME DOMAIN )
  --
  --      OUTPUT : THE FEJER TRANSFORM OF VECTOR A ( e.g. FREQUENCY DOMAIN )
  --               IS RETURNED IN PLACE FOR THE PROCEDURE,
  --               IS RETURNED AS THE VALUE OF THE FUNCTION
  --
  --      METHOD : N LOG N FFT ALGORITHM
  --               APPLY WEIGHTING OF (N-I)/N TO I TH FREQUENCY
  --               WHERE N IS THE LENGHT OF THE VECTOR

end GENERIC_COMPLEX_FAST_FOURIER_TRANSFORM ;

with GENERIC_ELEMENTARY_FUNCTIONS ;
with MATHEMATICAL_CONSTANTS ;

package body GENERIC_COMPLEX_FAST_FOURIER_TRANSFORM is

  --      PURPOSE : TO PERFORM THE BASIC FFT AND INVERSE FFT COMPUTATION
  --                THIS IS AN INTERNAL PROCEDURE USED BY THE ABOVE.

  package ELEMENTARY_FUNCTIONS is new GENERIC_ELEMENTARY_FUNCTIONS(REAL);



  procedure FFT_COMPUTATION ( A : in out COMPLEX_VECTOR ;
                              TRANS : REAL ) is

  --      PURPOSE : TO PERFORM THE BASIC FFT AND INVERSE FFT COMPUTATION
  --
  --
  --    WRITTEN BY : JON SQUIRE , 24 MAY 1983
  --    JSS REVISED to ISO standard packages 26 AUGUST 1990
  --    JSS Cleanup 19 October 1990
  --    JSS 15 June 1991, Make number of points not required to be a power of 2
  --    JSS 9 Sept 1996, Clean up generic formal parameters
  --
  --                                        LOCAL SPECIFICATIONS

    PI : REAL := MATHEMATICAL_CONSTANTS.PI;
    TMP : COMPLEX ;
    WX : COMPLEX ;
    W  : COMPLEX ;
    JS : INTEGER ;
    M : INTEGER ;
    IX : INTEGER ;
    ISP : INTEGER ;
    MMAX : INTEGER ;
    PH1 : REAL ;
    N : INTEGER := A'LENGTH ;
    OFFSET : INTEGER := A'FIRST - 1 ;
  begin
    JS := 1 ;
    for IX in 1 .. N loop  -- REORDER DATA
      if JS > IX then
        TMP := A ( JS + OFFSET ) ;
        A ( JS + OFFSET ) := A ( IX + OFFSET ) ;
        A ( IX + OFFSET ) := 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 ;
    MMAX := 1 ;
    while MMAX < N loop  -- COMPUTE TRANSFORM
      ISP := MMAX + MMAX ;
      PH1 := PI * TRANS / REAL ( MMAX ) ;
      WX := COMPOSE_FROM_CARTESIAN ( ELEMENTARY_FUNCTIONS.COS( PH1 ) ,
                                     ELEMENTARY_FUNCTIONS.SIN( PH1 )) ;
      W := COMPOSE_FROM_CARTESIAN (1.0, 0.0) ;
      for M in 1 .. MMAX loop
        IX := M ;
        while IX <= N loop
          JS := IX + MMAX ;
          TMP := W * A ( JS + OFFSET ) ;
          A ( JS + OFFSET ) := A ( IX + OFFSET ) - TMP ;  -- BASIC BUTTERFLY
          A ( IX + OFFSET ) := A ( IX + OFFSET ) + TMP ;
          IX := IX + ISP ;
        end loop ;
        W := W * WX ;
      end loop ;
      MMAX := ISP ;
    end loop ;
    if TRANS < 0.0 then -- ONLY DIVIDE BY N ON INVERSE TRANSFORM
      for IX in 1 .. N loop
        A ( IX + OFFSET ) := A ( IX + OFFSET ) / REAL(N);
      end loop ;
    end if ;
  end FFT_COMPUTATION ;

  procedure FFT ( A : in out COMPLEX_VECTOR ) is

    --      PURPOSE : COMPUTE DISCRETE FAST FOURIER TRANSFORM OF
    --                THE POINTS IN THE COMPLEX VECTOR
    --
    --      INPUT : THE COMPLEX VECTOR A ( e.g. TIME DOMAIN )
    --
    --      OUTPUT : THE TRANSFORM OF VECTOR A ( e.g. FREQUENCY DOMAIN )
    --               IS RETURNED IN PLACE
    --
    --      METHOD : N LOG N FFT ALGORITHM
    N : INTEGER := A'LENGTH ;
    NP : INTEGER := 2;
  begin
    if N = 0 then
      return;
    end if ;
    while NP < N loop
      NP := NP + NP ;
    end loop;
    if NP = N then -- 'LENGTH is a power of 2
      FFT_COMPUTATION ( A , 1.0 ) ;
    else
      declare
        B : COMPLEX_VECTOR(0..NP-1); -- make 'LENGTH a power of 2
      begin
        B := A & COMPLEX_VECTOR'((1..NP-N => COMPOSE_FROM_CARTESIAN(0.0,0.0)));
        FFT_COMPUTATION ( B , 1.0 ) ;
        A := B(0..N-1);
      end;
    end if;
  end FFT ;

  function  FFT ( A : COMPLEX_VECTOR ) return COMPLEX_VECTOR is
  N : INTEGER := A'LENGTH ;
    NP : INTEGER := 2;
  begin
    if N = 0 then
      return A;
    end if ;
    while NP < N loop
      NP := NP + NP ;
    end loop;
    declare
      B : COMPLEX_VECTOR(0..NP-1); -- make 'LENGTH a power of 2
    begin
      B := A & COMPLEX_VECTOR'((1..NP-N => COMPOSE_FROM_CARTESIAN(0.0,0.0)));
      FFT_COMPUTATION ( B , 1.0 ) ;
      return B(0..N-1) ;
    end;
  end FFT ;

  procedure INVERSE_FFT ( A : in out COMPLEX_VECTOR ) is

    --      PURPOSE : COMPUTE INVERSE DISCRETE FAST FOURIER TRANSFORM OF
    --                THE POINTS IN THE COMPLEX VECTOR
    --
    --      INPUT : THE COMPLEX VECTOR A ( e.g. FREQUENCY DOMAIN )
    --
    --      OUTPUT : THE INVERSE TRANSFORM OF VECTOR A ( e.g. TIME DOMAIN )
    --               IS RETURNED IN PLACE
    --
    --      METHOD : N LOG N FFT ALGORITHM
    N : INTEGER := A'LENGTH ;
    NP : INTEGER := 2;
  begin
    if N = 0 then
      return;
    end if ;
    while NP < N loop
      NP := NP + NP ;
    end loop;
    if NP = N then -- 'LENGTH is a power of 2
      FFT_COMPUTATION ( A , -1.0 ) ;
    else
      declare
        B : COMPLEX_VECTOR(0..NP-1); -- make 'LENGTH a power of 2
      begin
        B := A & COMPLEX_VECTOR'((1..NP-N => COMPOSE_FROM_CARTESIAN(0.0,0.0)));
        FFT_COMPUTATION ( B , -1.0 ) ;
        A := B(0..N-1);
      end;
    end if;
  end INVERSE_FFT ;

  function  INVERSE_FFT ( A : COMPLEX_VECTOR ) return COMPLEX_VECTOR is
    N : INTEGER := A'LENGTH ;
    NP : INTEGER := 2;
  begin
    if N = 0 then
      return A;
    end if ;
    while NP < N loop
      NP := NP + NP ;
    end loop;
    declare
      B : COMPLEX_VECTOR(0..NP-1); -- make 'LENGTH a power of 2
    begin
      B := A & COMPLEX_VECTOR'((1..NP-N => COMPOSE_FROM_CARTESIAN(0.0,0.0)));
      FFT_COMPUTATION ( B , -1.0 ) ;
      return B(0..N-1) ;
    end;
  end INVERSE_FFT ;

  procedure FEJER ( A : in out COMPLEX_VECTOR ) is
    --      PURPOSE : COMPUTE FEJER MODIFIED FOURIER TRANSFORM OF
    --                THE POINTS IN THE COMPLEX VECTOR
    --
    --      INPUT : THE COMPLEX VECTOR A ( e.g. TIME DOMAIN )
    --              A'LENGTH MUST BE A POWER OF 2
    --
    --      OUTPUT : THE FEJER TRANSFORM OF VECTOR A ( e.g. FREQUENCY DOMAIN )
    --               IS RETURNED IN PLACE FOR THE PROCEDURE,
    --               IS RETURNED AS THE VALUE OF THE FUNCTION
    --
    --      METHOD : N LOG N FFT ALGORITHM
    --               APPLY WEIGHTING OF (N-I)/N TO I TH FREQUENCY
    --               WHERE N IS THE LENGHT OF THE VECTOR
  begin
    FFT_COMPUTATION ( A , 1.0 ) ;
    for I in INTEGER(0)..A'LENGTH - 1 loop
      A(I+A'FIRST) := A(I+A'FIRST) * REAL(A'LENGTH-I) / REAL(A'LENGTH) ;
    end loop;
  end FEJER ;

  function  FEJER ( A : COMPLEX_VECTOR ) return COMPLEX_VECTOR is
    B : COMPLEX_VECTOR ( A'RANGE ) := A;
  begin
    FFT_COMPUTATION ( B , 1.0 ) ;
    for I in INTEGER(0)..B'LENGTH - 1 loop
      B(I+B'FIRST) := B(I+B'FIRST) * REAL(B'LENGTH-I) / REAL(B'LENGTH) ;
    end loop;
    return B ;
  end FEJER ;

end GENERIC_COMPLEX_FAST_FOURIER_TRANSFORM ;