基于FPGA的Cordic算法实现的设计与验证

CORDIC（Coordina     te   Rota         TI   on Digital Computer）算法即坐标旋转数字计算方法，是J.D.Volder1于1959年首次提出，主要用于三角函数、双曲线、指数、对数的计算。该算法通过基本的加和移位运算代替乘法运算，使得矢量的旋转和定向的计算不再需要三角函数、乘法、开方、反三角、指数等函数。

本文是基于FPGA实现Cordic算法的设计与验证，使用Verilog HDL设计，初步可实现正弦、余弦、反正切函数的实现。将复杂的运算转化成FPGA擅长的加减法和乘法，而乘法运算可以用移位运算代替。Cordic算法有两种模式，旋转模式和向量模式。可以在圆坐标系、线性坐标系、双曲线坐标系使用。本文线初步实现在圆坐标系下的两种模式的算法实现。

 Cordic算法简化

旋转模式，迭代位移算法。假设有一点P0（x0，y0），经过逆时针旋转角度θ，到达点Pm（xm，ym），我们根据数学运算可以得到公式如下：

xm = x0cosθ - y0sinθ = cosθ（x0 – y0tanθ）

ym = y0cosθ + x0sinθ = cosθ（y0 – x0tanθ）

如果不考虑旋转后的向量模值，只考虑旋转角度，即去掉cosθ，得到如下方程式。这里旋转的角度的正确的，但x和y的值增加。cosθ值是小于等于1的，值大于等于1，所以模值应该增大。我们不能通过适当的数学计算去掉cosθ，但是去掉cosθ项可以方便我们后面的坐标平面旋转的计算。这里称为伪旋转。

xm = x0 – y0tanθ

ym = y0 – x0tanθ

Cordic的方法核心就是伪旋转，将旋转角θ细化成若干个大小固定的角度θi，规定θi满足tanθi = 2^-i，通过一系列的迭代旋转，每次旋转θi，i为迭代次数，规定∑θi的范围即旋转角度θ的范围为［-99.7， 99.7］。如果θ的大于这个范围则可通过三角运算操作转化到该范围的角度。

我们通过事先将所有每次旋转的角度计算出来，由于每次旋转的角度是固定的，所以经过i次旋转的∑θi可能会超过θ，所以就必须设置一个方向值di，如果旋转角度之和已经小于θ，则di为1，下次旋转继续为顺时针旋转，如果旋转角度之和大于θ，则di为-1，下次旋转为逆时针。设置zi+1为旋转剩余角度，zi+1 = z0 – dizi，z0 = θ，随着i值得增大，zi+1会趋向于0时，即旋转结束。di与zi的符号位相同。

采用伪旋转的方法，每次提出一个cosθi，旋转结束后会产生一个∏cosθi的累乘，一旦我们确定了迭代次数，∏cosθi就是一个常数，迭代公式可写为。这是将cosθi提出、tanθi 替换成 2^-i后的结果。di与zi的符号位相同。

xi+1 = xi - di * yi * 2^-i

yi+1 = yi + di * xi * 2^-i

zi+1 = z0 - di * θi

设迭代i = n - 1，那么旋转n次后得到Pm的坐标应该为（xn * ∏cosθi， yn * ∏cosθi）。应为每次迭代都会提出一个cosθi，旋转n次后的xn和yn就会少乘一个∏cosθi，所以实际上最终的Pm坐标角度近似于（xn * ∏cosθi， yn * ∏cosθi）。

xn * ∏cosθi = x0cosθ - y0sinθ

yn * ∏cosθi = y0cosθ + x0sinθ

xn = 1/∏cosθi （x0cosθ – y0sinθ）

yn = 1/∏cosθi （y0cosθ – x0sinθ）

伸缩因子，KN = 1 / ∏cosθi，已知迭代次数，我们可以预先计算KN的值。如下这是博主使用MATLAB计算出的迭代结果数值。

xn =KN （x0cosθ – y0sinθ）

yn = KN （y0cosθ – x0sinθ）

从上表可以得出，我们预先计算出KN的值，然后令x0 = ∏cosθi，y0 = 0，则上述公式可化简为

xn = cosθ

yn = sinθ

即可实现正弦、余弦操作了。

旋转模式

总结一下，Cordic算法旋转模式使用Verilog HDL的实现流程

