fpga: sim: Add PkgComplex, PkgMath, and PkgRandom

PkgComplex adds functions for doing complex arithmetic in SystemVerilog
simulation.

PkgMath provides mathematical operations and constants that aren't
built into SystemVerilog, such as a constant for pi and the function
round().

PkgRandom adds randomization functions beyond what standard Verilog
supports but that don't require any special licenses or simulators.

5 files changed, 546 insertions, 0 deletions

diff --git a/fpga/usrp3/sim/packages/Makefile.srcs b/fpga/usrp3/sim/packages/Makefile.srcs
new file mode 100644
index 000000000..4f87f3615
--- /dev/null
+++ b/fpga/usrp3/sim/packages/Makefile.srcs
@@ -0,0 +1,11 @@
+#
+# Copyright 2021 Ettus Research, A National Instruments Brand
+#
+# SPDX-License-Identifier: LGPL-3.0-or-later
+#
+
+SIM_PACKAGES_SRCS = $(abspath $(addprefix $(BASE_DIR)/../sim/packages/, \
+PkgMath.sv \
+PkgComplex.sv \
+PkgRandom.sv \
+))
diff --git a/fpga/usrp3/sim/packages/PkgComplex.sv b/fpga/usrp3/sim/packages/PkgComplex.sv
new file mode 100644
index 000000000..9fc957b97
--- /dev/null
+++ b/fpga/usrp3/sim/packages/PkgComplex.sv
@@ -0,0 +1,230 @@
+//
+// Copyright 2021 Ettus Research, A National Instruments Brand
+//
+// SPDX-License-Identifier: LGPL-3.0-or-later
+//
+// Package: PkgComplex
+//
+// Description:
+//
+//   A package for doing complex arithmetic in SystemVerilog simulations.
+//   Fixed-point operations are implemented such that results clip to the range
+//   [-1.0, 1.0) and are rounded to the nearest ULP (half ULP is rounded away
+//   from zero, following SystemVerilog rounding behavior).
+//
+
+
+package PkgComplex;
+
+  //---------------------------------------------------------------------------
+  // Type Definitions
+  //---------------------------------------------------------------------------
+
+  // Define a signed 16-bit fixed-point type with 15 fractional bits (Q0.15).
+  typedef bit signed [15:0] s16_t;
+
+  // Signed complex 16-bit data type, the standard type used by UHD and RFNoC.
+  typedef struct packed {
+    s16_t re;
+    s16_t im;
+  } sc16_t;
+
+  // Complex floating point data type.
+  typedef struct {
+    real re;
+    real im;
+  } complex_t;
+
+  // Maximum and minimum allowed by the s16 type.
+  localparam s16_t MAX_S16 = 16'h7FFF;
+  localparam s16_t MIN_S16 = 16'h8000;
+
+
+  //---------------------------------------------------------------------------
+  // Conversion Functions
+  //---------------------------------------------------------------------------
+
+  // Create an sc16 value from two s16 values.
+  function sc16_t build_sc16(s16_t x = 0, s16_t y = 0);
+    sc16_t val;
+    val.re = x;
+    val.im = y;
+    return val;
+  endfunction : build_sc16
+
+  // Create a complex value from two real values.
+  function complex_t build_complex(real x = 0.0, real y = 0.0);
+    complex_t val;
+    val.re = x;
+    val.im = y;
+    return val;
+  endfunction : build_complex
+
+  // Convert s16 to real.
+  function real s16_to_real(s16_t x);
+    return real'(x) / (2.0**15);
+  endfunction : s16_to_real
+
+  // Convert real to s16.
+  function s16_t real_to_s16(real x);
+    real val;
+    val = x * (2.0**15);
+    val = (val > MAX_S16) ? MAX_S16 : val;
+    val = (val < MIN_S16) ? MIN_S16 : val;
+    return s16_t'(val);
+  endfunction : real_to_s16
+
+  // Convert complex to sc16.
+  function sc16_t complex_to_sc16(complex_t x);
+    sc16_t val;
+    val.re = real_to_s16(x.re);
+    val.im = real_to_s16(x.im);
+    return val;
+  endfunction : complex_to_sc16
+
+  // Convert sc16 to complex.
+  function complex_t sc16_to_complex(sc16_t x);
+    complex_t val;
+    val.re = s16_to_real(x.re);
+    val.im = s16_to_real(x.im);
+    return val;
+  endfunction : sc16_to_complex
+
+  // Convert polar coordinates to a complex number. The phase should be in
+  // radians.
+  function complex_t polar_to_complex(real mag, real phase);
+    complex_t val;
+    val.re = mag * $cos(phase);
+    val.im = mag * $sin(phase);
+    return val;
+  endfunction : polar_to_complex
+
+  // Convert polar coordinates to an sc16 complex number. The phase should be
+  // in radians.
+  function sc16_t polar_to_sc16(real mag, real phase);
+    return complex_to_sc16(polar_to_complex(mag, phase));
+  endfunction : polar_to_sc16
+
+
+  //---------------------------------------------------------------------------
+  // Floating Point Complex Arithmetic
+  //---------------------------------------------------------------------------
+
+  // Add two complex numbers: x + y
+  function complex_t add(complex_t x, complex_t y);
+    complex_t val;
+    val.re = x.re + y.re;
+    val.im = x.im + y.im;
+    return val;
+  endfunction : add
+
+  // Subtract two complex numbers: x - y
+  function complex_t sub(complex_t x, complex_t y);
+    complex_t val;
+    val.re = x.re - y.re;
+    val.im = x.im - y.im;
+    return val;
+  endfunction : sub
+
+  // Multiply two complex numbers: x * y
+  function complex_t mul(complex_t x, complex_t y);
+    complex_t val;
+    val.re = x.re*y.re - x.im*y.im;
+    val.im = x.re*y.im + x.im*y.re;
+    return val;
+  endfunction : mul
+
+  // Divide two complex numbers: x / y
+  function complex_t div(complex_t x, complex_t y);
+    complex_t z;
+    z.re = (x.re*y.re + x.im*y.im) / (y.re*y.re + y.im*y.im);
+    z.im = (x.im*y.re + x.re*y.im) / (y.re*y.re + y.im*y.im);
+    return z;
+  endfunction : div
+
+  // Compute the exponential: e^x
+  function complex_t exp(complex_t x);
+    complex_t val;
+    // exp(a+jb) = exp(a)*exp(jb) = exp(a)*(cos(b) + j*sin(b))
+    val.re = $exp(x.re)*$cos(x.im);
+    val.im = $exp(x.re)*$sin(x.im);
+    return val;
+  endfunction : exp
+
+  // Compute the sine: sin(x)
+  function complex_t sin(complex_t x);
+    complex_t val;
+    val.re = $sin(x.re)*$cosh(x.im);
+    val.im = $cos(x.re)*$sinh(x.im);
+    return val;
+  endfunction : sin
+
+  // Compute the cosine: cos(x)
+  function complex_t cos(complex_t x);
+    complex_t val;
+    val.re =      $cos(x.re)*$cosh(x.im);
+    val.im = -1.0*$sin(x.re)*$sinh(x.im);
+    return val;
+  endfunction : cos
+
+  // Compute the magnitude/modulus/absolute value: |x|
+  function real mag(complex_t x);
+    return $sqrt(x.re*x.re + x.im*x.im);
+  endfunction : mag
+
+  // Compute the phase/argument: arg(x)
+  function real arg(complex_t x);
+    return $atan2(x.im, x.re);
+  endfunction : arg
+
+
+  //---------------------------------------------------------------------------
+  // Fixed-Point Complex Arithmetic
+  //---------------------------------------------------------------------------
+
+  // Add two complex numbers: x + y
+  function sc16_t add_sc16(sc16_t x, sc16_t y);
+    return complex_to_sc16(add(sc16_to_complex(x), sc16_to_complex(y)));
+  endfunction : add_sc16
+
+  // Subtract two complex numbers: x - y
+  function sc16_t sub_sc16(sc16_t x, sc16_t y);
+    return complex_to_sc16(sub(sc16_to_complex(x), sc16_to_complex(y)));
+  endfunction : sub_sc16
+
+  // Multiply two complex numbers: x * y
+  function sc16_t mul_sc16(sc16_t x, sc16_t y);
+    return complex_to_sc16(mul(sc16_to_complex(x), sc16_to_complex(y)));
+  endfunction : mul_sc16
+
+  // Divide two complex numbers: x / y
+  function sc16_t div_sc16(sc16_t x, sc16_t y);
+    return complex_to_sc16(div(sc16_to_complex(x), sc16_to_complex(y)));
+  endfunction : div_sc16
+
+  // Compute the exponential: e^x
+  function sc16_t exp_sc16(sc16_t x);
+    return complex_to_sc16(exp(sc16_to_complex(x)));
+  endfunction : exp_sc16
+
+  // Compute the sine: sin(x)
+  function sc16_t sin_sc16(sc16_t x);
+    return complex_to_sc16(sin(sc16_to_complex(x)));
+  endfunction : sin_sc16
+
+  // compute the cosine: cos(x)
+  function sc16_t cos_sc16(sc16_t x);
+    return complex_to_sc16(cos(sc16_to_complex(x)));
+  endfunction : cos_sc16
+
+  // Compute the magnitude/modulus/absolute value: |x|
+  function s16_t mag_sc16(sc16_t x);
+    return real_to_s16(mag(sc16_to_complex(x)));
+  endfunction : mag_sc16
+
+  // Compute the phase/argument: arg(x)
+  function s16_t arg_sc16(sc16_t x);
+    return real_to_s16(arg(sc16_to_complex(x)));
+  endfunction : arg_sc16
+
+endpackage : PkgComplex
diff --git a/fpga/usrp3/sim/packages/PkgMath.sv b/fpga/usrp3/sim/packages/PkgMath.sv
new file mode 100644
index 000000000..db5c15840
--- /dev/null
+++ b/fpga/usrp3/sim/packages/PkgMath.sv
@@ -0,0 +1,157 @@
+//
+// Copyright 2021 Ettus Research, A National Instruments Brand
+//
+// SPDX-License-Identifier: LGPL-3.0-or-later
+//
+// Package: PkgMath
+//
+// Description:
+//
+//   SystemVerilog supports many Math functions. This adds a few that it
+//   doesn't have built in, as well as useful mathematical constants.
+//
+//   SystemVerilog has built-in support for the following:
+//
+//     $clog2   $asin
+//     $ln      $acos
+//     $log10   $atan
+//     $exp     $atan2
+//     $sqrt    $hypot
+//     $pow     $sinh
+//     $floor   $cosh
+//     $ceil    $tanh
+//     $sin     $asinh
+//     $cos     $acosh
+//     $tan     $atanh
+//
+
+package PkgMath;
+
+  //---------------------------------------------------------------------------
+  // Constants
+  //---------------------------------------------------------------------------
+
+  localparam real PI  = 2*$acos(0.0);
+  localparam real TAU = 4*$acos(0.0);
+
+  localparam real PHI = (1 + $sqrt(5.0)) / 2.0;
+  localparam real E   = $exp(1);
+
+  localparam byte     BYTE_MAX  =  8'sh7F;
+  localparam byte     BYTE_MIN  =  8'sh80;
+  localparam shortint SHORT_MAX = 16'sh7FFF;
+  localparam shortint SHORT_MIN = 16'sh8000;
+  localparam int      INT_MAX   = 32'sh7FFFFFFF;
+  localparam int      INT_MIN   = 32'sh80000000;
+  localparam longint  LONG_MAX  = 64'sh7FFFFFFFFFFFFFFF;
+  localparam longint  LONG_MIN  = 64'sh8000000000000000;
+
+  localparam byte     unsigned UBYTE_MAX  =  8'hFF;
+  localparam byte     unsigned UBYTE_MIN  =  8'h00;
+  localparam shortint unsigned USHORT_MAX = 16'hFFFF;
+  localparam shortint unsigned USHORT_MIN = 16'h0000;
+  localparam int      unsigned UINT_MAX   = 32'hFFFFFFFF;
+  localparam int      unsigned UINT_MIN   = 32'h00000000;
+  localparam longint  unsigned ULONG_MAX  = 64'hFFFFFFFFFFFFFFFF;
+  localparam longint  unsigned ULONG_MIN  = 64'h0000000000000000;
+
+
+  //---------------------------------------------------------------------------
+  // Functions (For real data types)
+  //---------------------------------------------------------------------------
+
+  // Return the absolute value
+  function automatic real abs(real num);
+    if (num < 0) return -1.0*num;
+    return num;
+  endfunction : abs
+
+  // Round a float to the nearest whole number, rounding away from zero for 0.5
+  // (same as C++ and default SystemVerilog behavior).
+  function automatic real round(real num);
+    if (num >= 0) begin
+      // Round toward +inf
+      if (num - $floor(num) < 0.5) return $floor(num);
+      return $ceil(num);
+    end else begin
+      // Round toward -inf
+      if (num - $floor(num) <= 0.5) return $floor(num);
+      return $ceil(num);
+    end
+
+  endfunction : round
+
+  // Round a float to the nearest value having bits to the right of the binary
+  // point. For example:
+  //
+  //   1.2265625 (0b1.0011101) --> 3 bits --> 1.25000 (0b1.0100000)
+  //   1.2265625 (0b1.0011101) --> 5 bits --> 1.21875 (0b1.0011100)
+  //
+  function automatic real round_bits(real num, int unsigned bits);
+    return round(num * 2.0**bits) / (2.0**bits);
+  endfunction : round_bits
+
+  // Return the sign of num as +1.0 or -1.0;
+  function automatic real sign(real num);
+    if (num < 0.0) return -1.0;
+    return 1.0;
+  endfunction : sign;
+
+  // Return the modulus (remainder) of a / b, with the sign of the numerator.
+  // This should match the C++ standard library std::fmod() behavior, as well
+  // as SystemVerilog % operator with integer values.
+  function automatic real fmod(real a, real b);
+    a = abs(a);
+    b = abs(b);
+    return sign(b) * (a - ($floor(a / b) * b));
+  endfunction : fmod
+
+  // Return the (remainder) of a / b, where the quotient is rounded to the
+  // nearest integer. This should approximate the C++ standard library
+  // std::remainder() behavior.
+  function automatic real remainder(real a, real b);
+    return a - round(a/b)*b;
+  endfunction : remainder
+
+  // Return the maximum of a and b.
+  function automatic real fmax(real a, real b);
+    if (a > b) return a;
+    return b;
+  endfunction : fmax
+
+  // Return the minimum of a and b.
+  function automatic real fmin(real a, real b);
+    if (a < b) return a;
+    return b;
+  endfunction : fmin
+
+
+  //---------------------------------------------------------------------------
+  // Template Functions (For any data type)
+  //---------------------------------------------------------------------------
+
+  class Math #(type T);
+
+    static function T abs(T num);
+      if (num < 0) return -num;
+      return num;
+    endfunction : abs
+
+    static function T sign(T num);
+      if (num < 0) return -1;
+      return 1;
+    endfunction : sign
+
+    static function T max(T a, T b);
+      if (a > b) return a;
+      else return b;
+    endfunction : max
+
+    static function T min(T a, T b);
+      if (a < b) return a;
+      else return b;
+    endfunction : min
+
+  endclass : Math
+
+endpackage : PkgMath
diff --git a/fpga/usrp3/sim/packages/PkgRandom.sv b/fpga/usrp3/sim/packages/PkgRandom.sv
new file mode 100644
index 000000000..94e08d447
--- /dev/null
+++ b/fpga/usrp3/sim/packages/PkgRandom.sv
@@ -0,0 +1,146 @@
+//
+// Copyright 2021 Ettus Research, A National Instruments Brand
+//
+// SPDX-License-Identifier: LGPL-3.0-or-later
+//
+// Package: PkgRandom
+//
+// Description:
+//
+//   SystemVerilog has great randomization support, but some features require a
+//   more expensive license or aren't supported by all tools. This package
+//   tries to fill that gap by providing some useful randomization functions
+//   beyond what's supported by standard Verilog.
+//
+
+package PkgRandom;
+
+  import PkgMath::*;
+
+  //---------------------------------------------------------------------------
+  // Functions
+  //---------------------------------------------------------------------------
+
+  // Return a real value in the range [0,max), where max is 1.0 by default.
+  function automatic real frand(real max = 1.0);
+    bit [63:0] real_bits;
+    real num;
+
+    // Build a double-precision floating point value per IEEE-754 standard,
+    // which SystemVerilog follows.
+
+    // Positive, with exponent 0
+    real_bits[63:52] = 12'h3FF;
+
+    // Mantissa in the range [1.0, 2.0). The leading 1 in the mantissa is
+    // implied by the floating point format.
+    real_bits[31: 0] = $urandom();
+    real_bits[51:32] = $urandom();
+
+    // Compensate for the implied leading 1 in the mantissa by subtracting 1.
+    num = $bitstoreal(real_bits) - 1.0;
+
+    // Scale the result to return a value in the desired range.
+    return num * max;
+  endfunction : frand
+
+
+  // Return a real value in the range [a,b), [b,a), or [0,a) depending on
+  // whether a or b is larger and whether b is provided. This matches the
+  // behavior of $urandom_range().
+  //
+  //   frand_range(1.0, 2.0) -> Random value in the range [1,2)
+  //   frand_range(2.0, 1.0) -> Random value in the range [1,2)
+  //   frand_range(1.0)      -> Random value in the range [0,1)
+  //
+  function automatic real frand_range(real a = 1.0, real b = 0.0);
+    if (a > b) return b + frand(a - b);
+    if (b > a) return a + frand(b - a);
+    return a;
+  endfunction : frand_range
+
+
+  // Return a real value with a normal distribution, having the mean value mu
+  // and standard deviation sigma.
+  function automatic real frandn(real sigma = 1.0, real mu = 0.0);
+    // Use the Box-Muller transform to convert uniform random variables to a
+    // Gaussian one.
+    return sigma*$sqrt(-2.0*$ln(frand())) * $cos(TAU*frand()) + mu;
+  endfunction : frandn
+
+
+  //---------------------------------------------------------------------------
+  // Template Functions
+  //---------------------------------------------------------------------------
+
+  class Rand #(WIDTH = 64);
+
+    // These are static class functions. They can be called directly, as in:
+    //
+    //   Rand#(N)::rand_bit()
+    //
+    // Or, you can declare an object reference, as in:
+    //
+    //   Rand #(N) rand;
+    //   rand.rand_bit();
+
+    typedef bit        [WIDTH-1:0] unsigned_t;
+    typedef bit signed [WIDTH-1:0] signed_t;
+
+
+    // Returns a WIDTH-bit random bit packed array.
+    static function unsigned_t rand_bit();
+      unsigned_t result;
+      int num_rand32 = (WIDTH + 31) / 32;
+      for (int i = 0; i < num_rand32; i++) begin
+        result = {result, $urandom()};
+      end
+      return result;
+    endfunction : rand_bit
+
+
+    // Returns a WIDTH-bit random number in the UNSIGNED range [a,b], [b,a], or
+    // [0,a] depending on whether a or b is greater and if b is provided. This
+    // is equivalent to $urandom_range() but works with any length.
+    static function bit [WIDTH-1:0] rand_bit_range(
+      unsigned_t a = {WIDTH{1'b1}},
+      unsigned_t b = 0
+    );
+      unsigned_t num;
+      int num_bits;
+      if (a > b) begin
+        // Swap a and b
+        unsigned_t temp;
+        temp = a;
+        a = b;
+        b = temp;
+      end
+      num_bits = $clog2(b - a + unsigned_t'{1});
+      do begin
+        num = a + (rand_bit() & ((unsigned_t'{1} << num_bits) - 1));
+      end while (num > b);
+      return num;
+    endfunction : rand_bit_range
+
+
+    // Returns a random number in the given SIGNED range. Behavior is the same
+    // as rand_bit_range(), bunsigned_t treats the range values as SIGNED numbers.
+    static function signed_t rand_sbit_range(
+      signed_t a = {1'b0, {WIDTH{1'b1}}},
+      signed_t b = 0
+    );
+      if (a > b) return b + $signed(rand_bit_range(0, a-b));
+      if (b > a) return a + $signed(rand_bit_range(0, b-a));
+      return a;
+    endfunction : rand_sbit_range
+
+
+    // Rand#(WIDTH)::rand_logic() returns a WIDTH-bit random logic packed
+    // array. Each bit will be 0 or 1 with equal probability (not X or Z).
+    static function logic [WIDTH-1:0] rand_logic();
+      return rand_bit();
+    endfunction : rand_logic
+
+  endclass : Rand
+
+endpackage : PkgRandom
diff --git a/fpga/usrp3/tools/make/viv_sim_preamble.mak b/fpga/usrp3/tools/make/viv_sim_preamble.mak
index afeca0b3c..514b06a0a 100644
--- a/fpga/usrp3/tools/make/viv_sim_preamble.mak
+++ b/fpga/usrp3/tools/make/viv_sim_preamble.mak
@@ -30,12 +30,14 @@ include $(BASE_DIR)/../sim/general/Makefile.srcs
 include $(BASE_DIR)/../sim/axi/Makefile.srcs
 include $(BASE_DIR)/../sim/control/Makefile.srcs
 include $(BASE_DIR)/../sim/rfnoc/Makefile.srcs
+include $(BASE_DIR)/../sim/packages/Makefile.srcs
 
 INC_SRCS = $(abspath \
 $(SIM_GENERAL_SRCS) \
 $(SIM_AXI_SRCS) \
 $(SIM_CONTROL_SRCS) \
 $(SIM_RFNOC_SRCS) \
+$(SIM_PACKAGES_SRCS) \
 )
 
 # Predeclare RFNOC_OOT_SRCS to make sure it's not recursively expanded