/************************************************************************
// Title of Application: Butterworth IIR LP filter using VpComplex 
//                       operands and finding optimal width/formats for 
//                       the datacontainers used.
//
// This example uses FinSimMath's VpComplex and VpPolar scalar and
// vector types in conjunction with $VpSin, $VpFft, $VpIfft system
// tasks to demonstrate the implementation of a low pas Butterworth
// filter and other DSP processing.
//
// VpComplex and VpPolar are types predefined in FinSimMath, to contain
// pairs of variable precision registers, having the information about
// their fomats associated via a descriptor at runtime.
// Assignments and arithmetic operators handle objects of these types
// as complex numbers in carthesian and polar coordinates representation,
// respectively. Their subfields can be accessed by name as Re and Im or
// Mag and Ang with the obvious meaning.
//
// FinSimMath also supports VpFComplex and VpFPolar types, which contain real
// numbers.
//
//
//Notes:
// i)   FinSimMath is an extension of the Verilog IEEE 1364 language
// having, as described in chapter 8 of FinSim's User's Guide
// available www.fintronic.com (click on Support, FAQ, download
// FinSim's Users Guide), the following features:
//  a) additional system functions, such as $VpSin, $VpCos, $VpTan,
//    $VpCtan, $VpAsin, $VpAcos, $VpAtan, $VpActan, $VpPow, $VpPow2,
//    $VpLog, $VpLn, $VpAbs, $VpFloor, $VpHypot, $Fft, $Ifft, etc.
//  b) the capability of declaring variable precision objects by using
//    the attributes varprec = data and varprec = descriptor using the
//    standard Verilog syntax.
//    Such objects can store values in fixed and floating point formats
//    and the size of the formats can be changed during the execution of
//    the simulation, by using $VpAssociateDescriptorToData and
//    $VpSetDescriptorInfo.
//  c) the capability of overloading arithmetic and logical operators
//    as well as assignment operators, i.e. the compiler will find, if
//    possible, the appropriate functions to implement the described
//    functionality, without needing conversion functions.
//  d) support for displaying variable precision objects (%y, %k, %h).
//  e) predefined types VpComplex, VpPolar, VpFComplex, VpFPolar with
//     predefined arithmetic operations
*************************************************************************/
module top;
`include "finsimmath.h"
parameter SIZE = 32 * 32;
parameter ORDER= 6;

VpComplex in[0:ORDER+SIZE-2], out[0:ORDER+SIZE-2], idealOut[0:ORDER+SIZE-2];
VpPolar in_polar[0:SIZE-1], polar_s;
VpReg [0:511] tmp;
VpReg [0:511] a [0:ORDER-1];
VpReg [0:511] b [0:ORDER-1];
VpDescriptor d1;
VpReg [0:511] t[0:ORDER-1];
VpReg [0:511] s[0:ORDER-3];

real acceptableDistance;
integer notDone, j, k, sizeInt, sizeDec, format;
real delta;
real distance;

initial begin
  $VpAssocDescrToData(s, d1);
  $VpAssocDescrToData(t, d1);
  $VpAssocDescrToData(in_polar, d1);
  $VpAssocDescrToData(polar_s, d1);
  $VpAssocDescrToData(a, d1);
  $VpAssocDescrToData(b, d1);
  $VpAssocDescrToData(in, d1);
  $VpAssocDescrToData(out, d1);
  $VpAssocDescrToData(idealOut, d1);

/*************************************************************************
// Because there is a distorsion in phase due to a delay between
// the filtered output and the ideal output, the distance depends on
// the number of samples per time unit, becoming smaller with more samples.
**************************************************************************/
  if (SIZE == 1024) acceptableDistance = 0.032;
  else if (SIZE == 4096) acceptableDistance = 0.012;
  else
    begin
     $display(" acceptableDistance is not yet known for SIZE=%d. Use operands of type real to determine accptableDistance \n", SIZE);
    end

  for (format = 0; format < 2; format = format + 1)
  begin
    if (format == 0) 
      begin
        $display("Try Floating point representation\n");
        sizeInt = 7;
        sizeDec = 14;
      end
    else 
      begin
        $display("Try Two's complement representation\n");
        sizeInt = 7;
        sizeDec = 14;
      end
    notDone = 1;
    while (notDone)
      begin
        if (format == 0) 
          begin
            $VpSetDefaultOptions(sizeInt, sizeDec, `FLOATING,
                  `TO_NEAREST_INTEGER_IF_TIE_TO_MINUS_INF,
                  `SATURATION, 1);
            $VpSetDescriptorInfo(d1, sizeInt, sizeDec, `FLOATING,
                  `TO_NEAREST_INTEGER_IF_TIE_TO_MINUS_INF,
                  `SATURATION, 1);
          end
        else 
          begin
            $VpSetDefaultOptions(sizeInt, sizeDec, `TWOS_COMPLEMENT,
                  `TO_NEAREST_INTEGER_IF_TIE_TO_MINUS_INF,
                  `SATURATION, 1);
            $VpSetDescriptorInfo(d1, sizeInt, sizeDec, `TWOS_COMPLEMENT,
                  `TO_NEAREST_INTEGER_IF_TIE_TO_MINUS_INF,
                  `SATURATION, 1);
          end
        $display("Trying sizeInt=%d, sizeDec=%d\n", sizeInt, sizeDec);

/******************************************************************************
//1. Load input in in.Re and load ideal output in idealOut.Re                  
// Notation:
//   a) sampling_rate: time passed between loading consecutive values in in.Re.  
//   b) SIZE: number of samples
//   c) delta: 2*Pi/SIZE is a constant chosen such that values $VpSin(n*delta*j)
//          with 0 <= j < SIZE represents a sinusoid as function of time
//          with frequency freq = ((sampling_rate/2)/SIZE) * n, in
//          other words within the time span of the collection of
//          all SIZE samples, there are n complete periods of the sinusoid.
********************************************************************************/

        delta = (2*$Pi) / SIZE;
        for (j = ORDER-1; j < ORDER+SIZE-1; j = j + 1)
          begin
            in[j].Re = $VpSin(delta * j) + $VpSin((SIZE/4)*delta*j)/10.0;
            idealOut[j].Re = $VpSin(delta * j);
            in[j].Im = 0.0;
            idealOut[j].Im = 0.0;
          end

/********************************************************************
//2. load coeficients of Butterworth IIR LP filter with
//  passband: 0-500 Hz (assuming a sample rate of 8000 samples /sec)
//  effective order= 5
*********************************************************************/
        $display("Loading coeficients\n");
        a[0] = 0.00016411125;     b[0] = 1.0;
        a[1] = 0.0008205562;      b[1] = -3.7314737;
        a[2] = 0.0016411124;      b[2] = 5.693888;
        a[3] = 0.0016411124;      b[3] = -4.420512;
        a[4] = 0.0008205562;      b[4] = 1.7411026;
        a[5] = 0.00016411125;     b[5] = -0.277753;

/**********************************************************************
//3. Initialize ORDER number of values of the history of 
//   in and out for the filter to operate in best conditions.
//   The values are chosen to be zero in this case. Other values may
//   be better in other circumstances.
***********************************************************************/
        $display("Initialize History\n");
        for (j = 0; j < ORDER; j = j + 1)
          begin
            in[j].Re = 0.0;
            out[j].Re = 0.0;
            idealOut[j].Re = 0.0;
            in[j].Im = 0.0;
            out[j].Im = 0.0;
            idealOut[j].Im = 0.0;
          end


/************************************************************
//4. Perform filtering according to the specified coeficients 
//   and initial values of histories of in and out
************************************************************/ 
       $InitM(out, 0.0, 0.0);
       $display("Perform Filtering\n");
       for (k=ORDER-1; k < ORDER+SIZE-1; k = k+1)
         begin
           t[ORDER-1] = a[ORDER-1] * in[k-ORDER+1].Re;
           for (j = ORDER-2; j >= 0; j = j - 1)
             t[j] = a[j] * in[k-j].Re + t[j+1];
           s[ORDER-3] = -b[ORDER-1]* out[k-ORDER+1].Re - 
                         b[ORDER-2] * out[k-ORDER+2].Re;
           for (j = ORDER-4; j >= 0; j = j - 1)
             s[j] = s[j+1] - b[j+1] * out[k-j-1].Re;
           out[k].Re = t[0] + s[0];
         end

       for (j = ORDER-1; j < ORDER+SIZE-1; j = j + 1)
         $display("filtered output[%d]=%y, idealOut[%d]=%y\n",
               j, out[j].Re,
               j, idealOut[j].Re);

       $display("Compute Distance\n");
       distance = $VpDistAbsSum(out, idealOut);
       distance = distance/SIZE;
       $display("distance between filtered out and ideal output = %e\n", distance);
         if (distance > acceptableDistance)
           begin
             $display("For sizeDec = %d the distance is %e, while acceptable is %e\n",
                      sizeDec, distance, acceptableDistance);
             sizeDec = sizeDec + 1;
          end
        else
          begin
            $display("sizeInt = %d\n sizeDec = %d\n lead to a distance of %e <= acceptable distance of %e", sizeInt,sizeDec, distance, acceptableDistance);
            notDone = 0;
          end
      end
    end
  end
endmodule