[学习笔记]走进CORDIC算法的世界（理论）

它是一种坐标数字计算的方法，这个方法在1959年被提出，主要用于三角函数、双曲线、指数、对数的计算。

该算法通过基本的加和移位运算代替乘法运算，使得矢量的旋转和定向的计算不再需要三角函数、乘法、开方、反三角、指数等函数。

DDS中需要正弦波的情况较多，相位是线性的，得到的相位然后获得其幅值，常用的做法：提前把sin的值和cos的值存放到ROM表（即mif表）里面，相位作为地址，然后寻址即可得到sin和cos的值。

mif表消耗存储资源大，对于采样深度为16M，宽度为16位的mif表，所需要的空间就是16M28=1Gbits。（大容量存储器限制系统工作频率，拉低成本）。

如图所示，初始向量（X1,Y1）旋转θ角度之后得到向量（X2,Y2），此向量有以下关系： X 2 = X 1 ∗ c o s ( θ ) − Y 1 ∗ s i n ( θ ) X2=X1*cos(θ)-Y1*sin(θ) X2=X1∗cos(θ)−Y1∗sin(θ) Y 2 = Y 1 ∗ c o s ( θ ) + X 1 ∗ s i n ( θ ) Y2=Y1*cos(θ)+X1*sin(θ) Y2=Y1∗cos(θ)+X1∗sin(θ)

注：θ为待求角

上面的方程组同样可以写成矩阵向量的形式：

[ x 2 y 2 ] = [ cos ⁡ θ − sin ⁡ θ sin ⁡ θ cos ⁡ θ ] [ x 1 y 1 ] \left[\begin{array}{l} 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}{l} x_{1} \\ y_{1} \end{array}\right] [x2​y2​​]=[cosθsinθ​−sinθcosθ​][x1​y1​​]

例如一个90°相移：

[ x 2 y 2 ] = [ 0 − 1 1 0 ] [ x 1 y 1 ] \left[\begin{array}{l} x_{2} \\ y_{2} \end{array}\right]=\left[\begin{array}{cc} 0 & -1 \\ 1 & 0 \end{array}\right]\left[\begin{array}{l} x_{1} \\ y_{1} \end{array}\right] [x2​y2​​]=[01​−10​][x1​y1​​]

通过提出因数cos(θ)，方程可以转换成下面的形式：

X 2 = c o s ( θ ) ( X 1 − Y 1 ∗ t a n ( θ ) ) X2=cos(θ)(X1-Y1*tan(θ)) X2=cos(θ)(X1−Y1∗tan(θ))

Y 2 = c o s ( θ ) ( Y 1 + X 1 ∗ t a n ( θ ) ) Y2=cos(θ)(Y1+X1*tan(θ)) Y2=cos(θ)(Y1+X1∗tan(θ))

如果去除cos(θ)项，得到伪旋转方程：

x 2 ^ = x 1 − y 1 t a n θ \hat{x_{2}}=x_{1}-y_{1}tan \theta x2​^​=x1​−y1​tanθ

y 2 ^ = y 1 + x 1 t a n θ \hat{y_{2}}=y_{1}+x_{1}tan \theta y2​^​=y1​+x1​tanθ

即旋转的角度是正确的，但是x与y的值增加 c o s − 1 θ cos^{-1} \theta cos−1θ

故模值变大。

去除cosθ项可以简化坐标平面旋转的计算操作，但是我们并不能通过适当的数学方法去除cosθ项。

经过伪旋转之后，向量R的模值将增加1/cosθ倍，向量旋转了正确的角度，但模值出现了错误。

CORDIC方法的核心是伪装旋转角度θ

假设初始向量经过N次旋转之后得到新的向量，且每次旋转角度θ正切值都为2的幂，则第i次旋转角度为：

t a n θ i = 2 − i tan \theta^{i} = 2 ^{-i} tanθi=2−i

下面的表格指出用于CORDIC算法中每个迭代（i）的旋转角度（精确到9位小数）：

由此表可以注意有三个方面的变化：

由于i是整数，故对应的角度值都是一一确定的，可以通过几个角度的加减组合来达到所需的角度值。

在这里，我们把变换改成了迭代算法。我们将各种可能的旋转角度加以限制，使得对任意角度 的旋转能够通过一系列连续小角度的旋转迭代 i 来完成。旋转角度遵循法则 :

t a n θ i = 2 − i tan \theta^{i} = 2 ^{-i} tanθi=2−i

遵循这样的法则，乘以正切项变成了移位操作。

前几次迭代的形式为 : 第 1 次迭代 : 旋转 45°; 第 2 次迭代 : 旋转 26.6° 第 3 次迭代 : 旋转 14° 等。

连续旋转13次，总的旋转的角度是99.7或者-99.7（假设每次都是角度相加或相减，也就是转的方向相同），每转一个角度a,对应的坐标值变为原来的1/cosa倍,连续旋转13次，那么最终的坐标的值就是原来的1/（cosa1*cosa2……）倍，也就是1.6467602。

即 c o s 45 ⋅ c o s 26.5 ⋅ c o s 14.03 ⋅ c o s 7.125 ⋯ ⋯ c o s 0.0139 = 0.607252941 cos 45 \cdot cos 26.5 \cdot cos 14.03 \cdot cos 7.125 \cdots \cdots cos 0.0139 = 0.607252941 cos45⋅cos26.5⋅cos14.03⋅cos7.125⋯⋯cos0.0139=0.607252941

