曲世隽，王翾

（中国传媒大学信息与通信工程学院，北京 100024）

1 引言

坐标旋转数字计算机（Coordinate Rotation Digital Computer,CORDIC）算法是一种用于高效计算基本函数的迭代算法。该算法通过一般的加减和移位操作实现了乘除法的相关计算，使得向量的的旋转和定向的计算不再依赖于三角函数、乘法、开方、反三角、指数等复杂函数。CORDIC 算法的硬件实现只需要移位器和加法器，简单的硬件就能有效的运算复杂函数。

在1959年，Volder开发了一类计算三角函数、双曲函数的算法，其中包括指数运算和对数运算，这就是CORDIC算法的雏形。在提出这个算法之后，Volder将其应用在导航领域之中。但是在当时这个算法并没有被广泛的讨论和应用。直到1980 年，Haviland 和Tuszynski[1]给出了一种全模式的CORDIC算法的处理芯片。在这之后，CORDIC算法才引起人们的关注，被应用于各种领域。Hu Y H[2]将CORDIC算法应用于基于VLSI的数字信号处理领域。应用包括离散傅里叶变换的计算、矩阵运算和三角函数发生器。在2009 年，Vachhani L,Sridharan K,Meher P K[3]通过FPGA实现高效CORDIC 算法并将其应用在机器人探测上。当然，CORDIC算法也可以应用于其他领域。例如，计算机制图中求点到线的距离、直角坐标与极坐标的相互变换及求多维向量的欧几里得范数等，可以预见，CORDIC算法将在未来有更广泛的应用。

目前人们对于CORDIC算法的关注点主要集中在硬件资源消耗和计算精度两个问题上。想要获得高计算精度就必然要消耗更多的硬件资源，因此如何平衡两者关系就成了一个重要问题。在文献[4]中，提出了一种在硬件资源和计算精度之间寻求折中的方法。文中介绍了两种解决CORDIC算法精度问题的方法，第一种建立在定点数CORDIC算法的运算基础上，文中认为输入变量的不规范会使误差变大，因此需要对输入的变量进行统一、规范的处理，以此提高CORDIC算法的精度，第二种方法是针对于浮点数的CORDIC算法，通过使用混合体系结构降低复杂度，以此提高精度。在实际应用中，为了硬件结构的简单和运算的高效，我们更加偏向于使用定点数进行CORDIC算法。但是文中介绍的提高定点数精度的CORDIC算法在硬件结构上还是过于复杂，因此需要提出一种硬件结构简单，更容易实现的方法来降低误差。

为了提出一种更好的方法实现硬件资源和计算精度的平衡，就需要对CORDIC 算法的误差来源进行详细的分析。CORDIC 算法的误差可以分为近似误差和截断误差，在任何坐标系的任何模式中，近似误差是由有限的迭代次数产生的，而截断误差是由有限的数据位位宽产生的[5]。Tze-Yun Sung,Yi-Hsun Sung[6]对所有坐标系下CORDIC算法的误差都进行了分析。包括圆坐标系，双曲坐标系，线性坐标系的旋转模式和矢量模式，用近似的方式给出了一个误差的最大值。但是文献主要是基于数学分析的角度，给出了近似误差和截断误差的分析结果。在实际CORDIC 算法的应用中，需要有硬件结构的设计和反正切、正余弦函数的计算结果，这些在文章中没有提到。

本文提出了一种降低CORDIC算法中截断误差的方法，将CORDIC算法的流程分为预处理、迭代移位、后处理三个部分，对每个部分进行不同的处理，旨在使用较少的硬件资源同时保持较高的计算精度。文中给出了基于此种方法的硬件架构设计，并做了基于FPGA的仿真实验，验证了结果的可行性和正确性。

2 CORDIC算法原理

CORDIC算法的基本公式如下所示[5]：

其中s(m,i)是非递减的移位序列，满足

是第i次旋转的角度。m= +1,- 1,0确定坐标系是圆，双曲或是线性。δ(i)确定旋转角的方向，+1和-1分别代表不同的方向。由于CORDIC 算法的实现方式很多，以圆坐标系为例，分旋转模式和矢量模式进行原理分析。

通常在旋转模式下，设定x(0) = 1/k1,其中km=，n 为迭代次数。y(0) = 0。z(0)为输入待计算的角度θ。旋转方向δ(i)由z(i)的正负决定，若z(i)为正，则δ(i)为+1，若z(i)为负，则δ(i)为-1。根据公式(1)(2)(3)(5)可以得到角度θ的正弦值y(n)和余弦值x(n)。

通常在矢量模式下，设定输入的值x 和y，通过迭代和旋转，使y(i)逐渐逼近于0。旋转方向δ(i)由y(i)的正负决定，若y(i)为正，则δ(i)为-1，若y(i)为负，则δ(i)为+1。同样根据公式(1)(2)(3)(5)，可以得到反正切函数z(n) = arctan(y/x)，得到向量(x,y)的长度d=x(n)*(1/k1)。

3 CORDIC算法的误差

在研究CORDIC算法的时候，我们需要在硬件资源的使用和计算精度上做出平衡，这就需要对CORDIC算法的误差进行分析。CORDIC算法的误差由迭代次数n，旋转角度小数位位数b，x和y的小数位位数bv确定。在本文中b和bv采用相同的大小。CORDIC算法的误差可以分为两部分：近似误差和截断误差。EA表示近似误差最大值，由有限的迭代次数造成，ER表示截断误差最大值，由有限的数据位位宽造成。

根据文献[6]，圆坐标系下旋转模式的近似误差为

截断误差为

矢量模式下的近似误差为

截断误差为

基于以上的理论推导，可以确定迭代次数和小数位位数对于CORDIC算法的误差的影响，如图1和图2所示。

图1 圆坐标系旋转模式误差

图2 圆坐标系矢量模式误差

4 降低截断误差的方法

4.1 方法简述

从前文的分析中可以得知，增加数据位位宽和迭代次数可以显著的提高CORDIC 算法的精度，但是同时也会消耗更多的硬件资源。降低近似误差，只需要增加迭代的次数，但是降低截断误差却不能只是简单的增加数据位位宽。过于长的数据位位宽不仅会消耗更多的硬件资源，同时也会对后续的数据处理造成影响。因此，需要有一个合适的方法来降低截断误差。

通过此种方式实现CORDIC 算法，完整的步骤可以分为预处理、迭代移位、后处理三个部分。三个步骤如图3所示。

图3 实现改进CORDIC算法完整步骤

为了有更高的精度，我们选择在预处理部分进行数据的规范化处理和扩充数据位位宽。在迭代移位的过程中如果改变数据位位宽，很可能会产生溢出问题，因此选择在后处理部分对迭代移位后的数据进行截断。这样就可以保证在迭代移位部分提高了计算精度，同时输出位宽又可以根据实际的需要进行截取。

4.2 预处理

预处理部分可以分为对输入数据进行规范化处理和对输入数据的位宽进行扩充两个模块。假设输入数据x_in和y_in为16bit 定点数，最高位为符号位。对输入数据进行规范化处理要计算除符号位之外的先导零的个数。比较x_in和y_in的先导零个数，j为两者先导零个数的较小值，根据公式

可以得到j的值。其中m为输入数据位宽，ε= 2-b为精度。在硬件设计中，可以通过一个15-4优先编码器和一个逐位比较的或门来得到j的值。x_in和y_in并行输入逐位比较的或门，可以得到两者之间较高位的1所在的位置。再通过15-4优先编码器得到较高位的1之前先导零的个数。在得到j之后，统一将x_in和y_in左移j位。为了防止迭代过程中发生溢出，还要在符号位后增加两个0 位。根据计算，可以得到k1的值在1.414 和1.646 之间，因此扩充两个0 位完全可以避免溢出。15-4优先编码器的真值表如表1所示。

表1 15-4优先编码器真值表

