icos.c

CHP 3 - Real-Time Digital Signal Processing: Implementations and Applications, Second Edition by Sen

// 
//  Project: Experiment 3.6.5.2 Function Approximaiton - Chapter 3 
//  File name: icos.c   
//
//  Description: This experiment program uses fixed-point C
//               to calculate cosine function in the range of [0, PI/2]
//               by approcimation equation: cos(x)=1-(1/2!)x^2+(1/4!)x^4-(1/6!)x^6
//
//  For the book "Real Time Digital Signal Processing: 
//                Implementation and Application, 2nd Ed"
//                By Sen M. Kuo, Bob H. Lee, and Wenshun Tian
//                Publisher: John Wiley and Sons, Ltd
//
//  Tools used: CCS v.2.12.07
//              TMS320VC5510 DSK Rev-C
//
//

#define PI (double)3.1415926
#define UNITQ15  0x7FFF
#define UNITQ14  0x3FFF

// Function approximation coefficients for cosine
short icosCoef[4]={
  (short)(UNITQ15), 
  (short)(-(UNITQ15/2.0)), 
  (short)(UNITQ15/(2.0*3.0*4.0)),
  (short)(-(UNITQ15/(2.0*3.0*4.0*5.0*6.0)))
};


// Fixed-point implementation of function approximation
short iCos1(short x)
{
  long  cosine,z;
  short x2;

  z = (long)x * x;
  x2 = (short)(z>>15);                       // x2 has x(Q14)*x(Q14)

  cosine = (long)icosCoef[3] * x2;
  cosine = cosine >> 13;                     // Scale back to Q15
  cosine = (short)(cosine + (long)icosCoef[2]) * (long)x2;
  cosine = cosine >> 13;                     // Scale back to Q15
  cosine = (cosine + (long)icosCoef[1]) * x2;
  cosine = cosine >> 13;                     // Scale back to Q15
  cosine = cosine + icosCoef[0];
  return((short)cosine);
}

// Simulated assembly implementation of function approximation
short iCos(short T0)
{
  long  AC0;
  short *ptr;

  ptr = &icosCoef[3];
  AC0 = (long)T0 * T0;
  T0  = (short)(AC0>>15);                  // AC0 has T0(Q14)*T0(Q14)

  AC0 = (long)T0 * *ptr--;
  AC0 = AC0 >> 13;                         // Scale back to Q15
  AC0 = (short)(AC0 + *ptr--) * (long)T0;
  AC0 = AC0 >> 13;                         // Scale back to Q15
  AC0 = (short)(AC0 + *ptr--) * (long)T0;
  AC0 = AC0 >> 13;                         // Scale back to Q15
  AC0 = AC0 +  *ptr;
  return((short)AC0);
}

void main(void)
{
  static double y[5];
  short  x[5],z[5];
	
  x[0] = (short)(0);                     // Q14 format input in rad
  x[1] = (short)(UNITQ14*PI/6.0);        // Q14 format input in rad
  x[2] = (short)(UNITQ14*PI/4.0);        // Q14 format input in rad
  x[3] = (short)(UNITQ14*PI/3.0);        // Q14 format input in rad
  x[4] = (short)(UNITQ14*PI/2.0);        // Q14 format input in rad
	
  z[0] = iCos1(x[0]);                    // Floating-point calculation
  z[1] = iCos1(x[1]);                    // Floating-point calculation
  z[2] = iCos1(x[2]);                    // Floating-point calculation
  z[3] = iCos1(x[3]);                    // Floating-point calculation
  z[4] = iCos1(x[4]);                    // Floating-point calculation

  y[0] = (double)z[0]/UNITQ15;
  y[1] = (double)z[1]/UNITQ15;
  y[2] = (double)z[2]/UNITQ15;
  y[3] = (double)z[3]/UNITQ15;
  y[4] = (double)z[4]/UNITQ15;  
    
  z[0] = iCos(x[0]);                     // Floating-point calculation
  z[1] = iCos(x[1]);                     // Floating-point calculation
  z[2] = iCos(x[2]);                     // Floating-point calculation
  z[3] = iCos(x[3]);                     // Floating-point calculation
  z[4] = iCos(x[4]);                     // Floating-point calculation

  y[0] = (double)z[0]/UNITQ15;
  y[1] = (double)z[1]/UNITQ15;
  y[2] = (double)z[2]/UNITQ15;
  y[3] = (double)z[3]/UNITQ15;
  y[4] = (double)z[4]/UNITQ15;                       
}