mathpack.vhd

fast fourior transform

-- History:
-- 	Version 0.1	Jose A. Torres	4/23/93	First draft
-- 	Version 0.2	Jose A. Torres	5/28/93	Fixed potentially illegal code
-------------------------------------------------------------
Library CUSTOM;

Package body MATH_REAL is

    -- 
    -- some constants for use in the package body only
    --
    constant  Q_PI :   real := MATH_PI/4.0;
    constant  HALF_PI :   real := MATH_PI/2.0;
    constant  TWO_PI :   real := MATH_PI*2.0;
    constant  MAX_ITER:  integer := 27; -- max precision factor for cordic

    --
    -- some type declarations for cordic operations
    -- 
    
    constant KC : REAL := 6.0725293500888142e-01; -- constant for cordic
    type REAL_VECTOR is array (NATURAL range <>) of REAL; 
    type NATURAL_VECTOR is array (NATURAL range <>) of NATURAL; 
    subtype REAL_VECTOR_N is REAL_VECTOR (0 to max_iter);
    subtype REAL_ARR_2 is REAL_VECTOR (0 to 1); 
    subtype REAL_ARR_3 is REAL_VECTOR (0 to 2); 
    subtype QUADRANT is INTEGER range 0 to 3; 
    type CORDIC_MODE_TYPE is (ROTATION, VECTORING); 
  
    --
    -- auxiliary functions for cordic algorithms
    --
    function POWER_OF_2_SERIES (d : NATURAL_VECTOR; initial_value : REAL;
                number_of_values : NATURAL) return REAL_VECTOR is 

      variable v : REAL_VECTOR (0 to number_of_values);  
      variable temp : REAL := initial_value; 
      variable flag : boolean := true; 
    begin
      for i in 0 to number_of_values loop 
         v(i) := temp; 
         for p in d'range loop 
            if i = d(p) then 
               flag := false;
            end if;                  
         end loop; 
         if flag then
            temp := temp/2.0;
         end if; 
         flag := true;
      end loop; 
      return v;
    end POWER_OF_2_SERIES; 


   constant two_at_minus : REAL_VECTOR := POWER_OF_2_SERIES(
                                               NATURAL_VECTOR'(100, 90),1.0,
                                                                  MAX_ITER);  

   constant epsilon : REAL_VECTOR_N := (
                                        7.8539816339744827e-01,
                                        4.6364760900080606e-01,
                                        2.4497866312686413e-01,
                                        1.2435499454676144e-01,
                                        6.2418809995957351e-02,
                                        3.1239833430268277e-02,
                                        1.5623728620476830e-02,
                                        7.8123410601011116e-03,
                                        3.9062301319669717e-03,
                                        1.9531225164788189e-03,
                                        9.7656218955931937e-04,
                                        4.8828121119489829e-04,
                                        2.4414062014936175e-04,
                                        1.2207031189367021e-04,
                                        6.1035156174208768e-05,
                                        3.0517578115526093e-05,
                                        1.5258789061315760e-05,
                                        7.6293945311019699e-06,
                                        3.8146972656064960e-06,
                                        1.9073486328101870e-06,
                                        9.5367431640596080e-07,
                                        4.7683715820308876e-07,
                                        2.3841857910155801e-07,
                                        1.1920928955078067e-07,
                                        5.9604644775390553e-08,
                                        2.9802322387695303e-08,
                                        1.4901161193847654e-08,
                                        7.4505805969238281e-09
                                       );

   function CORDIC ( x0 : REAL;  
                     y0 : REAL;  
                     z0 : REAL;  
                      n : NATURAL;                 --       precision factor 
            CORDIC_MODE : CORDIC_MODE_TYPE         --      rotation (z -> 0) 
                                                   --  or vectoring (y -> 0) 
                    ) return REAL_ARR_3 is 
     variable x : REAL := x0;
     variable y : REAL := y0;
     variable z : REAL := z0;
     variable x_temp : REAL; 
   begin
      if CORDIC_MODE = ROTATION then 
         for k in 0 to n loop 
            x_temp := x;
            if ( z >= 0.0) then
               x := x - y * two_at_minus(k);
               y := y + x_temp * two_at_minus(k); 
               z := z - epsilon(k);
            else 
               x := x + y * two_at_minus(k);
               y := y - x_temp * two_at_minus(k);
               z := z + epsilon(k);
            end if; 
         end loop;
      else 
         for k in 0 to n loop 
            x_temp := x;
            if ( y < 0.0) then
               x := x - y * two_at_minus(k);
               y := y + x_temp * two_at_minus(k); 
               z := z - epsilon(k);
            else 
               x := x + y * two_at_minus(k);
               y := y - x_temp * two_at_minus(k);
               z := z + epsilon(k);
            end if; 
         end loop;
      end if;
      return REAL_ARR_3'(x, y, z);
   end CORDIC;          

    --
    -- non-trascendental functions
    --
    function SIGN (X: real ) return real is
    	-- returns 1.0 if X > 0.0; 0.0 if X == 0.0; -1.0 if X < 0.0
    begin
	   if  ( X > 0.0 )  then
		return 1.0;
	   elsif ( X < 0.0 )  then
		return -1.0;
	   else
		return 0.0;
	   end if;
    end SIGN; 

    function CEIL (X : real ) return real is
    	-- returns smallest integer value (as real) not less than X
	-- No conversion to an integer type is expected, so truncate cannot 
	-- overflow for large arguments.

         variable large: real  := 1073741824.0;
         type long is range -1073741824 to 1073741824;
         -- 2**30 is longer than any single-precision mantissa
         variable rd: real;
   
      begin
         if abs( X) >= large then
   	    return X;
         else
   	    rd := real ( long( X));
   	    if X > 0.0 then
   	    	if rd >= X then
   	       		return rd;
   	    	else
   	       		return rd + 1.0;
   	    	end if;
   	    elsif  X = 0.0  then
			return 0.0;
	    else
   	    	if rd <= X then
   	       		return rd;
   	    	else
   	       		return rd - 1.0;
   	    	end if;
   	    end if;
         end if;
      end CEIL;

    function FLOOR (X : real ) return real is
    	-- returns largest integer value (as real) not greater than X
   	-- No conversion to an integer type is expected, so truncate 
	-- cannot overflow for large arguments.
   	-- 
         variable large: real  := 1073741824.0;
         type long is range -1073741824 to 1073741824;
         -- 2**30 is longer than any single-precision mantissa
         variable rd: real;
   
      begin
         if abs( X ) >= large then
   	 	return X;
         else
   	 	rd := real ( long( X));
   	 	if X > 0.0 then
   	    		if rd <= X then
   	       			return rd;
   	    		else
   	       			return rd - 1.0;
   	    		end if;
   	 	elsif  X = 0.0  then
			return 0.0;
		else
   	    		if rd >= X then
   	       			return rd;
   	    		else
   	       			return rd + 1.0;
   	    		end if;
   	 	end if;
         end if;
      end FLOOR;

    function ROUND (X : real ) return real is
    	-- returns integer FLOOR(X + 0.5) if X > 0;
    	-- return integer CEIL(X - 0.5) if X < 0
    begin
	   if  X > 0.0  then
		return FLOOR(X + 0.5);
	   elsif  X < 0.0  then
		return CEIL( X - 0.5);
	   else
		return 0.0;
	   end if;
    end ROUND;
    
    function FMAX (X, Y : real ) return real is
    	-- returns the algebraically larger of X and Y
    begin
	if X > Y then
	   return X;
	else
	   return Y;
	end if;
    end FMAX;

    function FMIN (X, Y : real ) return real is
    	-- returns the algebraically smaller of X and Y
    begin
	if X < Y then
	   return X;
	else
	   return Y;
	end if;
    end FMIN;

    --
    -- Pseudo-random number generators
    --

    procedure UNIFORM(variable Seed1,Seed2:inout integer;variable X:out real) is
	-- returns a pseudo-random number with uniform distribution in the 
	-- interval (0.0, 1.0).
	-- Before the first call to UNIFORM, the seed values (Seed1, Seed2) must
	-- be initialized to values in the range [1, 2147483562] and 
	-- [1, 2147483398] respectively.  The seed values are modified after 
	-- each call to UNIFORM.
	-- This random number generator is portable for 32-bit computers, and
	-- it has period ~2.30584*(10**18) for each set of seed values.
	--
	-- For VHDL-1992, the seeds will be global variables, functions to 
	-- initialize their values (INIT_SEED) will be provided, and the UNIFORM
	-- procedure call will be modified accordingly.  
	
	variable z, k: integer;
	begin
	k := Seed1/53668;
	Seed1 := 40014 * (Seed1 - k * 53668) - k * 12211;
	
	if Seed1 < 0  then
		Seed1 := Seed1 + 2147483563;
	end if;


	k := Seed2/52774;
	Seed2 := 40692 * (Seed2 - k * 52774) - k * 3791;
	
	if Seed2 < 0  then
		Seed2 := Seed2 + 2147483399;
	end if;

	z := Seed1 - Seed2;
	if z < 1 then
		z := z + 2147483562;
	end if;

	X :=  REAL(Z)*4.656613e-10;
    end UNIFORM;


    function SRAND (seed: in integer ) return integer is
    	--
	-- sets value of seed for sequence of 
    	-- pseudo-random numbers.   
	-- Returns the value of the seed.
    	-- It uses the foreign native C function srand().
    begin
    end SRAND;

    function RAND return integer is
    	--
	-- returns an integer pseudo-random number with uniform distribution.
	-- It uses the foreign native C function rand(). 
    	-- Seed for the sequence is initialized with the
    	-- SRAND() function and value of the seed is changed every
        -- time SRAND() is called,  but it is not visible.
    	-- The range of generated values is platform dependent.
    begin
    end RAND;

    function GET_RAND_MAX  return integer is
    	--
	-- returns the upper bound of the range of the
    	-- pseudo-random numbers generated by  RAND().
    	-- The support for this function is platform dependent, and
	-- it uses foreign native C functions or constants.
	-- It may not be available in some platforms.
    	-- Note: the value of (RAND / GET_RAND_MAX) is a
    	--       pseudo-random number distributed between 0 & 1.
    begin
    end GET_RAND_MAX;

    --
    -- trascendental and trigonometric functions
    --

    function SQRT (X : real ) return real is
    	-- returns square root of X;  X >= 0
	--
	-- Computes square root using the Newton-Raphson approximation:
	-- F(n+1) = 0.5*[F(n) + x/F(n)];
	--

	constant inival: real := 1.5;
	constant eps : real := 0.000001;
	constant relative_err : real := eps*X;

	variable oldval : real ;
	variable newval : real ;

    begin
	-- check validity of argument
	if ( X < 0.0 ) then
		assert false report "X < 0 in SQRT(X)" 
			severity ERROR;
		return (0.0);
	end if;

	-- get the square root for special cases
	if X = 0.0 then
	  	return 0.0;
	else
		if ( X = 1.0 ) then
			return 1.0; -- return exact value
		end if;
	end if;

	-- get the square root for general cases
	oldval := inival;
	newval := (X/oldval + oldval)/2.0;

	while ( abs(newval -oldval) > relative_err ) loop
		oldval := newval;
		newval := (X/oldval + oldval)/2.0;
	end loop;

	return newval;
    end SQRT;

    function CBRT (X : real ) return real is
    	-- returns cube root of X
	-- Computes square root using the Newton-Raphson approximation:
	-- F(n+1) = (1/3)*[2*F(n) + x/F(n)**2];
	--

		constant inival: real := 1.5;
		constant eps : real := 0.000001;
		constant relative_err : real := eps*abs(X);

		variable xlocal : real := X;
		variable negative : boolean := X < 0.0;
		variable oldval : real ;