对输入数据的位宽进行扩充就是在输入数据末尾补零。我们将输入数据扩充为32bit。也就是在进行完数据的规范化处理之后需要在末尾补14 位的0。预处理硬件设计如下图4所示。

图4 预处理硬件设计

4.3 迭代移位

CORDIC 算法的单层迭代移位结构如图5 所示。它需要1 个多路复用器，3 个加/减法器，2 个移位器。旋转方向由selx,sely,selz定义，多路复用器根据旋转模式或矢量模式来进行选择。当模式为旋转模式时，旋转方向由sign(zi)确定。当模式为矢量模式时，旋转方向由sign(yi)确定。要实现并行流水线的CORDIC 算法，只需要对单层的结构进行叠加就可以实现。叠加的次数取决于迭代次数N。

图5 单层迭代结构

4.4 后处理

后处理的部分实际上就是对迭代移位后输出的数据进行位宽的缩减，本质上也是对前导零的去除。缩减的原则是，首先从符号位后最高位开始判断，如果最高位为0，则向左移1 位，继续判断下一位，若还为0，则继续向左移1位，判断下一位，直到遇到为1的位数，则停止判断。对0的位数舍去完之后，整体保留16bit。设后处理部分总共左移的位数为num。后处理部分流程图如图6所示。

图6 后处理流程图

通过流程图可知，实现后处理部分需要通过循环结构来进行判断，为了节省硬件资源，可以使用计数器来实现循环结构。因此，后处理部分需要1 个累加器和1 个计数器，1 个寄存器。将需要处理的数据放入寄存器，在每个计数器+1时进行逐位判断。后处理结构如图7所示。

图7 后处理结构

最后可以总观整个流程中对数据的左移，可以得到数据在预处理部分左移的位数为14 +j,在后处理部分左移的位数为num- 16。所以整个流程对数据的放大为2j+num-2倍。也就是说得到的数据精度为ε=2j+num-2,在改进后的CORDIC 算法中，可以在预处理和后处理部分去除无意义的前导零，以此提升精度。

5 仿真与测试

在本设计中采用Spartan 3E 开发板，用Verilog 对CORDIC 算法在圆坐标系下的矢量模式进行了描述，采用流水线结构。对于100000 组输入范围在(0,10000)的独立均匀分布的16bit 数据进行基于FPGA的圆坐标系矢量模式下的CORDIC 算法计算。选择迭代次数为16 次，在预处理时将其扩充为32bit，在后处理时保留16bit。为了验证其正确性，对其实现过程用ModelSim10.5进行仿真，仿真结果如图8所示。

图8 ModelSim仿真结果

从仿真结果中可以看出输入和输出结果之间存在延时，这跟预处理和迭代移位过程有关。

提取仿真过程中激励文件随机生成的100000 组输入值，对比理想输出相位值、幅度值和实际输出相位值、幅度值之间的绝对误差。其绝对误差的分布如图9所示。

图9 使用截断方法的CORDIC算法的绝对误差分布

计算得到幅度值的平均绝对误差为1.67*10-2，方差为1.3956*10-4。相位值的平均绝对误差为1.1*10-3，方差为2.0246*10-6。对比使用相同数据位位宽的传统CORDIC算法，在精度上有了一定的提升。

6 结束语

该文简要介绍了CORDIC 算法的原理，从硬件实际应用角度出发，提出了一个降低CORDIC 算法截断误差的方法。实现了对复杂函数的高精度运算，并给出了圆坐标系矢量模式下的仿真结果。在实验过程中，可以通过改变迭代的次数或是数据位位宽来满足不同的精度要求。从仿真结果来看，提出的方法可以很好的在精度要求和硬件资源消耗上达到平衡。这个方法在本文中只给出了在圆坐标系矢量模式的应用分析，在以后的研究中，可以将这个方法推广到圆坐标系、线性坐标系、双曲坐标系的不同模式中。为CORDIC 算法计算开方、反三角等复杂函数问题提供了一个新的处理方法，具有一定的借鉴意义。