（1） 确定迭代次数，将每次迭代的角度计算出来，预先定义为参数，为了避免浮点运算，将角度值向左移位16位，取整数部分。

（2） 根据迭代公式进行迭代计算，本设计取16次迭代，从上表可以看出，当迭代次数越大时，1/∏cosθi会趋向于一个确定值。如果对结果精度要求更高，可以设置更高的迭代次数，根据迭代次数，可以将伸缩因子KN = 1/∏cosθi计算出来。同样将其左移16位。

xi+1 = xi - di * yi * 2^-i

yi+1 = yi + di * xi * 2^-i

zi+1 = z0 - di * θi

（3） 设置x0 = ∏cosθi，y0 = 0，则求出x16 = cosθ，y16 = sinθ。

这里需要注意的是，我们在进行迭代运算的时候，将2^-i变成移位运算，对于正余弦来说是有正负的，所以在一开始定义的时候，就应该定义成有符号数，Verilog中也可以定义有符号数，最高位表示符号位，定义如下

迭代     寄存器   定义为有符号数，那么我们移位运算就不能用》》逻辑右移《《逻辑左移或来移位了，而是用》》》算术右移和《《《算术左移。逻辑左移也就相当于算数左移，右边统一添0 ，逻辑右移，左边统一添0 ，算数右移，左边添加的数和符号有关。

例如1010_1010， ［］是添加的位

逻辑左移一位：0101_010［0］

算数左移一位：0101_010［0］

逻辑右移一位：［0］101_0101

算数右移一位：［1］101_0101

迭代运算采用16级流水线，进行运算，最终需要判断输出的正余弦值在哪个象限，前面讲旋转角度θ的范围为［-99.7，99.7］，不在这个范围我们要进行三角运算使其满足这个范围，当输入的角度小于90度即可进行计算，当输入角度大于90度小于180度，将输入角度减去90度并设定当前角度处于第二象限，然后进行计算，当输入角度大于180度小于270度，将输入的角度减去180度设置当前角度处于第三象限，进行计算，当输入的角度大于270度，减去270设置当前角度处于第四象限，进行计算。象限的设定通过quarant寄存器实现。

如果角度在第一象限，sin（x） = sin（a），cos（x） = sin（a）最后的结果x16 = cosθ， y16 = sinθ，这里我想起了那句口诀，一全正，二正弦，三正切，四余弦

如果角度在第二象限，sin（x） = sin（a+90） = cos（a），cos（x） = cos（a+90） = -sin（a）

如果角度在第三象限，sin（x） = sin（a+180） = -sin（a），cos（x） = cos（a+180） = -cos（a）

如果角度在第四象限，sin（x） = sin（a+270） = cos（a），cos（x） = cos（a+270） = -sin（a）

对于正数，我们直接赋值输出，负数，这里使用有符号数表示，将其取反加1即可。最终使用modelsim对算法进行仿真，从波形图上看已经初步实现了sin，cos函数。

向量模式

Cordic算法在向量模式下的计算方法和旋转模式基本上是类似的，设有一点P0（x0， y0），经过顺时针旋转角度到与轴重合，得到点Pm（xm， ym），即ym = 0。

xm = x0cosθ - y0sinθ = cosθ（x0 – y0tanθ）

ym = y0cosθ + x0sinθ = cosθ（y0 – x0tanθ） = 0

我们设置x0 = x， y0 = y， z0 = 0，迭代次数为16，经过16次迭代后得到zn = θ = arctan（y/x）和坐标所代表的向量的模值d = xm = xn * ∏cosθi，di与yi方向相反，即当时结束运算。实现方法为判断yi的符号位，符号位为1，di为1，符号位为0，di为-1。

xi+1 = xi - di * yi * 2^-i

yi+1 = yi + di * xi * 2^-i

zi+1 = z0 - di * θi

关于反正切函数，由于在［-99.7°，99.7°］范围内，所以我们输入向量P0（x0， y0）时，需要保证其在第一、四象限。

下面是使用MATLAB计算出来的数据和FPGA计算出来的数据进行比较。

从FPGA计算出的结果与MATLAB来比较，和实际结果之间的误差还是挺小的，毕竟是硬件计算出来的数据，向量的误差就比较大了，如果对于精度比较高的计算，我们可以通过提高迭代次数来提高精度。至此基于FPGA的Cordic算法就实现结束了。