CORDIC算法基本原理

引言

某些复杂的计算，例如三角函数和除法运算等涉及到大量浮点运算的计算任务，是数字电路天生的瓶颈所在。在某些场景下，可以使用查找表方法或者采用级数展开的方法来实现三角函数等运算功能。但是，这两种方法可能会占用大量的存储资源和硬件乘法计算单元，而想要节省资源，就要以牺牲精度为代价。

相对于前两种方法，CORDIC算法具有很大优势。首先，在计算过程中，它不使用任何的硬件乘法器单元，所涉及的只有移位和累加。然后，对于存储资源的占用，它仅仅需要少量的数据需要预先存储。在实际的数字电路设计中，可以将其设计为流水线方式或者是迭代复用方式，以提高运算速度或者是减少资源占用。

一、矢量旋转公式

CORDIC算法最最基本理论基础，是矢量旋转公式。即矢量$A\left(x,y\right)$顺时针旋转$\theta$之后，得到的矢量$\left(x^{\prime },y^{\prime }\right)$可以表示为 （式1）：
$$\begin{array}{l}{x}^{\prime }=x\cos \theta +y\sin \theta ,\\ y^{\prime }=y\cos \theta -x\sin \theta .\end{array}$$

下面用一张图来表示上述公式的原理。

如上图所示，矢量$A\left(x,y\right)$顺时针旋转$\theta$，等效为矢量坐标系逆时针旋转$\theta$。矢量$A$在旋转后的新的坐标系上的映射，就是我们所要求的$\left(x^{\prime },y^{\prime }\right)$。图中蓝色虚线为辅助线，在辅助线的帮助下，上述公式的原理不言自明。

二、矢量旋转公式的变换

考虑一种情景：矢量$A$不断地逆时针旋转，直至旋转到与横坐标平行。如下图所示：

可得 （式2），
$$\begin{array}{l}{x}_{n+1}=x_n\cos {\theta }_n-y_n\sin {\theta }_n=\cos {\theta }_n\left(x_n-y\tan {\theta }_n\right)\\ y_{n+1}=y_n\cos {\theta }_n+x_n\sin {\theta }_n=\cos {\theta }_n\left(y_n+x_n\tan {\theta }_n\right)\end{array}$$
每旋转一次，上述公式就迭代一次。CORDIC算法的核心就是矢量不断旋转，上述公式多次进行迭代的过程。 可见，此时的迭代运算，需要三角函数和乘法的计算，十分复杂，因此需要对其进行简化。

首先，如果把每次迭代过程中的 $\cos {\theta }_n$忽略 （式3），相当于在向量旋转过程中，对矢量进行了缩放，缩放因子为$\frac{1}{\cos {\theta }_n}$。
$$\begin{array}{l}{x}_{n+1}=x_n-y\tan {\theta }_n\\ y_{n+1}=y_n+x_n\tan {\theta }_n\end{array}$$
此时，迭代过程减少了一部分三角函数运算和乘法运算，但是还不够。

在每次旋转中，限定 ${\theta }_n$取某些特殊的值，使得 （式4）：
$$\begin{array}{l}\tan {\theta }_0=2^0,\\ \tan {\theta }_1=2^{-1},\\ \tan {\theta }_2=2^{-2},\\ \cdots \\ \tan {\theta }_n=2^{-n},\\ \cdots \end{array}$$
那么可以将迭代公式变成下面迭代公式所示 （式5） ：
$$\begin{array}{l}{x}_{n+1}=x_n-d_n\left(y_n\cdot 2^{-n}\right)\\ y_{n+1}=y_n+d_n\left(x_n\cdot 2^{-n}\right)\end{array}$$
可见，在每次迭代过程中，选定某些特殊的旋转角，使用移位运算来代替乘法运算。为了使得矢量可以不断接近目标位置，旋转的角度随着旋转次数的增加而减少的。同时，旋转方向是不断调整的，即如果越过了目标位置，则下一次旋转方向要取反。我们约定，逆时针方向时， $d_n=1$。

通过上述一系列的变换，CORDIC算法所使用的迭代过程，对于硬件实现来说就十分友好了。

三、COEDIC算法的实际应用

下面，我们使用COEDIC算法来求解具体的计算问题。

