git.gag.com Git - debian/gnuradio/blob - usrp2/fpga/models/math_real.v

   1 /*\r
   2  * This is a general recreation of the VHDL ieee.math_real package.\r
   3  */\r
   4  \r
   5 module math_real ;\r
   6    // Constants for use below and for general reference\r
   7    // TODO: Bring it out to 12 (or more???) places beyond the decimal?\r
   8    localparam MATH_E            =  2.7182818284;\r
   9    localparam MATH_1_OVER_E     =  0.3678794411;\r
  10    localparam MATH_PI           =  3.1415926536;\r
  11    localparam MATH_2_PI         =  6.2831853071;\r
  12    localparam MATH_1_OVER_PI    =  0.3183098861;\r
  13    localparam MATH_PI_OVER_2    =  1.5707963267;\r
  14    localparam MATH_PI_OVER_3    =  1.0471975511;\r
  15    localparam MATH_PI_OVER_4    =  0.7853981633;\r
  16    localparam MATH_3_PI_OVER_2  =  4.7123889803;\r
  17    localparam MATH_LOG_OF_2     =  0.6931471805;\r
  18    localparam MATH_LOG_OF_10    =  2.3025850929;\r
  19    localparam MATH_LOG2_OF_E    =  1.4426950408;\r
  20    localparam MATH_LOG10_OF_E   =  0.4342944819;\r
  21    localparam MATH_SQRT_2       =  1.4142135623;\r
  22    localparam MATH_1_OVER_SQRT_2=  0.7071067811;\r
  23    localparam MATH_SQRT_PI      =  1.7724538509;\r
  24    localparam MATH_DEG_TO_RAD   =  0.0174532925;\r
  25    localparam MATH_RAD_TO_DEG   =  57.2957795130;\r
  26 \r
  27    // The number of iterations to do for the Taylor series approximations\r
  28    localparam EXPLOG_ITERATIONS =  19;\r
  29    localparam COS_ITERATIONS    =  8;\r
  30    \r
  31    /* Conversion Routines */\r
  32     \r
  33    // Return the sign of a particular number.\r
  34    function real sign ;\r
  35       input real x ;\r
  36       begin\r
  37          sign = x < 0.0 ? 1.0 : 0.0 ;\r
  38       end\r
  39    endfunction\r
  40    \r
  41    // Return the trunc function of a number\r
  42    function real trunc ;\r
  43       input real x ;\r
  44       begin\r
  45          trunc = x - mod(x,1.0) ;\r
  46       end\r
  47    endfunction\r
  48    \r
  49    // Return the ceiling function of a number.\r
  50    function real ceil ;\r
  51       input real x ;\r
  52       real retval ;\r
  53       begin\r
  54          retval = mod(x,1.0) ;\r
  55          if( retval != 0.0 && x > 0.0 )  retval = x+1.0 ;\r
  56          else retval = x ;\r
  57          ceil = trunc(retval) ;\r
  58       end\r
  59    endfunction\r
  60    \r
  61    // Return the floor function of a number\r
  62    function real floor ;\r
  63       input real x ;\r
  64       real retval ;\r
  65       begin\r
  66          retval = mod(x,1.0) ;\r
  67          if( retval != 0.0 && x < 0.0 ) retval = x - 1.0 ;\r
  68          else retval = x ;\r
  69          floor = trunc(retval) ;\r
  70       end\r
  71    endfunction\r
  72    \r
  73    // Return the round function of a number\r
  74    function real round ;\r
  75       input real x ;\r
  76       real retval ;\r
  77       begin\r
  78          retval = x > 0.0 ? x + 0.5 : x - 0.5 ;\r
  79          round = trunc(retval) ;\r
  80       end\r
  81    endfunction\r
  82    \r
  83    // Return the fractional remainder of (x mod m)\r
  84    function real mod ;\r
  85       input real x ;\r
  86       input real m ;\r
  87       real retval ;\r
  88       begin\r
  89          retval = x ;\r
  90          if( retval > m ) begin\r
  91             while( retval > m ) begin\r
  92                retval = retval - m ;\r
  93             end\r
  94          end\r
  95          else begin\r
  96             while( retval < -m ) begin\r
  97                retval = retval + m ;\r
  98             end\r
  99          end\r
 100          mod = retval ;\r
 101       end\r
 102    endfunction\r
 103    \r
 104    // Return the max between two real numbers\r
 105    function real realmax ;\r
 106       input real x ;\r
 107       input real y ;\r
 108       begin\r
 109          realmax = x > y ? x : y ;\r
 110       end\r
 111    endfunction\r
 112    \r
 113    // Return the min between two real numbers\r
 114    function real realmin ;\r
 115       input real x ;\r
 116       input real y ;\r
 117       begin\r
 118          realmin = x > y ? y : x ;\r
 119       end\r
 120    endfunction\r
 121    \r
 122    /* Random Numbers */\r
 123    \r
 124    // Generate Gaussian distributed variables\r
 125    function real gaussian ;\r
 126       input real mean ;\r
 127       input real var ;\r
 128       real u1, u2, v1, v2, s ;\r
 129       begin\r
 130          s = 1.0 ;\r
 131          while( s >= 1.0  ) begin\r
 132             // Two random numbers between 0 and 1\r
 133             u1 = $random/4294967296.0 + 0.5 ;\r
 134             u2 = $random/4294967296.0 + 0.5 ;\r
 135             // Adjust to be between -1,1\r
 136             v1 = 2*u1-1.0 ;\r
 137             v2 = 2*u2-1.0 ;\r
 138             // Polar mag squared\r
 139             s = (v1*v1 + v2*v2) ;\r
 140          end\r
 141          gaussian = mean + sqrt((-2.0*log(s))/s) * v1 * sqrt(var) ;\r
 142          // gaussian2 = mean + sqrt(-2*log(s)/s)*v2 * sqrt(var) ;\r
 143       end\r
 144    endfunction\r
 145    \r
 146    /* Roots and Log Functions */\r
 147    \r
 148    // Return the square root of a number\r
 149    function real sqrt ;\r
 150       input real x ;\r
 151       real retval ;\r
 152       begin\r
 153          sqrt = (x == 0.0) ? 0.0 : powr(x,0.5) ;\r
 154       end\r
 155    endfunction\r
 156    \r
 157    // Return the cube root of a number\r
 158    function real cbrt ;\r
 159       input real x ;\r
 160       real retval ;\r
 161       begin\r
 162          cbrt = (x == 0.0) ? 0.0 : powr(x,1.0/3.0) ;\r
 163       end\r
 164    endfunction\r
 165    \r
 166    // Return the absolute value of a real value\r
 167    function real abs ;\r
 168       input real x ;\r
 169       begin\r
 170          abs = (x > 0.0) ? x : -x ;\r
 171       end\r
 172    endfunction \r
 173    \r
 174    // Return a real value raised to an integer power\r
 175    function real pow ;\r
 176       input real b ;\r
 177       input integer x ;\r
 178       integer absx ;\r
 179       real retval ;\r
 180       begin\r
 181          retval = 1.0 ;\r
 182          absx = abs(x) ;\r
 183          repeat(absx) begin\r
 184             retval = b*retval ;\r
 185          end\r
 186          pow = x < 0 ? (1.0/retval) : retval ;\r
 187       end\r
 188    endfunction\r
 189    \r
 190    // Return a real value raised to a real power\r
 191    function real powr ;\r
 192       input real b ;\r
 193       input real x ;\r
 194       begin\r
 195          powr = exp(x*log(b)) ;\r
 196       end\r
 197    endfunction\r
 198    \r
 199    // Return the evaluation of e^x where e is the natural logarithm base\r
 200    // NOTE: This is the Taylor series expansion of e^x\r
 201    function real exp ;\r
 202       input real x ;\r
 203       real retval ;\r
 204       integer i ;\r
 205       real nm1_fact ;\r
 206       real powm1 ;\r
 207       begin\r
 208          nm1_fact = 1.0 ;\r
 209          powm1 = 1.0 ;\r
 210          retval = 1.0 ;\r
 211          for( i = 1 ; i < EXPLOG_ITERATIONS ; i = i + 1 ) begin\r
 212             powm1 = x*powm1 ;\r
 213             nm1_fact = nm1_fact * i ;\r
 214             retval = retval + powm1/nm1_fact ;\r
 215          end\r
 216          exp = retval ;\r
 217       end\r
 218    endfunction\r
 219    \r
 220    // Return the evaluation log(x)\r
 221    function real log ;\r
 222       input real x ;\r
 223       integer i ;\r
 224       real whole ;\r
 225       real xm1oxp1 ;\r
 226       real retval ;\r
 227       real newx ;\r
 228       begin\r
 229          retval = 0.0 ;\r
 230          whole = 0.0 ;\r
 231          newx = x ;\r
 232          while( newx > MATH_E ) begin\r
 233             whole = whole + 1.0 ;\r
 234             newx = newx / MATH_E ;\r
 235          end\r
 236          xm1oxp1 = (newx-1.0)/(newx+1.0) ;\r
 237          for( i = 0 ; i < EXPLOG_ITERATIONS ; i = i + 1 ) begin\r
 238             retval = retval + pow(xm1oxp1,2*i+1)/(2.0*i+1.0) ;\r
 239          end\r
 240          log = whole+2.0*retval ;\r
 241       end\r
 242    endfunction\r
 243    \r
 244    // Return the evaluation ln(x) (same as log(x))\r
 245    function real ln ;\r
 246       input real x ;\r
 247       begin\r
 248          ln = log(x) ;\r
 249       end\r
 250    endfunction\r
 251    \r
 252    // Return the evaluation log_2(x)\r
 253    function real log2 ;\r
 254       input real x ;\r
 255       begin\r
 256          log2 = log(x)/MATH_LOG_OF_2 ;\r
 257       end\r
 258    endfunction\r
 259    \r
 260    function real log10 ;\r
 261       input real x ;\r
 262       begin\r
 263          log10 = log(x)/MATH_LOG_OF_10 ;\r
 264       end\r
 265    endfunction\r
 266    \r
 267    function real log_base ;\r
 268       input real x ;\r
 269       input real b ;\r
 270       begin\r
 271          log_base = log(x)/log(b) ;\r
 272       end\r
 273    endfunction\r
 274    \r
 275    /* Trigonometric Functions */\r
 276    \r
 277    // Internal function to reduce a value to be between [-pi:pi]\r
 278    function real reduce ;\r
 279       input real x ;\r
 280       real retval ;\r
 281       begin\r
 282          retval = x ;\r
 283          while( abs(retval) > MATH_PI ) begin\r
 284             retval = retval > MATH_PI ? \r
 285                      (retval - MATH_2_PI) : \r
 286                      (retval + MATH_2_PI) ;\r
 287          end\r
 288          reduce = retval ;\r
 289       end\r
 290    endfunction\r
 291    \r
 292    // Return the cos of a number in radians\r
 293    function real cos ;\r
 294       input real x ;\r
 295       integer i ;\r
 296       integer sign ;\r
 297       real newx ;\r
 298       real retval ;\r
 299       real xsqnm1 ;\r
 300       real twonm1fact ;\r
 301       begin\r
 302          newx = reduce(x) ;\r
 303          xsqnm1 = 1.0 ;\r
 304          twonm1fact = 1.0 ;\r
 305          retval = 1.0 ;\r
 306          for( i = 1 ; i < COS_ITERATIONS ; i = i + 1 ) begin\r
 307             sign = -2*(i % 2)+1 ;\r
 308             xsqnm1 = xsqnm1*newx*newx ;\r
 309             twonm1fact = twonm1fact * (2.0*i) * (2.0*i-1.0) ;\r
 310             retval = retval + sign*(xsqnm1/twonm1fact) ; \r
 311          end\r
 312          cos = retval ;\r
 313       end\r
 314    endfunction\r
 315    \r
 316    // Return the sin of a number in radians\r
 317    function real sin ;\r
 318       input real x ;\r
 319       begin\r
 320          sin = cos(x - MATH_PI_OVER_2) ;\r
 321       end\r
 322    endfunction\r
 323    \r
 324    // Return the tan of a number in radians\r
 325    function real tan ;\r
 326       input real x ;\r
 327       begin\r
 328          tan = sin(x) / cos(x) ;\r
 329       end\r
 330    endfunction\r
 331    \r
 332    // Return the arcsin in radians of a number\r
 333    function real arcsin ;\r
 334       input real x ;\r
 335       begin\r
 336          arcsin = 2.0*arctan(x/(1.0+sqrt(1.0-x*x))) ;\r
 337       end\r
 338    endfunction\r
 339    \r
 340    // Return the arccos in radians of a number\r
 341    function real arccos ;\r
 342       input real x ;\r
 343       begin\r
 344          arccos = MATH_PI_OVER_2-arcsin(x) ;\r
 345       end\r
 346    endfunction\r
 347    \r
 348    // Return the arctan in radians of a number\r
 349    // TODO: Make sure this REALLY does work as it is supposed to!\r
 350    function real arctan ;\r
 351       input real x ;\r
 352       real retval ;\r
 353       real y ;\r
 354       real newx ;\r
 355       real twoiotwoip1 ;\r
 356       integer i ;\r
 357       integer mult ;\r
 358       begin\r
 359          retval = 1.0 ;\r
 360          twoiotwoip1 = 1.0 ;\r
 361          mult = 1 ;\r
 362          newx = abs(x) ;\r
 363          while( newx > 1.0 ) begin\r
 364             mult = mult*2 ;\r
 365             newx = newx/(1.0+sqrt(1.0+newx*newx)) ;\r
 366          end\r
 367          y = 1.0 ;\r
 368          for( i = 1 ; i < 2*COS_ITERATIONS ; i = i + 1 ) begin\r
 369             y = y*((newx*newx)/(1+newx*newx)) ;\r
 370             twoiotwoip1 = twoiotwoip1 * (2.0*i)/(2.0*i+1.0) ;\r
 371             retval = retval + twoiotwoip1*y ;\r
 372          end\r
 373          retval = retval * (newx/(1+newx*newx)) ;\r
 374          retval = retval * mult ;\r
 375          \r
 376          arctan = (x > 0.0) ? retval : -retval ;\r
 377       end\r
 378    endfunction\r
 379    \r
 380    // Return the arctan in radians of a ratio x/y\r
 381    // TODO: Test to make sure this works as it is supposed to!\r
 382    function real arctan_xy ;\r
 383       input real x ;\r
 384       input real y ;\r
 385       real retval ;\r
 386       begin\r
 387          retval = 0.0 ;\r
 388          if( x < 0.0 ) retval = MATH_PI - arctan(-abs(y)/x) ;\r
 389          else if( x > 0.0 ) retval = arctan(abs(y)/x) ;\r
 390          else if( x == 0.0 ) retval = MATH_PI_OVER_2 ;\r
 391          arctan_xy = (y < 0.0) ? -retval : retval ;\r
 392       end\r
 393    endfunction\r
 394    \r
 395    /* Hyperbolic Functions */\r
 396    \r
 397    // Return the sinh of a number\r
 398    function real sinh ;\r
 399       input real x ;\r
 400       begin\r
 401          sinh = (exp(x) - exp(-x))/2.0 ;\r
 402       end\r
 403    endfunction\r
 404    \r
 405    // Return the cosh of a number\r
 406    function real cosh ;\r
 407       input real x ;\r
 408       begin\r
 409          cosh = (exp(x) + exp(-x))/2.0 ;\r
 410       end\r
 411    endfunction\r
 412    \r
 413    // Return the tanh of a number\r
 414    function real tanh ;\r
 415       input real x ;\r
 416       real e2x ;\r
 417       begin\r
 418          e2x = exp(2.0*x) ;\r
 419          tanh = (e2x+1.0)/(e2x-1.0) ;\r
 420       end\r
 421    endfunction\r
 422    \r
 423    // Return the arcsinh of a number\r
 424    function real arcsinh ;\r
 425       input real x ;\r
 426       begin\r
 427          arcsinh = log(x+sqrt(x*x+1.0)) ;\r
 428       end\r
 429    endfunction\r
 430    \r
 431    // Return the arccosh of a number\r
 432    function real arccosh ;\r
 433       input real x ;\r
 434       begin\r
 435          arccosh = ln(x+sqrt(x*x-1.0)) ;\r
 436       end\r
 437    endfunction\r
 438    \r
 439    // Return the arctanh of a number\r
 440    function real arctanh ;\r
 441       input real x ;\r
 442       begin\r
 443          arctanh = 0.5*ln((1.0+x)/(1.0-x)) ;\r
 444       end\r
 445    endfunction\r
 446    /*\r
 447     initial begin\r
 448     $display( "cos(MATH_PI_OVER_3): %f", cos(MATH_PI_OVER_3) ) ;\r
 449     $display( "sin(MATH_PI_OVER_3): %f", sin(MATH_PI_OVER_3) ) ;\r
 450     $display( "sign(-10): %f", sign(-10) ) ;\r
 451     $display( "realmax(MATH_PI,MATH_E): %f", realmax(MATH_PI,MATH_E) ) ;\r
 452     $display( "realmin(MATH_PI,MATH_E): %f", realmin(MATH_PI,MATH_E) ) ;\r
 453     $display( "mod(MATH_PI,MATH_E): %f", mod(MATH_PI,MATH_E) ) ;\r
 454     $display( "ceil(-MATH_PI): %f", ceil(-MATH_PI) ) ;\r
 455     $display( "ceil(4.0): %f", ceil(4.0) ) ;\r
 456     $display( "ceil(3.99999999999999): %f", ceil(3.99999999999999) ) ;\r
 457     $display( "pow(MATH_PI,2): %f", pow(MATH_PI,2) ) ;\r
 458     $display( "gaussian(1.0,1.0): %f", gaussian(1.0,1.0) ) ;\r
 459     $display( "round(MATH_PI): %f", round(MATH_PI) ) ;\r
 460     $display( "trunc(-MATH_PI): %f", trunc(-MATH_PI) ) ;\r
 461     $display( "ceil(-MATH_PI): %f", ceil(-MATH_PI) ) ;\r
 462     $display( "floor(MATH_PI): %f", floor(MATH_PI) ) ;\r
 463     $display( "round(e): %f", round(MATH_E)) ;\r
 464     $display( "ceil(-e): %f", ceil(-MATH_E)) ;\r
 465     $display( "exp(MATH_PI): %f", exp(MATH_PI) ) ;\r
 466     $display( "log2(MATH_PI): %f", log2(MATH_PI) ) ;\r
 467     $display( "log_base(pow(2,32),2): %f", log_base(pow(2,32),2) ) ;\r
 468     $display( "ln(0.1): %f", log(0.1) ) ;\r
 469     $display( "cbrt(7): %f", cbrt(7) ) ;\r
 470     $display( "cos(MATH_2_PI): %f", cos(20*MATH_2_PI) ) ;\r
 471     $display( "sin(-MATH_2_PI): %f", sin(-50*MATH_2_PI) ) ;\r
 472     $display( "sinh(MATH_E): %f", sinh(MATH_E) ) ;\r
 473     $display( "cosh(MATH_2_PI): %f", cosh(MATH_2_PI) ) ;\r
 474     $display( "arctan_xy(-4,3): %f", arctan_xy(-4,3) ) ;\r
 475     $display( "arctan(MATH_PI): %f", arctan(MATH_PI) ) ;\r
 476     $display( "arctan(-MATH_E/2): %f", arctan(-MATH_E/2) ) ;\r
 477     $display( "arctan(MATH_PI_OVER_2): %f", arctan(MATH_PI_OVER_2) ) ;\r
 478     $display( "arctan(1/7) = %f", arctan(1.0/7.0) ) ;\r
 479     $display( "arctan(3/79) = %f", arctan(3.0/79.0) ) ;\r
 480     $display( "pi/4 ?= %f", 5*arctan(1.0/7.0)+2*arctan(3.0/79.0) ) ;\r
 481     $display( "arcsin(1.0): %f", arcsin(1.0) ) ;\r
 482     $display( "cos(pi/2): %f", cos(MATH_PI_OVER_2)) ;\r
 483     $display( "arccos(cos(pi/2)): %f", arccos(cos(MATH_PI_OVER_2)) ) ;\r
 484     $display( "cos(0): %f", cos(0) ) ;\r
 485     $display( "cos(MATH_PI_OVER_4): %f", cos(MATH_PI_OVER_4) ) ;\r
 486     $display( "cos(MATH_PI_OVER_2): %f", cos(MATH_PI_OVER_2) ) ;\r
 487     $display( "cos(3*MATH_PI_OVER_4): %f", cos(3*MATH_PI_OVER_4) ) ;\r
 488     $display( "cos(MATH_PI): %f", cos(MATH_PI) ) ;\r
 489     $display( "cos(5*MATH_PI_OVER_4): %f", cos(5*MATH_PI_OVER_4) ) ;\r
 490     $display( "cos(6*MATH_PI_OVER_4): %f", cos(6*MATH_PI_OVER_4) ) ;\r
 491     $display( "cos(7*MATH_PI_OVER_4): %f", cos(7*MATH_PI_OVER_4) ) ;\r
 492     $display( "cos(8*MATH_PI_OVER_4): %f", cos(8*MATH_PI_OVER_4) ) ;\r
 493     end*/\r
 494    \r
 495 endmodule\r