Cordic算法

Cordic算法

​ CORDIC 算法是坐标旋转数字计算的缩写，它最初用于三角函数的坐标变换，经过一定的推广后也可用于计算线形函数和双曲线函数(开平方根)。CORDIC算法只由移位操作和加减操作，因此，非常适合于在硬件使用。

---

一、圆周系统

圆周系统下有两种模式，旋转模式和向量模式

1.1 旋转模式

1.1.1 原理

​ 图1.1是一个单位圆。对于P/Q两点坐标为式（1.1） （1.2）

图1.1 旋转模式

$$\{\begin{array}{l}{x}_P=cos\alpha \\ y_P=sin\alpha \end{array}\text{\hspace{1em}(1.1)}$$

$$\{\begin{array}{l}{x}_Q=cos(\alpha +\theta )=cos\alpha cos\theta -sin\alpha sin\theta \\ y_Q=sin(\alpha +\theta )=sin\alpha cos\theta +sin\theta cos\alpha \end{array}\text{\hspace{1em}(1.2)}$$

​ 对于式（1.2），提出cosθ，并带入式（1.1）得到式（1.3）
$$\{\begin{array}{l}{x}_Q=cos\theta (x_P-y_{P}tan\theta )\\ y_Q=cos\theta (y_P+x_{P}tan\theta )\end{array}\text{\hspace{1em}(1.3)}$$
​ 由伪旋转图1.2可知，OP逆时针旋转θ°之后，理应到达OQ位置，根据三角函数关系，可得OQ的长度（模）刚好是OR的cosθ倍,相位一致**。因此R的坐标如式（1.4）。我们称OR是伪旋转，并先计算伪旋转，再补偿模值。

图1.2 伪旋转

$$\{\begin{array}{l}{x}_R=x_P-y_{P}tan\theta \\ y_R=y_P+x_{P}tan\theta \end{array}\text{\hspace{1em}(1.4)}$$

​ 对于伪旋转的运算是算法的核心所在，首先分解成一次次的旋转迭代，先将θ分解如式（1.5）
$$\theta =\sum _{i=0}^{\infty }{\theta }_i\text{\hspace{1em}(1.5)}$$
​ 然后把坐标计算分解，第i+1次的计算公式（1.6）
$$\{\begin{array}{l}{x}_{i+1}=x_i-y_{i}tan{\theta }_i\\ y_{i+1}=y_i+x_{i}tan{\theta }_i\end{array}\text{\hspace{1em}(1.6)}$$
​ 由于tanθ并不是一个容易得到的值，因此把tanθ替换，如式（1.7）.这里又引入了一个z，是干什么用的呢？要理解z，首先理解旋转迭代的过程，因为旋转迭代是可以看成围绕最终目标（这里是OQ）的来回振荡，如果这次转到目标上方（也就是相位超出目标），下次就要往下转（顺时针），反之同理。因此需要记录本次旋转后与目标的距离和相位情况，而这就是z的作用。而d决定旋转方向，如果zi是负数，说明相位落后，要增加相位，因此取负然后负负得正，使zi+1变大，相反，如果zi是正数，说明相位超前，需要相位回来一点，因此取负让zi+1等于zi减去本次微旋转角度θi。
$$\begin{array}{c}tan{\theta }_i=d_i{2}^{-i}\\ z_{i+1}=z_i-{\theta }_i=z_i-ta{n}^{-1}(d_i{2}^{-i})\\ d_i=\{\begin{array}{ll}+1& z_i\ge 0\\ -1& z_i<0\end{array}\end{array}\text{\hspace{1em}(1.7)(1.8)(1.9)}$$
​ 伪旋转迭代完成后，需要进行补偿（也可以在旋转开始前补偿）。令补偿因子为K，则实际旋转迭代方程如（1.10）,由于x、y的次数相同，因此可以把补偿因子提出，结合式（1.7）累乘得到式（1.11）和式（1.12）。
$$\{\begin{array}{l}{x}_{i+1}=K_i(x_i-y_i{d}_i{2}^{-i})\\ y_{i+1}=K_i(y_i+x_i{d}_i{2}^{-i})\end{array}\text{\hspace{1em}(1.10)}$$

$$\begin{array}{c}{K}_i=cos{\theta }_i=\frac{1}{\sqrt{1^2+(2^{-i}{)}^2}}\\ K=\prod _0^{\infty }{K}_i\end{array}\text{\hspace{1em}(1.11)(1.12)}$$
​ 由于cosθi的值是固定的，如图1.3，因此可以提前把K算出来，取0.607252935.

图1.3 K取值表

​ 记输入为x0、y0、z0,分别对应图1.1的xp、yp、角度θ**注意，这里的xn、yn是伪旋转输出结果，因此尚未模值补偿，因此等式（1.13）右边的真实坐标需要乘K的倒数**。经过n次迭代，最终结果可以表示为

$$\begin{array}{c}\{\begin{array}{l}{x}_n=\frac{1}{K}(x_{0}cos{z}_0-y_{0}sin{z}_0)\\ y_n=\frac{1}{K}(y_{0}cos{z}_0+x_{0}sin{z}_0)\\ z_n=0\end{array}\end{array}\text{\hspace{1em}(1.13)}$$

1.1.2 应用

​ 旋转模式可以计算三角函数sin和cos，可以让z0为待求角度θ，y0为0，x0为K，输出为
$$\begin{array}{c}\{\begin{array}{l}{x}_n=cos\theta \\ y_n=sin\theta \\ z_n=0\end{array}\end{array}\text{\hspace{1em}(1.14)}$$

1.2 向量模式

1.2.1 原理

​ 向量模式与旋转模式基本原理相似，都是在单位圆做旋转。区别在于向量模式求得是模长和角度，本质是直角坐标向极坐标的转换。图1.4即为向量模式，只有一个点P，经过n次迭代，最终使OP顺时针旋转到x轴正半轴上，从而得到角度。

图1.4 向量模式

​ 伪旋转公式（1.15），由于目标是趋向x正半轴，因此决定旋转方向d的是y坐标，当y为正数，说明高于x轴，需要顺时针移动，则d为负，反之同理。

$$\{\begin{array}{l}{x}_{i+1}=x_i-y_i{d}_i{2}^{-i}\\ y_{i+1}=y_i+x_i{d}_i{2}^{-i}\\ z_{i+1}=z_i-{\theta }_i=z_i-ta{n}^{-1}(d_i{2}^{-i})\\ d_i=\{\begin{array}{ll}-1& y_i\ge 0\\ +1& y_i<0\end{array}\end{array}\text{\hspace{1em}(1.15)}$$
​ 同样，记输入为x0、y0、z0,分别对应图1.4的xp、yp、角度θ.伪旋转的最终结果为式（1.16），xn对应最终模值，yn对应差值，zn对应旋转角度。这里有个小思考，为什么zn是加θ呢？我的想法是，观察式1.15可知，第一次计算z1的值时，d0为负，因此实际z0是加上θ0（45°）。以此类推，每当高于x轴时（y>0），z是增加的，而OP又是从高于x轴的位置顺时针旋转，因此z最终结果必定是正。

$$\begin{array}{c}\{\begin{array}{l}{x}_n=\frac{1}{K}\sqrt{x_0^2+y_0^2}\\ y_n=0\\ z_n=z_0+\theta =z_0+ta{n}^{-1}(y_0/x_0)\end{array}\end{array}\text{\hspace{1em}(1.16)}$$

1.2.2 应用

• 求复数的模和相位

​ 令x0，y0为分别为复数x+yi的实部和虚部，并且对xn输出模值校正（乘K），结果x’n输出模值，z0输出相位。实际操作时也可把x0和y0先校正。

• 求反正切

  令x=1，y0为角度，则z0输出y0的反正切值。由于不涉及幅度，因此不需要校正。

CORDIC算法详解(一)-CORDIC 算法之圆周系统之旋转模式_cordic碎碎思-CSDN博客

CORDIC算法详解(二)-CORDIC 算法之圆周系统之向量模式_cordic向量模式-CSDN博客

二、双曲系统

2.1 旋转模式

2.1.1 原理

​ 在双曲系统，点的坐标落在双曲线正半轴上，分析方法与圆周系统一样，只不过方程和坐标略有不同。

图2.1 双曲系统

​ 由图得出PQ两点坐标，注意别跟正弦余弦函数搞混 $$ \begin{cases} x_P=cosh\alpha\\ y_P=sinh\alpha\\ \end{cases} \tag{2.1} $$

$$\{\begin{array}{l}{x}_Q=cosh(\alpha +\theta )=cosh\alpha cosh\theta +sinh\alpha sinh\theta \\ y_Q=sinh(\alpha +\theta )=sinh\alpha cosh\theta +sinh\theta cosh\alpha \end{array}\text{\hspace{1em}(2.2)}$$

​ 将式（2.1）带入式（2.2），并提出coshθ，得到式（2.3)
$$\{\begin{array}{l}{x}_Q=cosh\theta (x_P+y_{P}tanh\theta )\\ y_Q=cosh\theta (y_P+x_{P}tanh\theta )\end{array}\text{\hspace{1em}(2.3)}$$
​ 然后就是熟悉的伪旋转分解：
$$\{\begin{array}{l}{x}_{i+1}=x_i+y_i{d}_i{2}^{-i}\\ y_{i+1}=y_i+x_i{d}_i{2}^{-i}\\ z_{i+1}=z_i-{\theta }_i=z_i-tan{h}^{-1}(d_i{2}^{-i})\\ d_i=\{\begin{array}{ll}+1& z_i\ge 0\\ -1& z_i<0\end{array}\end{array}\text{\hspace{1em}(2.4)}$$
​ 当然，由于是cosh而不是cos，因此迭代过程略有不同。当迭代到4、13、40等满足i=3k+1形式，必须重复迭代一次，即1,2,3,4,4,5… 另外模补偿因子也有不同，见式（2.5），K≈1.20749706

$$\begin{array}{c}{K}_i=cosh{\theta }_i=\frac{1}{\sqrt{1^2-(2^{-i}{)}^2}}\\ K=\prod _1^{\infty }{K}_i\end{array}\text{\hspace{1em}(2.5)}$$
​ 由于i=0取不到，因此从1开始迭代，最终迭代结果如式（2.6）

$$\begin{array}{c}\{\begin{array}{l}{x}_n=\frac{1}{K}(x_{1}cosh{z}_1-y_{1}sinh{z}_1)\\ y_n=\frac{1}{K}(y_{1}cosh{z}_1+x_{1}sinh{z}_1)\\ z_n=0\end{array}\end{array}\text{\hspace{1em}(2.6)}$$

2.1.2 应用

​ 旋转模式下可以求双曲余弦函数和双曲正弦函数，令x1=1/K,y1=0,z1=θ,得到
$$\begin{array}{c}\{\begin{array}{l}{x}_n=cosh\theta \\ y_n=sinh\theta \\ z_n=0\end{array}\end{array}\text{\hspace{1em}(2.7)}$$

2.2 向量模式

2.2.1 原理

​ 在双曲系统，点的坐标落在双曲线正半轴上，分析方法与圆周系统一样，只不过方程和坐标略有不同。迭代分析过程不表，直接给出结果式（2.8）
$$\begin{array}{c}\{\begin{array}{l}{x}_n=\frac{1}{K}\sqrt{x_1^2-y_1^2}\\ y_n=0\\ z_n=z_1+tan{h}^{-1}(y_1/x_1)\end{array}\end{array}\text{\hspace{1em}(2.8)}$$

2.2.2 应用

​ 可以求开方和对数，令x1=θ+1，y1=θ-1，z1=0
$$\begin{array}{c}\{\begin{array}{l}{x}_n=\frac{2}{K}\sqrt{\theta }\\ y_n=0\\ z_n=z_1+tan{h}^{-1}(\frac{\theta -1}{\theta +1})=\frac{1}{2}\ln \theta \end{array}\end{array}\text{\hspace{1em}(2.9)}$$