1. 正弦函数$\sin \left(z\right)$和余弦函数$\cos \left(z\right)$的求解
   可知 （式6），
   $$\begin{array}{l}{x}_{n+1}=x_n-d_n\left(y_n\cdot \tan {\theta }_n\right)=\frac{1}{\cos {\theta }_n}\left[\cos {\theta }_n{x}_n-d_n\left(y_n\cdot \sin {\theta }_n\right)\right]\\ y_{n+1}=y_n+d_n\left(x_n\cdot \tan {\theta }_n\right)=\frac{1}{\cos {\theta }_n}\left[\cos {\theta }_n{y}_n+d_n\left(x_n\cdot \sin {\theta }_n\right)\right]\end{array}$$
   因此有 （式7）：
   $$\begin{array}{l}{x}_{n+1}=\underset{n}{\Pi }\frac{1}{\cos {\theta }_i}\left[\cos \left(\sum _{i=0}^n{d}_i{\theta }_i\right)x_0-\sin \left(\sum _{i=0}^n{d}_i{\theta }_i\right)y_0\right]\\ y_{n+1}=\underset{n}{\Pi }\frac{1}{\cos {\theta }_i}\left[\cos \left(\sum _{i=0}^n{d}_i{\theta }_i\right)y_0+\sin \left(\sum _{i=0}^n{d}_i{\theta }_i\right)x_0\right]\end{array}$$
   为了简单，令$k_n=\underset{n}{\Pi }\frac{1}{\cos {\theta }_i}$。
   设置迭代初始值 $x_0=\frac{1}{k_n}$，$y_n=0$。每次迭代的角度旋转方向取 $d_n=sign\left(z-\sum _{i=0}^{n-1}{d}_i\theta {}_i\right)$，可知，随着迭代的深入，旋转的角度$\sum _{i=0}^n{d}_i\theta {}_i$ 会逐渐的趋近于$z$。当迭代到第$\left(n+1\right)$次时，上述公式可化为 （式8）：
   $$\begin{array}{l}{x}_{n+1}=k_n\left[\cos \left(\sum _{i=0}^n{d}_i{\theta }_i\right)x_0-\sin \left(\sum _{i=0}^n{d}_i{\theta }_i\right)y_0\right]=\cos \left(\sum _{i=0}^n{d}_i{\theta }_i\right)\Rightarrow \cos z\\ y_{n+1}=k_n\left[\cos \left(\sum _{i=0}^n{d}_i{\theta }_i\right)y_0+\sin \left(\sum _{i=0}^n{d}_i{\theta }_i\right)x_0\right]=\sin \left(\sum _{i=0}^n{d}_i{\theta }_i\right)\Rightarrow \sin z\end{array}$$
   在这种初始值预设和旋转方向下，不断进行（式6）的计算，也就是（式5）的计算，最终计算出$\cos z$和$\sin z$的近似值。迭代次数$n$越大，近似的精度越高。
2. 反正切函数$\arctan (y/x)$的求解
   对于反正切函数来说，参照正弦余弦的模式，只需要初始值$\left(x_0,y_0\right)$ 设置为要求的正切参数的$\left(x,y\right)$，将迭代方向取为$d_n=-sign\left(x_n{y}_n\right)$，即不断朝着使得$y_n=0$的方向进行旋转，因此，最终$\sum _{i=0}^n{d}_i{\theta }_i$所累积的角度，即为所求的反正切角度。
3. 除法运算
   除法运算，因为不涉及到三角函数，只需要按照CORDIC算法的思想去求即可。除法运算的CODRDIC算法的迭代公式为：
   $$\begin{array}{l}{x}_{n+1}=x_0\\ y_{n+1}=y_n-d_n{2}^{-n}\cdot x_n\end{array}$$
   迭代方向 $d_n=sign\left(y_{n-1}\right)$，即不断地朝着 $y_{n+1}=0$ 的方向去迭代。
   这个是有问题的吧？如果 $x_0\ll y_0$，由于$\sum _{n=0}^{\infty }{2}^{-n}=2$，上述迭代永远也不会收敛，差距太大了。
   能否解决上述问题？个人认为，可以对n的值做预设，不能从0开始，需要根据 $x_0$ 和 $y_0$来决定n的初始值。根据$\sum _{n=i}^{\infty }{2}^{-n}=2^{1-i}$，可以根据 $x_0$ 和 $y_0$实际有效位数的差异，来决定，例如 $y_0$ 的有效位宽是5，$x_0$ 的有效位宽是2，因此他们的商可能在8左右，此时n的初始值，应该取-2开始。

================================================
参考资料：《通信IC设计》