硬件三角函数算法CORDIC¶

背景介绍¶

出于工作需要，希望能够在硬件上实现三角函数等超越函数的计算，于是查找资料来找寻实现的一般方法，搜索给出的结果出奇的一致，大家普遍都是使用cordic算法实现的。尽管具体的实现会有差异，但核心的思想是一致的。尽管最后出于各种原因，没有把这个工作做到最后，但对于cordic算法的理论已经基本掌握，特此记录。

值得说明的是，网络上这方面的文章已经很多了，所以写这篇文章的目的是把它作为一个思路引导，而不是再把网络上的资料抄一遍。

CORDIC算法¶

据了解，CORDIC(Coordinate Rotation Digital Computer)坐标旋转数字计算算法1959年就被提出了，用来解决导航系统中三角函数运算的问题，用这种方法能够满足实时性要求。那么在2021年我们使用这个算法，显然也能够满足实时性要求。

首先，要了解到其基本原理，用文字描述就是已知一个角度的正弦值a和余弦值b，然后这个角度再旋转一个角度c,旋转之后的角度的正弦值和余弦值能够通过a,b,c这三个量求出来。这提供了逐次逼近的基础。在逐次逼近时，比如我们求sin(30°)，实质上是用一系列已知正余弦值的角度组合来逐次逼近到30°。(当然，这些角度都是一些特定角度)。

举例来讲，30.0° ≈ 45.0° - 26.6° + 14.0° - 7.1° + 3.6° + 1.8° - 0.9° + 0.4° - 0.2° + 0.1° = 30.1°

那么根据上面的原理，知道45.0度的正余弦值，它旋转了26.6°，就求出了旋转之后的(45.0° - 26.6°)的正余弦值，这算一次迭代，当迭代了10次之后，就计算出了30.1°的正余弦值，也就近似是30度的正余弦值。

那么肯定会问，为什么是这些特定角度？因为这些角度都会使计算后的纵坐标缩减为原来的½,这显然是逐次逼近能够快速接近的基础。

当然你也会疑惑，旋转后的角度真的能够通过a,b,c三个量求出来吗？我从上学就养成了一个很不好的习惯，当我不明白原理，只知道结论的时候，就很难下手做一件事情。但其实很多时候，对于快速搞定一件事情来讲，原理并不那么重要，比如说考试的时候。既然现在不是考试，那么我们还是来看一下这里面的原理，先不要怕，初三水平足够了。

原理参见链接，其中它的图的标注存在问题，sin和cos写反了，但除此之外已经非常简单明了https://www.cnblogs.com/touchblue/p/3535968.html

代码实现¶

github中已经有很多版本的实现，这里附带几个链接，方便参考，代码普遍比较简单，一看就懂

cordic_wrapped CORDIC-all-in-one-verilog

关于预处理¶

通过熟悉原理和代码，想必我们还发现了一个细节，就是cordic对于输入的角度是有要求的，一般的应用场景，我们不会对角度做限制，但cordic只能处理[-99.7°,99.7°]之间的角度，因此对于更大的角度，就需要对2Π取余。

这里面涉及除法和Π的精度问题，比如下面的处理

int shift = (int) floor(a/(M_PI/2.0));
double b = a - shift*(M_PI/2.0);

其中a为输入角度，b为取余之后的角度，假如使用fp16来表示M_PI，a数值较大时，必然会带来精度问题,最简单的做法是内部计算采用fp32,计算完成后再用fp16进行表示。

关于象限¶

对于计算不同象限的角度，在符号上是不同的，详细的处理可以参考给出的github代码。