/************************************************************************
// Title of Application: Butterworth IIR LP filter using VpFComplex 
//                       operands
//
// This example uses FinMath's VpFComplex and VpFPolar 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.
//
// VpFComplex and VpFPolar are types predefined in FinMath, to contain
// pairs of real values represented on 64 bits. Assignments and
// arithmetic operators (including power) 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.
//
// FinMath also supports VpComplex and VpPolar types, which contain pairs
// of variable precison registers, where the format, size and other
// descriptor information can be associated to the pair at runtime.
//
*************************************************************************/
module top;

parameter SIZE = 32 * 32;
parameter ORDER= 6;

VpFComplex in[0:ORDER+SIZE-2], out[0:ORDER+SIZE-2], idealOut[0:ORDER+SIZE-2];
VpFPolar in_polar[0:SIZE-1], polar_s;

real a [0:ORDER-1];
real b [0:ORDER-1];
real t[0:ORDER-1];
real s[0:ORDER-3];

integer i, j, k;
real delta, distance;

initial begin

/******************************************************************************
//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;
       idealOut[j].Re = $VpSin(delta * j);
       in[j].Im = 0;
       idealOut[j].Im = 0;
     end

/**********************************************************************
//2. 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

/****************************************************************************** 
//3. Load coeficients of Butterworth IIR LP filter with
//      passband: 0-500 Hz (assuming a sample rate of 8000 samples /sec)
//      effective order= 5
*******************************************************************************/
a[0] = 1.6411125E-4;            b[0] = 1.0;
a[1] = 8.205562E-4;             b[1] = -3.7314737;
a[2] = 0.0016411124;            b[2] = 5.693888;
a[3] = 0.0016411124;            b[3] = -4.420512;
a[4] = 8.205562E-4;             b[4] = 1.7411026;
a[5] = 1.6411125E-4;            b[5] = -0.277753;


/************************************************************
//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
/**************************************
//5. Display sampled values - in[].Re 
**************************************/
   for (j = ORDER-1; j < ORDER+SIZE-1; j = j + 1)
     begin
       $display("sampled value[j]=%e\n", in[j].Re);
     end

/*******************************************
//6. Display ideal output values - idealOut 
********************************************/
   for (j = ORDER-1; j < ORDER+SIZE-1; j = j + 1)
     begin
       $display("ideal output[%d]=%e\n", j, idealOut[j].Re);
     end

/***********************************
//7. Display filtered values - out 
***********************************/
   for (j = ORDER-1; j < ORDER+SIZE-1; j = j + 1)
     begin
       $display("filtered output[%d]=%e\n", j, out[j].Re);
     end

/**********************************************************************
//8. Compute distance between filtered output and ideal output vectors  
***********************************************************************/
   $display("Compute Distance\n");
   distance = $VpDistAbsSum(out, idealOut);
   distance = distance/SIZE;
   $display("distance between filtered out and ideal out = %e\n", distance);

end /*initial*/
endmodule