用cordic实现的反正切.c

用CORDIC算法实现的2参数反正切。结果的精度与CORDIC的迭代次数有关

//用CORDIC实现的反正切
 

//      用CORDIC算法实现的2参数反正切。结果的精度与CORDIC的迭代次数有关，迭代次数越多，精度越高。本例子中精确到小数点后4位。要提高迭代次数，还得把增加1QN格式的位数，比如32位long，程序多处需要修改，有需要的话自己改吧。
//      测试时间：10M晶振，51单片机(西门子C515C)，库函数里提供的actan2函数1.3~6.9ms,x,y越相近时间越短，各个运算时间因参数不同而不同。cordic_myactan函数2.9ms,不因参数不同而变化。
 //     偶尔一次看到DSP实现乘法的资料，才实现了本函数，定义QN为C的整型，其加减乘除和整型一样。希望对想在单片机里实现CORDIC的朋友有帮助。如有建议或提供更好的实现代码，感激不尽………… [e-mail]。本人也搜集了些CORDIC的资料，需要的话可以找我。
//代码如下：
typedef  short int QN1;
typedef short  int QN2 ;
//1QN 2QN数据格式均16位，1QN第一位位符号位，第二位后小数点,定点数。即0x2000=0.5,0xe000=-0.5
//2QN第一位符号位，第3位后小数点，定点数。即0x2000=1,0x1000=0.5
//1QN和2QN是XILINX CORDIC v2.0说明文档里提到的表示方法，见cordic.pdf
//1QN即相当于定点数里的Q14表示法，2QN相当于Q13,参见《DSP芯片的原理与开发应用(第3版)》(电子工业出版社
//张雄伟等著)第3章 DSP芯片的定点运算
QN2 cordic_myactan(QN1 x, QN1 y);//输入参数为1QN格式，输出结果为2QN格式
#define CORDIC_LENGTH 14 //迭代次数
/*#define PI       3.141592653589793
#define HALF_PI  (PI * 0.5)
double CORDIC_TABLE[CORDIC_LENGTH]={0.785398,0.463648,0.244979,0.124355,0.062419,0.031240,
0.015624,0.007812,0.003906,0.001953,0.000977,0.000488,0.000244,0.000122,0.000061,0.000031,
0.000015,0.000008,0.000004,0.000002,0.000001};*/
//上面的浮点数对应的2QN格式的定点数为
#define PI       0x6487
#define HALF_PI  0x3243
QN2 CORDIC_TABLE[CORDIC_LENGTH]={0x1921,0xED6,0x7D6,0x3FA,0x1FF,0xFF,0x7F,0x3F,0x1F,0xF,
0x8,0x3,0x1,0x0/*,0x0,0x0,0x0,0x0,0x0,0x0,0x0*/};

QN2 cordic_myactan(QN1 x, QN1 y)
{
  int L;
  QN1 tmp_I;
  QN2   phase_rads, acc_phase_rads;
  
  if (x < 0) {
    /* rotate by an initial +/- 90 degrees */
    tmp_I = x;
    if (y > 0) {
       x = y;  /* subtract 90 degrees */
       y = -tmp_I;
       acc_phase_rads = -HALF_PI;
    } else {
       x = -y;  /* add 90 degrees */
       y = tmp_I;
       acc_phase_rads = HALF_PI;
    }
  } else {
    acc_phase_rads = 0;
  }
  
  
  /* rotate using "1 + jK" factors */
  for (L = 0; L <= CORDIC_LENGTH; L++) {
   // K = mp_cordic_table[L].K;
    phase_rads = CORDIC_TABLE[L];
    tmp_I = x;
    if (y >= 0) {
      /* phase is positive: do negative roation */
      x += y>>L;//*pow(2,-L); // K  K=2的-L次方,见CORDIC例子代码;
      y -= tmp_I>>L;// *pow(2,-L);
      acc_phase_rads -= phase_rads;
    } else {
      /* phase is negative: do positive rotation */
      x -= y >>L;//*pow(2,-L);
      y += tmp_I>>L;//*pow(2,-L);
      acc_phase_rads += phase_rads;
    }
  }
  return -acc_phase_rads;
}