CORDIC算法:一种高效计算三角函数值的方法

By Long Luo

任何一款科学计算器上都可以计算三角函数,三角函数应用于生活工作中的方方面面,但计算机是如何计算三角函数值的呢?

三角函数本质上是直角三角形的3条边的比例关系,计算并没有很直观,那么计算机是如何计算三角函数值的呢?

在微积分中我们学习过 泰勒公式 ,我们知道可以使用泰勒展开式来对某个值求得其近似值,例如:

sin18=5140.3090169943\sin \angle 18^{\circ} = \frac{\sqrt {5} - 1}{4} \approx 0.3090169943

利用泰勒公式,取前 44 项:

sinxxx33!+x55!x77!\sin x \approx x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!}

代入可得:

sin18=sinπ10=π1016(π10)3+1120(π10)515040(π10)70.30903399538\sin \angle 18^{\circ} = \sin \frac{\pi}{10} = \frac{\pi}{10} - \frac{1}{6} (\frac{\pi}{10})^3 + \frac{1}{120} (\frac{\pi}{10})^5 - \frac{1}{5040} (\frac{\pi}{10})^7 \approx 0.30903399538

可见取前 44 项时精度已经在 10510^{-5} 之下,对于很多场合精度已经足够高了。

在没有了解 CORDIC(Coordinate Rotation Digital Computer) Algorithm [1] 算法之前,我一直以为计算器是利用泰勒公式去求解,最近才了解到原来在计算机中,CORDIC 算法远比泰勒公式高效。

下面我们就来了解下泰勒公式的不足之处和 CORDIC 算法是怎么做的。

泰勒公式的缺点

上一节我们使用泰勒公式去计算三角函数值时,需要做很多次乘法运算,而计算器中乘法运算是很昂贵的,其缺点如下所示:

  1. 展开过小则会导致展开精度不够,展开过大则会导致计算量过大;
  2. 幂运算需要使用乘法器,存在很多重复计算;
  3. 需要很多变量来存储中间值。

在之前的文章 矩阵乘法的 Strassen 算法 ,就是通过减少乘法,多做加法,从而大大降低了运算量,那么我们可以用相同的思想来优化运算吗?

当然可以,让我们请出今天的主角:CORDIC 算法

解析 CORDIC 算法

我们知道单位圆上一点 PP ,其坐标为:(cosθ,sinθ)(\cos \theta , \sin \theta ) ,如下图所示:

Unit Cycle

如果我们接收一个旋转向量 MinM_{in} 逆时针旋转 θ\theta ,将点 P(xin,yin)P(x_{in} , y_{in}) 旋转到 P(xR,yR)P'(x_{R} , y_{R}) , 如下图所示:

CORDIC

很容易得到如下公式:

xR=xincos(θ)yinsin(θ)x_R = x_{in} cos(\theta) - y_{in} sin(\theta)

yR=xinsin(θ)+yincos(θ)y_R = x_{in} sin(\theta) + y_{in} cos(\theta)

实际上由 复数运算 ,我们知道复数乘法就是幅角相加,模长相乘。我们可以将上式写成下列矩阵运算形式:

[xRyR]=[cos(θ)sin(θ)sin(θ)cos(θ)][xinyin]\begin{aligned} \begin{bmatrix} x_{R} \\ y_{R} \end{bmatrix} = \begin{bmatrix} \cos (\theta) & -\sin (\theta) \\ \sin (\theta) & \cos (\theta ) \end{bmatrix} \begin{bmatrix} x_{in} \\ y_{in} \end{bmatrix} \end{aligned}

但上式运算时,只是对向量 vin=[xinyin]v_{in} = {\begin{bmatrix} x_{in} \\ y_{in} \end{bmatrix}} 进行了线性变换,乘以一个旋转向量 MinM_{in} ,得到了旋转后的结果:向量 vR=[xRyR]v_{R} = {\begin{bmatrix} x_{R} \\ y_{R} \end{bmatrix}}

但是上式仍然需要 44 次乘法和 22 次加减法操作, 复杂度没有任何降低,那怎么办呢?

当当…当!

通过上述分析,我们已经知道可以使用有限次旋转操作来避免复杂的乘法操作,我们修改矩阵运算公式,提取 cos(θ)\cos (\theta ) ,则公式可以修改为:

[xRyR]=cos(θ)[1tan(θ)tan(θ)1][xinyin]\begin{bmatrix} x_{R} \\ y_{R} \end{bmatrix} = cos(\theta) \begin{bmatrix} 1 & -tan(\theta) \\ tan(\theta) & 1 \end{bmatrix} \begin{bmatrix} x_{in} \\ y_{in} \end{bmatrix}

如果我们选择合适的角度值 θi\theta_i,使得

tan(θi)=2ii=0,1,,n\tan (\theta_{i}) = 2^{-i} \quad i=0, 1,\dots , n

这样和 tan(θi)\tan (\theta_{i}) 乘法操作就变成了移位操作,我们知道计算机中移位操作是非常快的,就可以大大加快计算速度了。

但这里仍然有 33 个问题需要解决:

  1. 对于任意角度 ,可以通过满足条件的角度累加来得到在数学上相同的结果吗?
  2. 每次旋转得到的结果仍然需要乘以 cos(θ)\cos(\theta ) ,这部分的计算成本如何?如何计算?
  3. 因为每次旋转角度 θ=arctan(2i)\theta = \arctan(2^{-i}) ,朝着目标角度进行旋转时,可能会出现没有超过目标角度的情况,也会存在超过目标角度的情况,这种情况如何解决呢?

CORDIC Expand

对于第一个问题,答案是否定的。可以从数学上证明只有 45\angle 45^{\circ} 的倍数角才可以得到完全一致的结果。但是在工程应用中,我们只需要满足一定精度即可,可以增加迭代次数无限逼近原始角度,如下所示提高 nn 值以无限逼近原始角度。

θd=i=0nθitan(θi)=2i\theta_{d} = \sum_{i=0}^{n} \theta_{i} \quad \forall \tan(\theta_{i}) = 2^{-i}

对于第二个问题,我们先来个例子,以 57.53557.535^{\circ} 为例来看看求解过程:

57.535=45+26.56514.0357.535^{\circ} = 45^{\circ}+26.565^{\circ}-14.03^{\circ}

那么第一次旋转:

[x0y0]=cos(45)[1111][xinyin]\begin{bmatrix} x_{0} \\ y_{0} \end{bmatrix} = cos(45^{\circ}) \begin{bmatrix} 1 & -1 \\ 1 & 1 \end{bmatrix} \begin{bmatrix} x_{in} \\ y_{in} \end{bmatrix}

第二次旋转:

[x1y1]=cos(26.565)[121211][x0y0]\begin{bmatrix} x_{1} \\ y_{1} \end{bmatrix} = cos(26.565^{\circ}) \begin{bmatrix} 1 & -2^{-1} \\ 2^{-1} & 1 \end{bmatrix} \begin{bmatrix} x_{0} \\ y_{0} \end{bmatrix}

第三次旋转:

[x2y2]=cos(14.03)[122221][x1y1]\begin{bmatrix} x_{2} \\ y_{2} \end{bmatrix} = cos(-14.03^{\circ}) \begin{bmatrix} 1 & 2^{-2} \\ -2^{-2} & 1 \end{bmatrix} \begin{bmatrix} x_{1} \\ y_{1} \end{bmatrix}

综合可得:

[x2y2]=cos(45)cos(26.565)cos(14.03)[1111][121211][122221][xinyin]\begin{bmatrix} x_{2} \\ y_{2} \end{bmatrix} = cos(45^{\circ}) cos(26.565^{\circ}) cos(-14.03^{\circ}) \begin{bmatrix} 1 & -1 \\ 1 & 1 \end{bmatrix} \begin{bmatrix} 1 & -2^{-1} \\ 2^{-1} & 1 \end{bmatrix} \begin{bmatrix} 1 & 2^{-2} \\ -2^{-2} & 1 \end{bmatrix} \begin{bmatrix} x_{in} \\ y_{in} \end{bmatrix}

因为 tan(θi)=2i\tan (\theta_{i}) = 2^{-i} ,由三角公式可以计算出:

cos(θi)=11+tan2(θi)=11+22i\cos(\theta_{i}) = \frac {1}{\sqrt {1 + \tan ^{2}(\theta_{i})}} = \frac {1}{\sqrt {1 + 2^{-2i}}}

Ki=cos(θi)K_i = \cos(\theta_{i}) ,则当进行 nn 次迭代之后:

K(n)=i=0n1Ki=i=0n111+22i K(n) = \prod _{i=0}^{n-1}K_{i} = \prod _{i=0}^{n-1}{\frac {1}{\sqrt {1 + 2^{-2i}}}}

θi\theta_{i} 越来越小时, cosθ\cos \theta 也越来越逼近 11 ,当迭代次数 nn \to \inftyK(n)K(n) 极限存在,求解可得:

K=limnK(n)0.6072529350088812561694K = \lim _{n \to \infty }K(n) \approx 0.6072529350088812561694

KK 我们实际上可以得到最终的向量 vRv_R 的模长极限为:

A=1K=limni=0n11+22i1.64676025812107A = {\frac {1}{K}} = \lim _{n \to \infty } \prod _{i=0}^{n - 1}{\sqrt {1 + 2^{-2i}}} \approx 1.64676025812107

实际上当迭代次数为 66 时,可以计算出缩放比例 KK ,就已经精确到 0.60720.6072 了,如下所示:

Kcos(45)cos(26.565)××cos(0.895)=0.6072K \approx cos(45^{\circ}) cos(26.565^{\circ}) \times \dots \times cos(0.895^{\circ}) = 0.6072

实际上,任意角度只要迭代次数超过 66 ,我们可以直接使用 K=0.6072K = 0.6072 这个值。

对于第三个问题,稍微有点复杂,我们在下一节继续讲解!

角度累加

上一节遗留的问题是迭代旋转角度时,旋转角度不一定会落在目标角度内,我们需要引入一个角度误差,用来衡量旋转角度和目标角之间距离,如下所示:

θerror=θdi=0nθi\theta_{error} = \theta_d - \sum_{i=0}^{n} \theta_{i}

θerror>0\theta_{error} > 0 时,我们应该逆时针旋转,而 θerror<0\theta_{error} < 0 ,则顺时针旋转。根据精度需要,当 θerrorϵ\left | \theta_{error} \right | \le \epsilon 即可退出迭代。

同时我们修改之前的公式,引入 σi{+1,1}\sigma_{i} \in \left \{ +1, -1 \right \} ,于是可以得到最终公式:

x[i+1]=x[i]σi2iy[i]x \left [ i+1 \right ] = x \left [ i \right ] - \sigma_{i} 2^{-i} y \left [ i \right ]

y[i+1]=y[i]+σi2ix[i]y \left [ i+1 \right ] = y \left [ i \right ] + \sigma_{i} 2^{-i} x \left [ i \right ]

z[i+1]=z[i]σitan1(2i)z \left [ i+1 \right ] = z \left [ i \right ] - \sigma_{i} tan^{-1} ( 2^{-i} )

举个例子

上面讲了这么多,来个实例吧,练习巩固下知识,看看自己是否真的懂了?

计算 sin70\sin 70^{\circ}cos70\cos 70^{\circ} 的值。

xin=1,yin=0x_{in}=1, y_{in} = 0 开始,迭代 66 次结果如下:

ii 次迭代 σi\sigma_{i} x[i]x \left[ i \right ] y[i]y \left[ i \right ] z[i]z \left[ i \right ]
- - 11 00 7070^{\circ}
00 11 11 11 2525^{\circ}
11 11 0.50.5 1.51.5 1.5651-1.5651^{\circ}
22 1-1 0.8750.875 1.3751.375 12.471112.4711^{\circ}
33 11 0.70310.7031 1.48441.4844 5.34615.3461^{\circ}
44 11 0.61030.6103 1.52831.5283 1.76981.7698^{\circ}
55 11 0.56250.5625 1.54741.5474 0.0201-0.0201^{\circ}
66 1-1 0.58670.5867 1.53861.5386 0.87510.8751^{\circ}

迭代到第 66 次时,角度误差已经小于 11^{\circ} 了, 通过表格可知:

xR=0.6072×0.5867=0.3562x_{R} = 0.6072 \times 0.5867 = 0.3562

yR=0.6072×1.5386=0.9342y_{R} = 0.6072 \times 1.5386 = 0.9342

通过计算器可知,cos(70)=0.34202\cos(70^{\circ}) = 0.34202sin(70)=0.93969\sin(70^{\circ}) = 0.93969 ,误差已经在 1100\frac {1}{100} 之下了,实际应用中我们会迭代 1616 次,误差会非常小

在线 CORDIC 算法Demo

通过上面分析,我们已经知道了 CORDIC 算法的原理,下面就开始编程吧!

用 JavaScript 写了一个在线互动版本,传送门 →

  1. 可以调整不同迭代次数,和系统库函数 Math.cos\textit{Math.cos}Math.sin\textit{Math.sin} 进行比较:

Cordic Results

  1. 可以查看每次迭代的结果,掌握 Cordic 算法迭代原理:

Cordic Iteration Results

CORDIC 算法的优点

CORDIC 算法相比其他算法的优点,体现在以下几个方面:

  1. 简化运算:CORDIC 算法主要使用位移、加法和减法等简单的运算,避免了复杂的乘法操作,从而提高了计算速度;
  2. 并行计算:CORDIC 算法的迭代操作是相互独立的,可以进行并行计算,利用现代计算机的多核优势,进一步提升计算效率;
  3. 硬件优化:CORDIC 算法适用于硬件实现,可以通过专用硬件电路(如FPGA)进行加速,极大地提高计算速度;
  4. 低存储需求:CORDIC 算法只需存储一小组预先计算好的旋转角度,节省了存储空间;
  5. 迭代控制:通过控制迭代次数,可以平衡计算精度和计算速度,根据需求进行调整。

CORDIC 算法的不足

几个月这篇文章发布之后,陆陆续续得到了不少读者的宝贵意见,有读者反馈说没有写 CORDIC 算法的缺点,也有读者反馈在他在嵌入式处理器STM32上还不如 泰勒展开式(Taylor Series) 快。

对于这 2\large {2} 个问题,我又查阅了一些 CORDIC 算法[2] 的资料,可以回答下这 2\large {2} 个问题。

CORDIC 算法的不足之处在于下面几点:

收敛速度慢

CORDIC 算法的收敛速度较慢,比如当输入值 θ\large {\theta } 很接近 0\large {0}π2\large {\frac {\pi}{2}} 时,这个时候就需要进行多轮迭代才能达到足够的精度。

θ0\large {\theta \approx 0 } 时,实际上只需要取前 2\large {2} 项就有足够精度了,而且收敛速度快多了!

sin(x)=x13!x3\sin(x) = x - \frac {1}{3!}x^3

cos(x)=112!x2\cos(x) = 1 - \frac {1}{2!}x^2

精度

如上一小节,当输入值 θ\large {\theta } 是小角度或者大角度的情况时,CORDIC 算法的固定迭代次数可能也无法提供足够的精度。

输入参数范围

CORDIC 算法在实际应用中需要考虑输入参数的限制 ,特别是幅值和角度范围,这是因为 CORDIC 算法是一种迭代算法,对于较大的角度,需要更多的迭代次数才能达到所需的精度,从而导致计算效率的下降。

CORDIC 算法适用于处理相对较小范围的角度,确切来说是在 [99.88,99.88][−99.88^{\circ} , 99.88^{\circ}] 的范围内表现最佳,下面我们来证明下:

由 CORDIC 算法计算过程可知:

z[i]j=iarctan2j\left | z[i] \right | \leq \sum_{j=i}^{\infty } \arctan 2^{-j}

那么可得:

θmax=z[0]max=j=0arctan2j1.7429(99.88)σj=1,z[j]>0\theta_{max} = z[0]_{max} = \sum_{j=0}^{\infty } \arctan 2^{-j} \approx 1.7429(99.88^{\circ}) \quad \quad \sigma_j = 1, z[j] > 0

假设 θ<θmax\large {\theta < \theta_{max}} , 则有:

z[i]arctan(2(i1))\left | z[i] \right | \leq \arctan (2^{−(i−1)})

因此

arctan(2i)j=i+1arctan(2j)i\arctan (2^{−i}) \leq \sum_{j=i+1}^{\infty } \arctan (2^{−j}) \quad \forall i

所以 CORDIC算法的输入范围是:[99.88,99.88][−99.88^{\circ} , 99.88^{\circ}]

CORDIC 算法的FPGA实现

对于第 2\large {2} 个问题,我不太清楚那位读者是如何测试的!如之前所述,在某些情况下,泰勒展开式是优于 CORDIC 算法的!

CORDIC 算法是在硬件中实现的,是 FPGA 设计[3] 中的经典例题。

总结

CORDIC 算法是一种高效计算三角函数值的方法。相比传统的泰勒展开式,它具有简单高效、低存储需求和迭代控制等优势。在需要计算三角函数值的应用中,CORDIC 算法更快速、更节省资源。

参考文献


  1. CORDIC ↩︎

  2. Digital Circuits/CORDIC ↩︎

  3. Using a CORDIC to calculate sines and cosines in an FPGA ↩︎