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