CORDIC算法详解

一、什么是CORDIC

CORDIC是一个漂亮的算法, 全称Coordinate Rotation Digital Computer, 即“数字坐标轴旋转计算”算法。它可以用于高精度计算各种三角函数和根号，且效率高，使用方便，不需要硬件乘法。

CORDIC出现最早是在20世纪50年代，最初是用来计算需要进行旋转而又不需要特别高的精度的数值。20世纪末，CORDIC又开始被计算机、数字信号处理器、运算放大器等各个领域广泛使用。

计算结果的精度由所采用的迭代次数决定，多迭代几次计算精度就越高，迭代次数可以根据需要取值。

二、CORDIC基本思想

CORDIC算法可以把一个向量递归旋转直到旋转向量的终点在一个期望的位置（当旋转到期望的位置时旋转操作停止并返回）。这个向量的实部和虚部对应于三角函数计算中的余弦和正弦。

CORDIC算法的基本原理是使用数字位移操作，而不是复杂的乘法和除法操作。所以CORDIC算法可以在没有FPU的微处理器中被轻松地实现。

为了把坐标从(x,y)旋转到(x,y+z)，根据旋转公式COS(a-b)=COS(a)COS(b)+SIN(a)SIN(b)，对于初始角度为0的单位坐标向量(x, y)，通过转化(把坐标(x,y)旋转到(x,y+z))，就可以将它转化成一个旋转过z度的新向量(x1,y1)。它的公式是x1=x*Cos(z)-y*Sin(z), y1=x*Sin(z)+y*Cos(z)。每次的旋转角度是要处理的余弦值和正弦值的估计，所以CORDIC在每一次的迭代中增加角度以及调整放大系数用来进一步改善精度。

三、CORDIC数学推导

以下是CORDIC算法的数学推导，包括旋转方程的推导、旋转角度计算公式的推导和相应的三角函数估计值的计算方法。

CORDIC旋转方程推导

x1=y*sinh(z)+x*cosh(z)
y1=x*sinh(z)+y*cosh(z)

CORDIC旋转角度计算公式推导

从x轴正向到y轴正向旋转的角度为pi/2, 从x轴正向到x轴负向旋转的角度为-pi。为了找到导致x趋近于0和y的符号变化的最小的角度z, 我们采用迭代的方法（角度分解法），如下公式所示：

kd = 1
for (k = 0; k  0) {
    x1 = x + y/left_shift(k) /*将y缩小k位，然后除以kd进行余弦运算*/
    y1 = y - x/left_shift(k)
    z = z + atan(left_shift(-k))/kd
  } else {
    x1 = x - y/left_shift(k)
    y1 = y + x/left_shift(k)
    z = z - atan(left_shift(-k))/kd
  }
  x = x1
  y = y1
  kd = kd * sqrt(1 + left_shift(-2*k))
}

CORDIC三角函数的估计值计算方法

在CORDIC迭代的过程中计算cos(theta) 和 sin(theta)，第k次迭代后给出cos(theta) 和 sin(theta) 的估计值:

cos(k) = P(i=0, k-1)(cos(atan(2^-i)))
sin(k) = P(i=0, k-1)( d(i)*sin(atan(2^-i)) )

其中，P()表示π的积分，d(i)表示第i位的符号。

四、示例代码

CORDIC在C语言中的实现

#include
#include

int cordic_ctab [] = { // atan(2^-i)
  0x3243F6A9, 0x1DAC6705, 0x0FADBDDF, 0x07F56EA6,
  0x03FEAB76, 0x01FFD55B, 0x00FFFAAA, 0x007FFF55,
  0x003FFFEA, 0x001FFFFD, 0x000FFFFF, 0x0007FFFF,
  0x0003FFFF, 0x0001FFFF, 0x0000FFFF, 0x00007FFF,
  0x00003FFF, 0x00001FFF, 0x00000FFF, 0x000007FF,
  0x000003FF, 0x000001FF, 0x000000FF, 0x0000007F,
  0x0000003F, 0x0000001F, 0x0000000F, 0x00000007,
  0x00000003, 0x00000001, 0x00000000, 0x00000000
};

void cordic(int *x, int *y, int *z) {
  int i, tx, ty, tz, d;
  tz = 0;
  for (i = 0; i < 32; i++) {
    d = (*y > i);
    ty = (*y+((d * (*x))>>i));
    tz += d * cordic_ctab[i];
    *x = tx, *y = ty;
  }
  *z = tz;
}

void print_cordic(int x, int y, int z, double dt) {
  printf("CORDIC: cos(%f) = %f, sin(%f) = %f\n", dt, (double)x / (double)(1<<15), dt, (double)y / (double)(1<<15));
  printf("         atan2(sin(%f), cos(%f)) = %f\n\n", dt, dt, (double)z / (double)(1<<15));
}

int main() {
  int x, y, z;
  double dt;
  for (dt = 0.0; dt < 2.0 * M_PI; dt += 0.1) {
    x = 1<<15; y = 0;
    cordic(&x, &y, &z);
    print_cordic(x, y, z, dt);
  }
}

CORDIC在Verilog中的实现

module cordic_top(input clk, input rst, input [15:0] theta, output [15:0] x, output [15:0] y, output [15:0] z);
  integer i;
  reg [15:0] x_reg, y_reg, z_reg;
  reg [15:0] temp;
  reg [31:0] kd;

  const [31:0] cordic_ctab [0:31] = { // atan(2^-i)
    16'd3217, 16'd1609, 16'd804, 16'd402,
    16'd201, 16'd100, 16'd50, 16'd25,
    16'd12, 16'd6, 16'd3, 16'd1,
    16'd0, 16'd0, 16'd0, 16'd0,
    16'd0, 16'd0, 16'd0, 16'd0,
    16'd0, 16'd0, 16'd0, 16'd0,
    16'd0, 16'd0, 16'd0, 16'd0,
    16'd0, 16'd0, 16'd0, 16'd0
  };

  always @ (posedge clk) begin
    if (rst) begin
      x_reg <= 16'd16384;
      y_reg <= 16'd0;
      z_reg <= 16'd0;
      kd <= 33'd133997;
    end else begin
      temp <= (kd <= 33'd273296) ? cordic_ctab[31] : cordic_ctab[31-$clog2(kd-$clog2(kd))]; // choose table index
      if (y_reg[15]) begin
        {x_reg, y_reg} <= {y_reg[14:0], x_reg[15]};
        z_reg <= z_reg - temp;
      end else begin
        {x_reg, y_reg} <= {y_reg[14:0], ~(x_reg[15])};
        z_reg <= z_reg + temp;
      end
      kd > $clog2(kd * sqrt(1+(2**(-2*i))));
    end
  end

  assign x = x_reg;
  assign y = y_reg;
  assign z = z_reg;
endmodule