CN103713878B - 一种应用补码方法的正余弦cordic算法在fpga实现的方法 - Google Patents

Abstract

本发明公开了一种应用补码方法的正余弦CORDIC算法在FPGA实现的方法，该方法包括：采用数学归纳法由CORDIC算法推导正余弦CORDIC算法；采用补码法判断正余弦CORDIC算法坐标旋转方向，根据坐标旋转方向求出旋转坐标量，获取旋转坐标运算结果；通过对正余弦CORDIC算法的仿真实验，对旋转坐标运算结果进行分析。本发明能有效解决迭代法不停迭代过程所消耗很多时间资源与流水线所消耗的逻辑资源的缺点，减少比较次数，减少时间延迟，能更好地解决实时性问题，在信号实时计算中具有重要的实际意义。

Description

一种应用补码方法的正余弦CORDIC算法在FPGA实现的方法

技术领域

本发明涉及一种应用补码方法的正余弦CORDIC算法，尤其涉及一种应用补码方法的正余弦CORDIC算法在FPGA实现的方法。

背景技术

坐标旋转数字计算（Coordinate Rotation Digital Computer,CORDIC）算法是1959年由Volder J提出，通过一系列固定与运算基数有关角度不断偏摆迭代以逼近所需旋转角度，由简单加法、移位操作即可完成三角函数、指数函数、对数函数等超越函数的计算。1971年Walther J D把圆周旋转、双曲旋转和直线旋转统一成标准CORDIC迭代方程，进一步增加算法运算功能。CORDIC算法具有规则化、结构简单的特点，容易通过硬件方法实现。流水线法、迭代法是目前CORDIC算法在FPGA实现的主要方法。其中流水线法基于并行结构，处理速度快，而带来结构复杂、占用硬件资源多的问题，初始时延较大；迭代法基于串行结构，结构简单、占用资源少，则存在迭代次数多、运算速度慢的问题，初始时延较大。为提高CORDIC算法在FPGA实现性能，国内外开展很多有益研究，从不同角度对CORDIC算法进行结构改进，减少存储单元、降低功耗，但部分仍存在算法复杂的问题。为此，本发明针对CORDIC算法在FPGA实现存在问题，研究通项式正余弦CORDIC算法，并在FPGA实现本算法。

发明内容

为解决上述技术问题，本发明的目的是提供一种应用补码方法的正余弦CORDIC（Coordinate Rotation Digital Computer，坐标旋转数字计算）算法在FPGA实现的方法，该方法降低结构复杂度、减少迭代次数和时间延迟，提高实时性。

本发明的目的通过以下的技术方案来实现：

一种应用补码方法的正余弦CORDIC算法在FPGA实现的方法，包括：

A采用数学归纳法由CORDIC算法推导正余弦CORDIC算法；

B采用补码法判断正余弦CORDIC算法坐标旋转方向，根据坐标旋转方向求出旋转坐标量，获取旋转坐标运算结果；

C通过对正余弦CORDIC算法的仿真实验，对旋转坐标运算结果进行分析。

与现有技术相比，本发明的一个或多个实施例可以具有如下优点：能有效解决迭代法不停迭代过程所消耗很多时间资源与流水线所消耗的逻辑资源的缺点，减少比较次数，减少时间延迟，能更好的解决实时性问题，在信号实时计算中有重要的实际意义。

附图说明

附图用来提供对本发明的进一步理解，并且构成说明书的一部分，与本发明的实施例共同用于解释本发明，并不构成对本发明的限制。在附图中：

图1是应用补码方法的正余弦CORDIC算法在FPGA实现的方法流程图；

图2是直角坐标系旋转向量图；

图3是基于通项法的正余弦函数CORDIC算法流程图；

图4是传统流水线法、传统迭代法和通项法的仿真模拟波形图；

图5是传统流水线法、传统迭代法和通项法的仿真延时情况图。

具体实施方式

容易理解，根据本发明的技术方案，在不变更本发明的实质精神下，本领域的一般技术人员可以提出本发明的多个结构方式和制作方法。因此以下具体实施方式以及附图仅是本发明的技术方案的具体说明，而不应当视为本发明的全部或者视为本发明技术方案的限定或限制。

下面结合实施例及附图对本发明作进一步详细的描述。

图1是应用补码方法的正余弦CORDIC算法在FPGA实现的方法，该包括：

步骤10将CORDIC算法采用数学归纳法推导正余弦CORDIC算法；

在直角坐标系旋转向量图中（如图2所示），单位向量p₀(x₀,y₀)旋转角θ到p_n(x_n,y_n)，可表示为：

$\left[\begin{array}{c}{x}_n\\ y_n\end{array}\right]=\cos \mathrm{\θ}\left[\begin{array}{cc}1& -\tan \mathrm{\θ}\\ \tan \mathrm{\θ}& 1\end{array}\right]\left[\begin{array}{c}{x}_0\\ y_0\end{array}\right]---\left(1\right)$

将角度θ分成n个递减小角度θ_i=arctan(2^-i)之和，即方向系数d_i顺时针、逆时针分别取-1、1。令k_i＝cosθ_i，引入角度中间变量z_i，可由(1)式，可推得CORDIC公式：

$\left\{\begin{array}{c}{x}_{i+1}=k_i(x_i-y_i{d}_i\tan {\mathrm{\θ}}_i)\\ y_{i+1}=k_i(y_i+x_i{d}_i\tan {\mathrm{\θ}}_i)\\ z_{i+1}=\mathrm{\θ}-\underset{0}{\overset{i}{\mathrm{\Σ}}}{d}_i{\mathrm{\θ}}_i\end{array}\right.$

由于经过n次旋转后变量z_n+1→0，令校模因子则有：

$\left\{\begin{array}{c}{x}_{n+1}=K[x_0\cos \left(z_0\right)-y_0\sin \left(z_0\right)]\\ y_{n+1}=K[y_0\cos \left(z_0\right)+x_0\sin \left(z_0\right)]\\ z_{n+1}=0\end{array}\right.---\left(2\right)$

况且有n→∞，K=0.607253。在初值x₀=1/K，y₀=0，z₀=θ情况下，则通过n次迭代后CORDIC公式输出：

[x_n+1,y_n+1,z_n+1]＝[cosθ,sinθ,0] (3)

（2）在 $\left(\begin{array}{c}{x}_0\\ y_0\end{array}\right)=\left(\begin{array}{c}1\\ 0\end{array}\right)$ 情况下，基于CORDIC算法递推公式，有：

$\begin{array}{c}\left(\begin{array}{c}{x}_{i+1}\\ y_{i+1}\end{array}\right)=k_i\left(\begin{array}{cc}1& -\frac{d_i}{2^i}\\ \frac{d_i}{2^i}& 1\end{array}\right)\left(\begin{array}{c}{x}_i\\ y_i\end{array}\right)\\ \⇒\left(\begin{array}{c}{x}_{i+1}\\ y_{i+1}\end{array}\right)=\underset{i=0}{\overset{i}{\mathrm{\Π}}}{k}_i\left(\begin{array}{cc}1& -\frac{d_i}{2^i}\\ \frac{d_i}{2^i}& 1\end{array}\right)...\left(\begin{array}{cc}1& -\frac{d_0}{2^0}\\ \frac{d_0}{2^0}& 1\end{array}\right)\left(\begin{array}{c}1\\ 0\end{array}\right)\end{array}$

可推得：

$\left(\begin{array}{c}{x}_{n+1}\\ y_{n+1}\end{array}\right)=\underset{i=0}{\overset{n}{\mathrm{\Π}}}{k}_i\left(\begin{array}{c}1-\underset{N=1}{\overset{\mathrm{int}\left(\frac{n+1}{2}\right)}{\mathrm{\Σ}}}\left\{{(-1)}^N\right[\underset{i=1}{\overset{n-2N+2}{\mathrm{\Σ}}}...\underset{j=i+2N-1}{\overset{n}{\mathrm{\Σ}}}(\frac{d_i}{2^i}...\frac{d_j}{2^j})\left]\right\}\\ \underset{N=0}{\overset{\mathrm{int}\left(\frac{n+1}{2}\right)}{\mathrm{\Σ}}}\left\{{(-1)}^N\right[\underset{i=1}{\overset{n+1-2N}{\mathrm{\Σ}}}...\underset{j=i+2N}{\overset{n}{\mathrm{\Σ}}}(\frac{d_i}{2^i}...\frac{d_j}{2^j})\left]\right\}\end{array}\right)---\left(4\right)$

该式可用数学归纳法加以证明。

①i=1，由式(4)有：

$\left(\begin{array}{c}{x}_1\\ y_1\end{array}\right)=k_0\left(\begin{array}{cc}1& -\frac{d_0}{2^0}\\ \frac{d_0}{2^0}& 1\end{array}\right)\left(\begin{array}{c}1\\ 0\end{array}\right)=k_0\left(\begin{array}{c}1\\ \frac{d_0}{2^0}\end{array}\right)$

②假设i=n-1时等式成立，即：

$\left(\begin{array}{c}{x}_n\\ y_n\end{array}\right)=\underset{i=0}{\overset{n-1}{\mathrm{\Π}}}{k}_i\·\left(\begin{array}{c}1-\underset{N=1}{\overset{\mathrm{int}(n/2)}{\mathrm{\Σ}}}\left\{{(-1)}^N\right[\underset{i=0}{\overset{n-2N+1}{\mathrm{\Σ}}}...\underset{j=i+2N-1}{\overset{n-1}{\mathrm{\Σ}}}(\frac{d_i}{2^i}...\frac{d_j}{2^j})\left]\right\}\\ \underset{N=0}{\overset{\mathrm{int}(n/2)}{\mathrm{\Σ}}}\left\{{(-1)}^N\right[\underset{i=0}{\overset{n-2N}{\mathrm{\Σ}}}...\underset{j=i+2N}{\overset{n-1}{\mathrm{\Σ}}}(\frac{d_i}{2^i}...\frac{d_j}{2^j})\left]\right\}\end{array}\right)$

当i=n有:

$\left(\begin{array}{c}{x}_{n+1}\\ y_{n+1}\end{array}\right)=\underset{i=0}{\overset{n-1}{\mathrm{\Π}}}{k}_i\left(\begin{array}{cc}1& -\frac{d_n}{2^n}\\ \frac{d_n}{2^n}& 1\end{array}\right)\left(\begin{array}{c}1-\underset{N=1}{\overset{\mathrm{int}(n/2)}{\mathrm{\Σ}}}\left\{{(-1)}^N\right[\underset{i=0}{\overset{n-2N+1}{\mathrm{\Σ}}}...\underset{j=i+2N-1}{\overset{n-1}{\mathrm{\Σ}}}(\frac{d_i}{2^i}...\frac{d_j}{2^j})\left]\right\}\\ \underset{N=0}{\overset{\mathrm{int}(n/2)}{\mathrm{\Σ}}}\left\{{(-1)}^N\right[\underset{i=0}{\overset{n-2N}{\mathrm{\Σ}}}...\underset{j=i+2N}{\overset{n-1}{\mathrm{\Σ}}}(\frac{d_i}{2^i}...\frac{d_j}{2^j})\left]\right\}\end{array}\right)$

$=\underset{i=0}{\overset{n}{\mathrm{\Π}}}{k}_i\left(\begin{array}{c}1-\underset{N=1}{\overset{\mathrm{int}\left(\frac{n+1}{2}\right)}{\mathrm{\Σ}}}\left\{{(-1)}^N\right[\underset{i=0}{\overset{n-2N+2}{\mathrm{\Σ}}}...\underset{j=i+2N-1}{\overset{n}{\mathrm{\Σ}}}(\frac{d_i}{2^i}...\frac{d_j}{2^j})\left]\right\}\\ \underset{N=0}{\overset{\mathrm{int}\left(\frac{n+1}{2}\right)}{\mathrm{\Σ}}}\left\{{(-1)}^N\right[\underset{i=0}{\overset{n+1-2N}{\mathrm{\Σ}}}...\underset{j=i+2N}{\overset{n}{\mathrm{\Σ}}}(\frac{d_i}{2^i}...\frac{d_j}{2^j})\left]\right\}\end{array}\right)$

故式(4)成立。并且在n足够大情况下，有：

$\left(\begin{array}{c}{x}_{n+1}\\ y_{n+1}\end{array}\right)=\left(\begin{array}{c}\cos \mathrm{\θ}\\ \sin \mathrm{\θ}\end{array}\right)\≈K\left(\begin{array}{c}1-\underset{N=1}{\overset{\mathrm{int}\left(\frac{n+1}{2}\right)}{\mathrm{\Σ}}}\left\{{(-1)}^N\right[\underset{i=0}{\overset{n-2N+2}{\mathrm{\Σ}}}...\underset{j=i+2N-1}{\overset{n}{\mathrm{\Σ}}}(\frac{d_i}{2^i}...\frac{d_j}{2^j})\left]\right\}\\ \underset{N=0}{\overset{\mathrm{int}\left(\frac{n+1}{2}\right)}{\mathrm{\Σ}}}\left\{{(-1)}^N\right[\underset{i=0}{\overset{n+1-2N}{\mathrm{\Σ}}}...\underset{j=i+2N}{\overset{n}{\mathrm{\Σ}}}(\frac{d_i}{2^i}...\frac{d_j}{2^j})\left]\right\}\end{array}\right)---\left(5\right)$

由CORDIC算法通项公式可看出，计算转化为旋转方向d_i的移位、加法运算，若能先确定每次旋转方向，则可省去每次判断，用一个时钟即可同时计算任何一项迭代结果，非常适合FPGA实现。若先确定旋转次数n，借助z_i计算，d_i就可确定。

步骤20采用补码法判断正余弦CORDIC算法坐标旋转方向，根据坐标旋转方向求出旋转坐标量，获取旋转坐标运算结果；

（1）判断下次旋转方向d_i，要根据角度变量z_i大小判断，设固定旋转角度位宽为w，z_i符号位以及扩展到w位分别表示为C_zi、|C_zi|_w，第i次旋转固定角度的原码、补码及反码分别表示为那么有：

①对于z_i>0情况，C_zi=0，需顺时针旋转再通过求出下次z_i+1，由|C_zi|_w与异或求得；

②对于z_i<0情况，C_zi=1，需逆时针旋转由求出z_i+1，由先|C_zi|_w与异或得再加C_zi得到。

综合以上分析，可得如下表达式：

其物理意义是：在先知道一系列以及初始值z₀、迭代次数n情况下，可借助求z_i+1及C_zi+1，获得下次d_i+1，通过一条数据链可连续求出n位旋转方向d₁、d₂、…、d_n。

（2）设坐标(x_i、y_i)位宽w，d_i符号位以及扩展到w位分别表示为C_di、|C_di|_w；x_i、y_i右移i位后分别表示为x_i(>>i)、y_i(>>i)，它们的反码、补码分别表示为[x_i(>>i)]_反和[y_i(>>i)]_反、[x_i(>>i)]_补和[y_i(>>i)]_补，那么式(5)计算可用下式实现，这里省略k_i，在n足够大下，校模因子用K=0.607253统一处理。

$\left\{\begin{array}{c}{x}_{i+1}=x_i-[{\left|C_{\mathrm{di}}\right|}_w\&CirclePlus;y_i(>>i)+C_{\mathrm{di}}]\\ y_{i+1}=y_i+[{\left|C_{\mathrm{di}}\right|}_w\&CirclePlus;x_i(>>i)+C_{\mathrm{di}}]\end{array}\right.---\left(7\right)$

通过一条数据链可连续求出n次旋转坐标(x₁、y₁)、(x₂、y₂)…、(x_n+1、y_n+1)即可获得运算结果x_n+1＝cosθ，y_n+1＝sinθ。

步骤30通过对正余弦CORDIC算法的仿真实验，对旋转坐标运算结果进行分析，如图3仿真试验的实现步骤如下：

（1）初始化旋转坐标变量初值和初始角度值，确定输入角度位宽、迭代次数等参数。

（2）输入角度范围转换和数据浮定点格式转换，根据输入角度的大小，转化为到之间，如果则如果 则如果则如果则θ_d＝θ，最终转化为到的角度值；并根据位宽把角度浮点数转化为定点数。本算法采用13位位宽，则以4096做为基数表示1，则角度θ_d＝π*4096*θ/180。

（3）应用补码法求解角度变量z_i+1，确定旋转方向d_i，求解旋转坐标变量(x_i+1,y_i+1)，按公式连续迭代n次，角度变量z_i＝0，输出待求角度的正弦和余弦值，即x_n+1＝cosθ，y_n+1＝sinθ。

（4）通过仿真，对计算速度、初始时延、存储资源和运算精度进行分析。

计算速度分析，如图4所示，设置时钟频率为100MHz，迭代次数为8次，输入角度采用13位位宽，添加测试文件，输入角度θ在为0°～360°，添加cosθ输出作为观察信号(sinθ类似)，波形显示方式。可以看出，通项法在采取仅流水线所需1/6左右寄存器情况下，通项法、流水线输出波形几乎同步，而迭代法则慢得非常多。

初始时延分析，参见图5，输出直接采用数字表示方式，选择θ＝0°作为观察对象，图中4097处表示cosθ=1、时钟周期为T。可以看出，通项法的初始时延3T明显优于传统迭代法、流水线法的初始时延11T。

存储资源分析，参见表1，输入角度采用13位位宽迭代法、流水线、通项法所需的逻辑单元(ALUTs)以及寄存器（registers)分别为330、369、370和168、322、52，这三种方法所需逻辑单元比较接近，但通项法使用寄存器要少得多，仅流水线所需的1/6左右，硬件资源明显减少。

表1

w=13下硬件开销、初始时延比较表

运算精度分析参见表2，采用13位运算精度，量化误差为1/4096，cosθ迭代值与cosθ函数真值的差值最大为0.0076，约0.8%的计算误差，在一定要求下，这个精度是可以接受的。

表2

任意角度输入cosθ运算值与计算值的比较表

虽然本发明所揭露的实施方式如上，但所述的内容只是为了便于理解本发明而采用的实施方式，并非用以限定本发明。任何本发明所属技术领域内的技术人员，在不脱离本发明所揭露的精神和范围的前提下，可以在实施的形式上及细节上作任何的修改与变化，但本发明的专利保护范围，仍须以所附的权利要求书所界定的范围为准。

Claims (5)

1.一种应用补码方法的正余弦CORDIC算法在FPGA实现的方法，该包括：

A采用数学归纳法由CORDIC算法推导正余弦CORDIC算法；

B采用补码法判断正余弦CORDIC算法坐标旋转方向，根据坐标旋转方向求出旋转坐标量，获取旋转坐标运算结果；

C通过对正余弦CORDIC算法的仿真实验，对旋转坐标运算结果进行分析；

所述步骤B中采用补码法判断坐标旋转方向具体包括：

根据角度中间变量z_i大小判断，z_i符号位为C_zi，设固定旋转角度位宽为w，z_i扩展到w位时的符号位为|C_zi|_w，第i次旋转固定角度的原码、补码及反码分别表示为所述i＝0，1，2……n，其公式为：

由C_zi+1可获得下一次d_i+1，通过一条数据链可连续求出n位旋转方向d₁、d₂、…、d_n。

2.根据权利要求1所述的应用补码方法的正余弦CORDIC算法在FPGA实现的方法，其特征在于，所述步骤A具体包括：

在直角坐标系旋转向量图中，确定单位向量p₀(x₀，y₀)旋转角θ到p_n(x_n，y_n)的关系式；

根据所述单位向量p₀(x₀，y₀)旋转角θ到p_n(x_n，y_n)的关系式得到CORDIC公式，从而计算转化为旋转方向的移位和加法运算，确定正余弦CORDIC算法坐标旋转次数，通过角度变量确定旋转坐标方向。

3.根据权利要求1所述的应用补码方法的正余弦CORDIC算法在FPGA实现的方法，其特征在于，所述步骤B具体包括：

用一条数据链连续求出n位正余弦CORDIC算法坐标旋转方向；

通过补码法求解n次旋转坐标值，根据所述坐标值获得正余弦函数运算结果。

4.根据权利要求1所述的应用补码方法的正余弦CORDIC算法在FPGA实现的方法，其特征在于，所述步骤C具体包括：

初始化旋转坐标变量初值和初始角度值，确定输入角度位宽及迭代次数参数；

转换输入角度的范围，并根据位宽将角度浮点数转化为定点数；

应用补码法求解角度变量，确定坐标旋转方向，求解旋转坐标变量，输出待求角度的正弦和余弦值；

对待求角度的正弦和余弦值的计算速度、初始时延、存储资源和运算精度进行分析。

5.根据权利要求1所述的应用补码方法的正余弦CORDIC算法在FPGA实现的方法，其特征在于，所述旋转坐标量的求解具体包括：

设旋转坐标(x_i、y_i)位宽w，第i次旋转方向为d_i，d_i的符号位为C_di，d_i扩展到w位时的符号位为|C_di|_w，所述i＝0，1，2……n；x_i、y_i右移i位后分别表示为x_i(>>i)、y_i(>>i)，它们的反码、补码分别表示为[x_i(>>i)]_反和[y_i(>>i)]_反、[x_i(>>i)]_补和[y_i(>>i)]_补，其公式为：

$\left\{\begin{array}{c}{x}_{i+1}=x_i-\&lsqb;|C_{di}{|}_w\&CirclePlus;y_i(>>i)+C_{di}\&rsqb;\\ y_{i+1}=y_i+\&lsqb;|C_{di}{|}_w\&CirclePlus;x_i(>>i)+C_{di}\&rsqb;\end{array}\right.$

通过一条数据链可连续求出n次旋转坐标(x₁、y₁)、(x₂、y₂)…、(x_n+1、y_n+1〕即可获得运算结果x_n+1＝cosθ，y_n+1＝sinθ。