发几个用CORDIC快速求三角函数的定点程序

只看该作者 · 2009-12-27 21:53

最近在做FFT，发几个用到的三角函数计算程序

#define DSP32_Q(x) ((int) ((x)*(((unsigned) (1 << (20))))))

const int atanLUT[20] = {
DSP32_Q(0.78539816340),
DSP32_Q(0.46364760900),
DSP32_Q(0.24497866313),
DSP32_Q(0.12435499455),
DSP32_Q(0.06241881000),
DSP32_Q(0.03123983343),
DSP32_Q(0.01562372862),
DSP32_Q(0.00781234106),
DSP32_Q(0.00390623013),
DSP32_Q(0.00195312252),
DSP32_Q(0.00097656219),
DSP32_Q(0.00048828121),
DSP32_Q(0.00024414062),
DSP32_Q(0.00012207031),
DSP32_Q(0.00006103516),
DSP32_Q(0.00003051758),
DSP32_Q(0.00001525879),
DSP32_Q(0.00000762939),
DSP32_Q(0.00000381470),
DSP32_Q(0.00000190735)
};

int MySin(int zi)
{
  int x, y, z, x_new;
  int i, sign;

  x = DSP32_Q(0.60725293500888);
  y = 0;
  z = zi;

  for(i=0; i<20; i++)
  {
if(z >= 0)
{
   sign = 1;
}
else
{
   sign = -1;
}

x_new = x - ((sign*y)>>i);
y = y + ((sign*x)>>i);
x = x_new;
z = z - sign*atanLUT[i];
  }

  return (y);
}

int MyCos(int zi)
{
  int x, y, z, x_new;
  int i, sign;

  x = DSP32_Q(0.60725293500888);
  y = 0;
  z = zi;

  for(i=0; i<20; i++)
  {
if(z >= 0)
{
   sign = 1;
}
else
{
   sign = -1;
}

x_new = x - ((sign*y)>>i);
y = y + ((sign*x)>>i);
x = x_new;
z = z - sign*atanLUT[i];
  }

  return (x);
}

int MyArctan2(int xi, int yi)
{
  int x, y, z, x_new;
  int i, sign;

  x = xi;
  y = yi;
  z = 0;

  for(i=0; i<20; i++)
  {
if(y >= 0)
{
   sign = -1;
}
else
{
   sign = 1;
}

x_new = x - ((sign*y)>>i);
y = y + ((sign*x)>>i);
x = x_new;
z = z - sign*atanLUT[i];
  }

  return (z);
}