系列文章目录

---

一、什么是CORDIC算法？

  CORDIC（COordinate Rotation DIgital Computer）算法是一种用于计算三角函数、双曲函数、对数、指数和平方根等数学函数的迭代算法。对于硬件电路来说，计算三角函数不像CPU那样方便，如果使用大量的乘法器和加法器，那么整个硬件电路会消耗很多资源，因此cordic算法就有效地解决了乘法器和加法器的问题，它的主要优点是：可以使用简单的加法和位移操作，而不需要乘法和除法，从而适合在硬件中实现三角函数的计算，特别是数字信号处理（DSP）和嵌入式系统中。

二、CORDIC算法原理推导

  CORDIC算法的基本原理是：将旋转角度$\theta$分成多个连续的小偏转角，通过逐次摇摆逼近目标旋转角度来完成旋转过程。

  如上图所示：在直角坐标系内将点（$x_1，y_1$）逆时针旋转$\theta$至点（$x_2$，$y_2$），那么怎么计算$cos\theta$和$sin\theta$的值呢？
  我们知道：
$$\{\begin{array}{rl}{x}_1& =\cos \phi \\ y_1& =\sin \phi \end{array}$$

$$\{\begin{array}{rl}{x}_2& =\cos (\theta +\phi )\\ y_2& =\sin (\theta +\phi )\end{array}$$

  根据三角函数和差公式，我们可以得到：
$$\{\begin{array}{rl}{x}_2& =\cos \theta cos\phi -sin\theta sin\phi \\ y_2& =sin\theta cos\phi +cos\theta sin\phi \end{array}$$

  再用$x_1$替换掉$cos\phi$，$y_1$替换掉$sin\phi$可以得到：
$$\{\begin{array}{rl}{x}_2& =cos\theta \ast x_1-sin\theta \ast y_1\\ y_2& =sin\theta \ast x_1+cos\theta \ast y_1\end{array}$$

  也可以写成矩阵形式：

$$\begin{array}{c}\left[\begin{array}{c}{x}_2\\ y_2\end{array}\right]=\left[\begin{array}{cc}cos\theta & -sin\theta \\ sin\theta & cos\theta \end{array}\right]\left[\begin{array}{c}{x}_1\\ y_1\end{array}\right]\end{array}$$
  提取公因子$cos\theta$：
$$\begin{array}{c}\left[\begin{array}{c}{x}_2\\ y_2\end{array}\right]=cos\theta \left[\begin{array}{cc}1& -tan\theta \\ tan\theta & 1\end{array}\right]\left[\begin{array}{c}{x}_1\\ y_1\end{array}\right]\end{array}$$

$$\{\begin{array}{rl}{x}_2& =cos\theta (x_1-tan\theta \ast y_1)\\ y_2& =cos\theta (tan\theta \ast x_1+y_1)\end{array}$$

  这个公式就能将点（$x_1，y_1$）逆时针旋转$\theta$至点（$x_2$，$y_2$）。因为旋转的θ是个固定值，因此$cos\theta$是一个常数，我们暂时去掉$cos\theta$可以得到伪旋转方程：
$$\{\begin{array}{rl}{x}_2& =x_1-tan\theta \ast y_1\\ y_2& =y_1+tan\theta \ast x_1\end{array}$$

  伪旋转方程意思就是，通过伪旋转方程我们可以将点（$x_1，y_1$）逆时针旋转$\theta$角，但是模长不是我们想要的模长。简而言之就是伪旋转方程负责旋转到正确的角度，去掉的$cos\theta$能使得模长为想要的长度，示意图如下所示：

  因此经过伪旋转后到了点$(\widehat{x_2},\widehat{y_2})$，向量R的模值增加了$1/cos\theta$倍。即：向量旋转到了正确的角度，但是模值不对。CORDIC的核心思想就是：先旋转到正确的角度，最后再对模长进行补偿。 我们前面也说了硬件电路不擅长乘除法运算，但是擅长加减法运算。这个伪旋转方程里还是有乘法以及$tan\theta$，我们就要做一种处理:
  如果$tan\theta$刚好等于1/2，即θ角=45°，那么这个伪旋转方程就变成了
$$\{\begin{array}{rl}{x}_2& =x_1-y_1/2\\ y_2& =y_1+x_1/2\end{array}$$
  在硬件电路里，÷2可以用右移一位表示。这样整个方程就变成了只有加减法和移位操作了，因此，我们可以让θ为一系列固定的角度，使得$tan{\theta }^i$ = $2^{-i}$，这样经过一系列旋转后，总能找到一个角度无限逼近与我们想要的角度，这就是CORDIC算法的最核心的思想，故伪旋转方程变为：
$$\{\begin{array}{rl}\widehat{x_2}& =x_1-y_1\ast tan\theta =x_1-y_1\ast 2^{-i}\\ \widehat{y_2}& =y_1+x_1\ast tan\theta =y_1+x_1\ast 2^{-i}\end{array}$$

  其中$tan{\theta }^i$对应的角度值表如下所示：

  这里把旋转角度变换改成了迭代运算。将各种可能的旋转角度进行一些条件限制，可以使得对任意旋转角度都能够通过连续的小角度的迭代来完成。旋转迭代顺序为：

1. 第一次迭代旋转的角度是45.0°
2. 第二次迭代旋转的角度是26.555°
3. 第三次迭代旋转的角度是14.036°
4. 第四次迭代旋转的角度是7.125°
5. 第五次迭代旋转的角度是3.576°
6. 第六次迭代旋转的角度是1.789°
     以此类推，迭代次数越多，精度越高，越能逼近想要的角度θ。在迭代到一定次数后，我们角度已经旋转到正确位置了，现在回头来看我们舍掉的$cos\theta$，这是确定我们模长的因子，由上表也能看出，当θ值越小时，$cos\theta$无限接近1。所以我们迭代过程中舍掉的$$cos{\theta }^1\ast cos{\theta }^2\ast cos{\theta }^3\ast cos{\theta }^4\ast cos{\theta }^5\ast cos{\theta }^6\ast cos{\theta }^7\ast cos{\theta }^8........一定是收敛某个数的$$
     经过计算，在迭代16次以后，上面的乘积≈0.60725941；1/0.60725941 ≈1.64676024187 。;在进行伪旋转之后，需要对伪旋转后的幅度进行标定，所以伪旋转后的坐标还需要乘以一个系数0.60725941才能得到我们想要的正确坐标。

  最终CORDIC算法伪旋转的公式为：
$$\{\begin{array}{rl}{x}^{i+1}& =x^i-d_i\ast y^i\ast \ast 2^{-i}\\ y^{i+1}& =y^i+d_i\ast x^i\ast \ast 2^{-i}\\ z^{i+1}& =z^i-d_i{\theta }^i\end{array}$$
  其中$d_i$表示旋转因子，$d_i$=±1，+1表示逆时针旋转，-1表示顺时针旋转。因为伪旋转不是最终的旋转，还有前面的$cos\theta$的累乘称为缩放因子，伪旋转后的坐标伸缩了目标坐标的Kn倍，其中Kn：

  当确定迭代次数n时，可以预先计算出缩放因子Kn。当n趋近于正无穷时，Kn趋近于1.64676024187，1/Kn趋近于0.60725941。前面通过省略$cos\theta$来进行伪旋转，使得在n次迭代后得到
$$\{\begin{array}{rl}{x}_{伪}& =K_n(x_{0}cos\theta -y_{0}sin\theta )，x_n=x_{伪}\ast 0.60725941\\ y_{伪}& =K_n(y_{0}cos\theta +x_{0}sin\theta )，y_n=y_{伪}\ast 0.60725941\\ z_n& =0\end{array}$$

  也可以设置$x_0=\mathrm{0.0.60725941}，y_0=0$，来消除伸缩因子，从而得到：
$$\{\begin{array}{r}{x}_n=cos\theta \\ y_n=sin\theta \end{array}$$

三、CORDIC模式

3.1 旋转模式

  旋转模式就是，给定一个角度θ和初始坐标，算出旋转后的坐标。旋转模式下，旋转方向是以剩余角度$z^{i+1}$的符号作为每次迭代旋转方向的依据。在旋转模式下，等迭代完成后，$z^{i+1}趋近于0$$，x^{i+1}=cos\theta$，$y^{i+1}=sin\theta$。
  例如，输入向量为$(x_0，y_0)$输入角度为$z$，经过迭代旋转后，剩余角度为0，则输出目标向量为$(xcosz-ysinz)/K_n$，$(ycosz+xsinz)/K_n$其中nK为缩放因子。当旋转角度为30°时，其迭代过程如下:

  $sin30°=0.5006，cos30°=0.8657$

3.2 向量模式

  向量模式就是给一个坐标(x，y)，让向量旋转到x轴上，当y趋近于0时候，向量就无限接近x轴，此时x的坐标就是向量的模长，θ等于arctan(y/x) 。CORDIC算法的向量模式可以计算出输入向量的幅度值。并且使用向量模式就是使旋转后的向量与x轴重合。因此，向量的幅度就是旋转向量的x值，幅度结果由Kn增益标定。向量模式主要是用来求反正切、反正弦、反余弦以及向量幅值等

四、Verilog实现CORDIC

4.1 判断象限

  由迭代公式可以看出，每次迭代的范围是±（45 + 26.6 + 14 +… ）≈ ±99.8°，因此我们要将输入的象限通过三角变换到第一和第四象限，例如
$$cos150°=-cos30°，sin150°=sin30°$$
  所以首先要对输入的角度值判断：

always @(*) beginif((w_act == 1'b1) && (i_angle >=32'd0) && (i_angle < 32'd900000))      //[0:90]第一象限w_quadrant <= 2'd0;else if((w_act == 1'b1) && (i_angle >=32'd900000) && (i_angle <= 32'd1800000)) //[90:180]第二象限w_quadrant <= 2'd1;else if((w_act == 1'b1) && (i_angle >=-32'd1800000) && (i_angle < -32'd900000))[-180:-90]第三象限w_quadrant <= 2'd2;else if(w_act == 1'b1)  [-90:0]第四象限w_quadrant <= 2'd3;
end 

  然后再根据象限对角度进行处理，±180°：

		if(w_quadrant == 2'd1)  //第二象限就减180°r_zi <= i_angle - 32'd1800000;else if(w_quadrant == 2'd2) //第三象限就+180°r_zi <= i_angle + 32'd1800000;elser_zi <= i_angle;

  最后根据象限判断输出的正负值

assign  o_sin   =   ((w_quadrant == 2'd1)||(w_quadrant == 2'd2)) ? -ro_sin :  ro_sin;//根据象限输出正确的数值
assign  o_cos   =   ((w_quadrant == 2'd2)||(w_quadrant == 2'd3)) ? -ro_cos :  ro_cos;//根据象限输出正确的数值

4.2 定义角度表

  由上面公式可以知道，我们必须提前将45°，26.565°…等等存储起来，然后每次按照顺序旋转角度：

//角度扩大了10000倍
wire signed [31:0]  w_arctan[0:15]  ;  //角度值表assign  w_arctan[0 ] = 32'd450000   ;   //arctan45°     = 1
assign  w_arctan[1 ] = 32'd265650   ;   //arctan26.565° = 1/2
assign  w_arctan[2 ] = 32'd140362   ;   //arctan14.036° = 1/4
assign  w_arctan[3 ] = 32'd71250    ;   //arctan7.1250° = 1/8
assign  w_arctan[4 ] = 32'd35763    ;   //arctan3.5763° = 1/16
assign  w_arctan[5 ] = 32'd17899    ;   //arctan1.7899° = 1/32
assign  w_arctan[6 ] = 32'd8951     ;   //arctan0.8951° = 1/64
assign  w_arctan[7 ] = 32'd4476     ;   //arctan0.4476° = 1/128
assign  w_arctan[8 ] = 32'd2238     ;   //arctan0.2238° = 1/256
assign  w_arctan[9 ] = 32'd1119     ;   //arctan0.1119° = 1/512
assign  w_arctan[10] = 32'd559      ;   //arctan0.0559° = 1/1024
assign  w_arctan[11] = 32'd279      ;   //arctan0.0279° = 1/2048
assign  w_arctan[12] = 32'd139      ;   //arctan0.0139° = 1/4096
assign  w_arctan[13] = 32'd69       ;   //arctan0.0069° = 1/8192    
assign  w_arctan[14] = 32'd34       ;   //arctan0.0034° = 1/16384
assign  w_arctan[15] = 32'd17       ;   //arctan0.0017° = 1/32768    

4.3 迭代公式

  最核心的迭代公式就是通过加减和右移实现：

assign w_xi_1 = (w_di == 1'b1) ? (r_xi - (r_yi >>> r_cal_cnt)) : (r_xi + (r_yi >>> r_cal_cnt));
assign w_yi_1 = (w_di == 1'b1) ? (r_yi + (r_xi >>> r_cal_cnt)) : (r_yi - (r_xi >>> r_cal_cnt));
assign w_zi_1 = (w_di == 1'b1) ? (r_zi - w_arctan[r_cal_cnt])  : (r_zi + w_arctan[r_cal_cnt]);

五、仿真验证

5.1 matlab打印各角度的正余弦值

  我们仿真从-180 ° 一直到180°，步进为5°，我们先来用matlab打印出 -180° 、-175°、-170°、、、、一直到180°的sin值和cos值，matlab代码如下：

% 创建一个从-90到90的角度数组
angles = -180:5:180;% 计算正弦和余弦值
sine_values = sin(deg2rad(angles));  % 将角度转换为弧度
cosine_values = cos(deg2rad(angles));% 打印结果
fprintf('角度\t正弦值\t\t余弦值\n');
fprintf('---------------------------\n');
for i = 1:length(angles)fprintf('%d\t%.4f\t%.4f\n', angles(i), sine_values(i), cosine_values(i));
end

  打印结果如下：

5.2 Verilog仿真结果观察

  先来看-175°，在第三象限，sin（-175°）= -0.0872，cos（-175°）= -0.9961

  与matlab一致。

  再来看-140°，在第三象限，sin（-140°）= -0.6427，cos（-140°）= -0.7659

  与matlab一致。

  再来看-55°，在第四象限，sin（-55°）= -0.8191，cos（-55°）= -0.5737

  与matlab一致。

  再来看45°，在第一象限，sin（45°）= 7069，cos（45°）= 7071

  与matlab一致。其他的角度值也是一致的，至此我们用verilog实现了CORDIC算法，接下来我们调用Xilinx 的 CORDIC IP核来使用验证一下。

六、CORDIC IP核介绍