而 1 0.607252941 = 1.6467602 \frac{1}{0.607252941}= 1.6467602 0.6072529411​=1.6467602

这个系数是相对于连续旋转13次（按同一方向转）得到的要相乘的系数。

很显然，每次旋转的方向都影响到最终要旋转的累积角度。在–99.7° ≤ θ ≤ 99.7°的范围内的任意角度都可以旋转。

满足法则的所有角度的总和为 99.7° ，对于该范围之外的角度， 可使用三角恒等式转化成该范围内的角度。 当然，角分辨率的数据位数与最终的精度有关。

前面所示的伪旋转现在可以表示为（对每次迭代）：

x i + 1 = x i − d i ( 2 − i y i ) x^{i+1}=x^{i}-d_{i}(2^{-i}y^{i}) xi+1=xi−di​(2−iyi) y i + 1 = y i + d i ( 2 − i x i ) y^{i+1}=y^{i}+d_{i}(2^{-i}x^{i}) yi+1=yi+di​(2−ixi)

在这里我们引入第三个方程，被称为角度累加器，用来在每次迭代过程中追踪累加的旋转角度 :

z i + 1 = z i − d i θ i z^{i+1} = z^{i}-d_{i}\theta^{i} zi+1=zi−di​θi

符号 di 是一个判决算子，用于确定旋转的方向。

因此，原始的算法现在已经被简化为使用向量的伪旋转来表示的迭代移位-相加算法：

x i + 1 = x i − d i ( 2 − i y i ) x^{i+1}=x^{i}-d_{i}(2^{-i}y^{i}) xi+1=xi−di​(2−iyi)

y i + 1 = y i + d i ( 2 − i x i ) y^{i+1}=y^{i}+d_{i}(2^{-i}x^{i}) yi+1=yi+di​(2−ixi)

z i + 1 = z i − d i θ i z^{i+1} = z^{i}-d_{i}\theta^{i} zi+1=zi−di​θi

因此每个迭代需要：

前面提到的去除cosθ的原因是显而易见的。当该项去除时，转换公式已经被简化成伪旋转的迭代移位相加计算。

CORDIC硬件：

伸缩因子是伪旋转的副产物，当简化算法允许伪旋转时，cosθ项被忽略。

由三角函数变换公式：

1 c o s 2 a = 1 + t a n 2 a \frac {1} {cos^{2}a}=1 + tan^{2}a cos2a1​=1+tan2a

可以知道，这样输出xn,yn被伸缩K_n倍，其中： K n = ∏ n 1 ( cos ⁡ θ ( i ) ) = ∏ n ( 1 + 2 ( − 2 i ) ) . K_{n}=\prod_{n} \frac{1} {\left(\cos \theta^{(i)}\right)}=\prod_{n}\left(\sqrt{1+2^{(-2 i)}} \right). Kn​=n∏​(cosθ(i))1​=n∏​(1+2(−2i) ​).

然后如果迭代次数可知，则我们可以预先计算伸缩因子Kn。

同样，1/Kn也可被预先计算得到xn,yn。

旋转的角度和次数不同，Kn的值是不同的，但是是可以确定的，当依次旋转的角度增加，最大达到99.7°时，Kn接近1.6476，i越大2的-i次方越小，旋转的角度越小。

如果我们已知了将被执行的迭代次数，我们便可预先计算出1/Kn的值，并通过1/Kn与xn,yn相乘来校正xn,yn的最终值。

工作模式决定了控制算子di的条件；

CORDIC方法有两种模，分别是旋转模式和向量模式。

旋转模式中选择di=sign(zi)，并最终使zi趋近于0；

n次迭代后我们得到：

X n = K n ( X 0 ∗ c o s ( Z 0 ) − Y 0 ∗ s i n ( Z 0 ) ) X^n=K_n(X^0*cos(Z^0)-Y^0*sin(Z^0)) Xn=Kn​(X0∗cos(Z0)−Y0∗sin(Z0)) Y n = K n ( Y 0 ∗ c o s ( Z 0 ) + X 0 ∗ s i n ( Z 0 ) ) Y^n=K_n(Y^0*cos(Z^0)+X^0*sin(Z^0)) Yn=Kn​(Y0∗cos(Z0)+X0∗sin(Z0)) Z n = 0 Z^n=0 Zn=0

注意：这里Z^0是指最终要旋转的角度之和，其中的角度有正有负，但最终的角度之和是确定的。y0,x0标志没有旋转时候的x,y的坐标值

通过预先设置X0=1/Kn 和 y0=0来计算cos(Z0)和sin(Z0)，这样最终的坐标值就只和初始的没旋转之前的坐标值和最终旋转要达到的角度值有关了。

旋转模式中，判决算子di=sign(zi)，因此，我们输入x(0)和z0（y0=0），然后通过迭代使得z0取值趋近于0。

例如我们这里，当z0=30°时，计算得到的sin z0，cos z0：

当zi 接近于0表示角度的累加是30°了，正如下图所示：

在向量模式中选择：di= -sign(xiyi)，并最终使yi趋近于0；

经过n次迭代后： 通过设定x0=1和z0=0来计算tan-1y0。 向量模式中，判决算子满足下列条件：

d i = − s i g n ( x i y i ) d_{i}=-sign(x^{i}y^{i}) di​=−sign(xiyi)

因此，我们输入x0和y0（z0=0），并通过迭代使y0取值趋近于0；

例如：当y0=0.414210并且x0=1时，计算tan-1（y0/x0）: