CORDIC算法:一种高效计算三角函数值的方法 |
您所在的位置:网站首页 › 三角函数cos2x公式 › CORDIC算法:一种高效计算三角函数值的方法 |
By Long Luo 任何一款科学计算器上都可以计算三角函数,三角函数应用于生活工作中的方方面面,但计算机是如何计算三角函数值的呢? 三角函数本质上是直角三角形的3条边的比例关系,计算并没有很直观,那么计算机是如何计算三角函数值的呢? 在微积分中我们学习过 泰勒公式 ,我们知道可以使用泰勒展开式来对某个值求得其近似值,例如: sin∠18∘=5−14≈0.3090169943\sin \angle 18^{\circ} = \frac{\sqrt {5} - 1}{4} \approx 0.3090169943 sin∠18∘=45−1≈0.3090169943 利用泰勒公式,取前 444 项: sinx≈x−x33!+x55!−x77!\sin x \approx x - \frac{x^3}{3!} + \frac{x^5}{5!} - \frac{x^7}{7!} sinx≈x−3!x3+5!x5−7!x7 代入可得: sin∠18∘=sinπ10=π10−16(π10)3+1120(π10)5−15040(π10)7≈0.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 sin∠18∘=sin10π=10π−61(10π)3+1201(10π)5−50401(10π)7≈0.30903399538 可见取前 444 项时精度已经在 10−510^{-5}10−5 之下,对于很多场合精度已经足够高了。 在没有了解 CORDIC(Coordinate Rotation Digital Computer) Algorithm [1] 算法之前,我一直以为计算器是利用泰勒公式去求解,最近才了解到原来在计算机中,CORDIC 算法远比泰勒公式高效。 下面我们就来了解下泰勒公式的不足之处和 CORDIC 算法是怎么做的。 泰勒公式的缺点上一节我们使用泰勒公式去计算三角函数值时,需要做很多次乘法运算,而计算器中乘法运算是很昂贵的,其缺点如下所示: 展开过小则会导致展开精度不够,展开过大则会导致计算量过大; 幂运算需要使用乘法器,存在很多重复计算; 需要很多变量来存储中间值。在之前的文章 矩阵乘法的 Strassen 算法 ,就是通过减少乘法,多做加法,从而大大降低了运算量,那么我们可以用相同的思想来优化运算吗? 当然可以,让我们请出今天的主角:CORDIC 算法。 解析 CORDIC 算法我们知道单位圆上一点 PPP ,其坐标为:(cosθ,sinθ)(\cos \theta , \sin \theta )(cosθ,sinθ) ,如下图所示: 如果我们接收一个旋转向量 MinM_{in}Min 逆时针旋转 θ\thetaθ ,将点 P(xin,yin)P(x_{in} , y_{in})P(xin,yin) 旋转到 P′(xR,yR)P'(x_{R} , y_{R})P′(xR,yR) , 如下图所示: 很容易得到如下公式: xR=xincos(θ)−yinsin(θ)x_R = x_{in} cos(\theta) - y_{in} sin(\theta) xR=xincos(θ)−yinsin(θ) yR=xinsin(θ)+yincos(θ)y_R = x_{in} sin(\theta) + y_{in} cos(\theta) yR=xinsin(θ)+yincos(θ) 实际上由 复数运算 ,我们知道复数乘法就是幅角相加,模长相乘。我们可以将上式写成下列矩阵运算形式: [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} [xRyR]=[cos(θ)sin(θ)−sin(θ)cos(θ)][xinyin] 但上式运算时,只是对向量 vin=[xinyin]v_{in} = {\begin{bmatrix} x_{in} \\ y_{in} \end{bmatrix}}vin=[xinyin] 进行了线性变换,乘以一个旋转向量 MinM_{in}Min ,得到了旋转后的结果:向量 vR=[xRyR]v_{R} = {\begin{bmatrix} x_{R} \\ y_{R} \end{bmatrix}}vR=[xRyR] 。 但是上式仍然需要 444 次乘法和 222 次加减法操作, 复杂度没有任何降低,那怎么办呢? 当当…当! 通过上述分析,我们已经知道可以使用有限次旋转操作来避免复杂的乘法操作,我们修改矩阵运算公式,提取 cos(θ)\cos (\theta )cos(θ) ,则公式可以修改为: [xRyR]=cos(θ)[1−tan(θ)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} [xRyR]=cos(θ)[1tan(θ)−tan(θ)1][xinyin] 如果我们选择合适的角度值 θi\theta_iθi,使得 tan(θi)=2−i,i=0,1,…,n\tan (\theta_{i}) = 2^{-i} , i=0, 1,\dots , n tan(θi)=2−i,i=0,1,…,n 这样和 tan(θi)\tan (\theta_{i})tan(θi) 乘法操作就变成了移位操作,我们知道计算机中移位操作是非常快的,就可以大大加快计算速度了。 但这里仍然有 333 个问题需要解决: 对于任意角度 ,可以通过满足条件的角度累加来得到在数学上相同的结果吗? 每次旋转得到的结果仍然需要乘以 cos(θ)\cos(\theta )cos(θ) ,这部分的计算成本如何?如何计算? 因为每次旋转角度 θ=arctan(2−i)\theta = \arctan(2^{-i})θ=arctan(2−i) ,朝着目标角度进行旋转时,可能会出现没有超过目标角度的情况,也会存在超过目标角度的情况,这种情况如何解决呢?对于第一个问题,答案是否定的。可以从数学上证明只有 ∠45∘\angle 45^{\circ}∠45∘ 的倍数角才可以得到完全一致的结果。但是在工程应用中,我们只需要满足一定精度即可,可以增加迭代次数无限逼近原始角度,如下所示提高 nnn 值以无限逼近原始角度。 θd=∑i=0nθi,∀tan(θi)=2−i\theta_{d} = \sum_{i=0}^{n} \theta_{i} , \forall \tan(\theta_{i}) = 2^{-i} θd=i=0∑nθi,∀tan(θi)=2−i 对于第二个问题,我们先来个例子,以 57.535∘57.535^{\circ}57.535∘ 为例来看看求解过程: 57.535∘=45∘+26.565∘−14.03∘57.535^{\circ} = 45^{\circ}+26.565^{\circ}-14.03^{\circ} 57.535∘=45∘+26.565∘−14.03∘ 那么第一次旋转: [x0y0]=cos(45∘)[1−111][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} [x0y0]=cos(45∘)[11−11][xinyin] 第二次旋转: [x1y1]=cos(26.565∘)[1−2−12−11][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} [x1y1]=cos(26.565∘)[12−1−2−11][x0y0] 第三次旋转: [x2y2]=cos(−14.03∘)[12−2−2−21][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(−14.03∘)[1−2−22−21][x1y1] 综合可得: [x2y2]=cos(45∘)cos(26.565∘)cos(−14.03∘)[1−111][1−2−12−11][12−2−2−21][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} [x2y2]=cos(45∘)cos(26.565∘)cos(−14.03∘)[11−11][12−1−2−11][1−2−22−21][xinyin] 因为 tan(θi)=2−i\tan (\theta_{i}) = 2^{-i}tan(θi)=2−i ,由三角公式可以计算出: cos(θi)=11+tan2(θi)=11+2−2i\cos(\theta_{i}) = \frac {1}{\sqrt {1 + \tan ^{2}(\theta_{i})}} = \frac {1}{\sqrt {1 + 2^{-2i}}} cos(θi)=1+tan2(θi)1=1+2−2i1 令 Ki=cos(θi)K_i = \cos(\theta_{i})Ki=cos(θi) ,则当进行 nnn 次迭代之后: K(n)=∏i=0n−1Ki=∏i=0n−111+2−2i K(n) = \prod _{i=0}^{n-1}K_{i} = \prod _{i=0}^{n-1}{\frac {1}{\sqrt {1 + 2^{-2i}}}} K(n)=i=0∏n−1Ki=i=0∏n−11+2−2i1 当 θi\theta_{i}θi 越来越小时, cosθ\cos \thetacosθ 也越来越逼近 111 ,当迭代次数 n→∞n \to \inftyn→∞ , K(n)K(n)K(n) 极限存在,求解可得: K=limn→∞K(n)≈0.6072529350088812561694K = \lim _{n \to \infty }K(n) \approx 0.6072529350088812561694 K=n→∞limK(n)≈0.6072529350088812561694 由 KKK 我们实际上可以得到最终的向量 vRv_RvR 的模长极限为: A=1K=limn→∞∏i=0n−11+2−2i≈1.64676025812107A = {\frac {1}{K}} = \lim _{n \to \infty } \prod _{i=0}^{n - 1}{\sqrt {1 + 2^{-2i}}} \approx 1.64676025812107 A=K1=n→∞limi=0∏n−11+2−2i≈1.64676025812107 实际上当迭代次数为 666 时,可以计算出缩放比例 KKK ,就已经精确到 0.60720.60720.6072 了,如下所示: K≈cos(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 K≈cos(45∘)cos(26.565∘)×⋯×cos(0.895∘)=0.6072 实际上,任意角度只要迭代次数超过 666 ,我们可以直接使用 K=0.6072K = 0.6072K=0.6072 这个值。 对于第三个问题,稍微有点复杂,我们在下一节继续讲解! 角度累加上一节遗留的问题是迭代旋转角度时,旋转角度不一定会落在目标角度内,我们需要引入一个角度误差,用来衡量旋转角度和目标角之间距离,如下所示: θerror=θd−∑i=0nθi\theta_{error} = \theta_d - \sum_{i=0}^{n} \theta_{i} θerror=θd−i=0∑nθi 当 θerror>0\theta_{error} > 0θerror>0 时,我们应该逆时针旋转,而 θerror |
今日新闻 |
推荐新闻 |
CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3 |