/************************************************************************
// 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[-ORDER+1:SIZE-1], out[-ORDER+1:SIZE-1], idealOut[0:SIZE-1];
VpFPolar in_polar[0:SIZE-1], polar_s;
real a [0:ORDER-1];
real b [0:ORDER-1];
real t[0:ORDER-1], s[0:ORDER-3];
real delta;
integer i, j, k;
real distance;
/**************************************************
//
//Compute distance between vectors out and idealOut
//
//
**************************************************/
task computeDistance;
begin
distance = 0;
for (j = 0; j < SIZE; j = j + 1)
begin
distance = distance + $VpAbs(out[j].Re - idealOut[j].Re);
end
distance = distance / SIZE;
end
endtask
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*$VpGetPi()) / SIZE;
for (j = 0; j < SIZE; 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.
***********************************************************************/
for (j = 0; j < ORDER; j = j + 1)
begin
in[j-ORDER+1].Re = 0;
out[j-ORDER+1].Re = 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
************************************************************/
for (k=0; k < SIZE; k = k+1)
begin
t[ORDER-1] = a[ORDER-1] * in[k-ORDER+1].Re;
for (j = ORDER-2; j >= 0; j = j - 1)
begin
t[j] = a[j] * in[k-j].Re + t[j+1];
end
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)
begin
s[j] = s[j+1] - b[j+1] * out[k-j-1].Re;
end
out[k].Re = t[0] + s[0];
end
/**************************************
//5. Display sampled values - in[].Re
**************************************/
for (j = 0; j < SIZE; j = j + 1)
begin
$display("sampled value[j]=%e\n", in[j].Re);
end
/*******************************************
//6. Display ideal output values - idealOut
********************************************/
for (j = 0; j < SIZE; j = j + 1)
begin
$display("ideal output[%d]=%e\n", j, idealOut[j].Re);
end
/***********************************
//7. Display filtered values - out
***********************************/
for (j = 0; j < SIZE; 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
***********************************************************************/
computeDistance;
$display("distance between filtered out and ideal out = %e\n", distance);
/**********************************************************
//9.Display frequency spectrum of input/sampled values - in
***********************************************************/
$VpFft(in, 0, SIZE-1);
for (j = 0; j < SIZE/2; j = j + 1)
begin
$display("in[%d] Re=%e, Im=%e\n", j, in[j].Re, in[j].Im);
end
$display("finished display of freq dom of input\n");
/********************************************************************
//10.a Display magnitude and phase of spectrum of input/sampled values
// using array assignment with implicit conversion from carthesian
// to polar coordinates
*********************************************************************/
in_polar = in;
for (j = 0; j < SIZE/2; j = j + 1)
begin
$display("in_polar[%d] Mag=%e, Ang=%e\n",
j, in_polar[j].Mag, in_polar[j].Ang);
end
$display("finished display of frequency spectrum of input\n");
/*****************************************************************************
//10.b Display magnitude and phase of spectrum of input/sampled values
// using assignment to scalar with implicit conversion from carthesian to
// polar coordinates
******************************************************************************/
for (j = 0; j < SIZE/2; j = j + 1)
begin
polar_s = in[j];
$display("polar_s[%d] Mag=%e, Ang=%e\n", j, polar_s.Mag, polar_s.Ang);
end
$display("finished display of frequency spectrum of input\n");
/*************************************************************************
//11. Perform $VpIfft on the content of in vector. The result must be
// close to the sampled input. This is just a check for the
// accuracy of $VpFft and $VpIfft system tasks for the given number
// of bits used (precision).
**************************************************************************/
$VpIfft(in, 0, SIZE-1);
for (j = 0; j < SIZE/2; j = j + 1)
begin
$display("should be close to in[%d] Re=%e, Im=%e\n",
j, in[j].Re, in[j].Im);
end
/***********************************************************
//12. Display frequency spectrum of ideal output - idealOut
***********************************************************/
$VpFft(idealOut, 0, SIZE-1);
for (j = 0; j < SIZE/2; j = j + 1)
begin
$display("idealOut[%d] Re=%e, Im=%e\n",
j, idealOut[j].Re, idealOut[j].Im);
end
$display("finished display of freq dom of ideal\n");
/*******************************************************
//13. Display frequency spectrum of filtered output - out
********************************************************/
$VpFft(out, 0, SIZE-1);
for (j = 0; j < SIZE/2; j = j + 1)
begin
$display("out[%d] Re=%e, Im=%e\n", j, out[j].Re, out[j].Im);
end
end /*initial*/
endmodule