CORDIC(Coordinate Rotation Digital Computer)算法即坐标旋转数字计算方法,是J.D.Volder1于1959年首次提出,主要用于三角函数、双曲线、指数、对数的计算。该算法通过基本的加和移位运算代替乘法运算,使得矢量的旋转和定向的计算不再需要三角函数、乘法、开方、反三角、指数等函数,计算向量长度并能把直角坐标系转换为极坐标系。因为Cordic 算法只用了移位和加法,很容易用纯硬件来实现,非常适合FPGA实现。
CORDIC算法完成坐标或向量的平面旋转(下图以逆时针旋转为例)。
旋转后,可得如下向量:
旋转的角度θ经过多次旋转得到的(步步逼近,接近二分查找法),每次旋转一小角度。单步旋转定义如下公式:
公式(2)提取cosθ,可修改为:
修改后的公式把乘法次数从4次改为3次,剩下的乘法运算可以通过选择每次旋转的角度去除,将每一步的正切值选为2的指数(二分查找法),除以2的指数可以通过右移操作完成(verilog)。
每次旋转的角度可以表示为:
所有迭代角度累加值等于最终需要的旋转角度θ:
这里Sn为1或者-1,根据旋转方向确定(后面有确定方法,公式(15)),顺时针为-1,逆时针为1。
可以得到如下公式:
结合公式(3)和(7),得到公式(8):
到这里,除了余弦值这个系数,算法只要通过简单的移位和加法操作完成。而这个系数可以通过预先计算最终值消掉。首先重新重写这个系数如下:
第二步计算所有的余弦值并相乘,这个值K称为增益系数。
由于K值是常量,我们可以先忽略它。
到这里我们发现,算法只剩下移位和加减法,这就非常适合硬件实现了,为硬件快速计算三角函数提供了一种新的算法。在进行迭代运算时,需要引入一个新的变量Z,表示需要旋转的角度θ中还没有旋转的角度。
这里,我们可以把前面提到确定旋转方向的方法介绍了,就是通过这个变量Z的符号确定。
通过公式(5)和(15),将未旋转的角度变为0。
一个类编程风格的结构如下,反正切值是预先计算好的。
旋转模式下,CORDIC算法驱动Z变为0,结合公式(13)和(16),算法的核心计算如下:
一种特殊情况是,另初始值如下:
因此,旋转模式下CORDIC算法可以计算一个输入角度的正弦值和余弦值。
向量模式下,有两种特例:
因此,向量模式下CORDIC算法可以用来计算输入向量的模和反正切,也能开方计算,并可以将直角坐标转换为极坐标。
算法介绍:http://en.wikipedia.org/wiki/Cordic,http://blog.csdn.net/liyuanbhu/article/details/8458769
根据算法原理,利用维基百科中给的程序,在matlab中跑了一遍,对算法有了一定程度的了解。
程序如下:
实现主要参考了相关作者的代码,然后对其进行了修改,最终实现了16级的流水线,设计完成旋转模式下正弦值和余弦值的计算。
http://www.cnblogs.com/qiweiwang/archive/2010/07/28/1787021.html,http://www.amobbs.com/forum.php?mod=viewthread&tid=5513050&highlight=cordic
下面分段介绍下各部分代码:
首先是角度的表示,进行了宏定义,360读用16位二进制表示2^16,每一度为2^16/360。
1 //360°--2^16,phase_in = 16bits (input [15:0] phase_in) 2 //1°--2^16/360 3 `define rot0 16'h2000 //45 4 `define rot1 16'h12e4 //26.5651 5 `define rot2 16'h09fb //14.0362 6 `define rot3 16'h0511 //7.1250 7 `define rot4 16'h028b //3.5763 8 `define rot5 16'h0145 //1.7899 9 `define rot6 16'h00a3 //0.895210 `define rot7 16'h0051 //0.447611 `define rot8 16'h0028 //0.223812 `define rot9 16'h0014 //0.111913 `define rot10 16'h000a //0.056014 `define rot11 16'h0005 //0.028015 `define rot12 16'h0003 //0.014016 `define rot13 16'h0001 //0.007017 `define rot14 16'h0001 //0.003518 `define rot15 16'h0000 //0.0018
然后是流水线级数定义、增益放大倍数以及中间结果位宽定义。流水线级数16,为了满足精度要求,有文献指出流水线级数必须大于等于角度位宽16(针对正弦余弦计算的CORDIC算法优化及其FPGA实现)。增益放大2^16,为了避免溢出状况中间结果(x,y,z)定义为17为,最高位作为符号位判断,1为负数,0为正数。
还需要定义memory型寄存器数组并初始化为0,用于寄存输入角度高2位的值。
接着,是对输入角度象限处理,将角度都转换到第一象限,方便处理。输入角度值最高两位赋值0,即转移到第一象限[0°,90°]。此外,完成x0,y0和z0的初始化,并增加一位符号位。
接下来根据剩余待旋转角度z的符号位进行16次迭代处理,即完成16级流水线处理。迭代公式:x(n+1) <= x(n) + {{n{y(n)[16]}},y(n)[16:n]},n为移位个数。右移之后高位补位,这里补位有一些不理解。移位可能存在负数,且没有四舍五入。按理说第一象限不存在负数,但后续仿真汇总确实有负数出现,但仿真结果良好。
由于进行了象限的转换,最终流水结果需要根据象限进行转换为正确的值。这里寄存17次高2位角度输入值,配合流水线结果用于象限判断,并完成转换。
最后,根据寄存的高2位角度输入值,利用三角函数关系,得出最后的结果,其中负数进行了补码操作。
仿真结果应该还是挺理想的。后续需要完成的工作:1.上述红色出现的问题的解决;2.应用cordic算法,完成如FFT的算法。
后记:
在3中,迭代公式:x(n+1) <= x(n) + {{n{y(n)[16]}},y(n)[16:n]},上述右移操作都是手动完成:首先最高位增加1位符号位(1为负,0为正),然后手动添加n位符号位(最高位)补齐,即实际上需要完成的算术右移(>>>)。本设计定义的reg为无符号型,在定义时手动添加最高位为符号位。verilog-1995中只有integer为有符号型,reg和wire都是无符号型,只能手动添加扩展位实现有符号运算。而在verilog-2001中reg和wire可以通过保留字signed定义为有符号型。另外,涉及有符号和无符号型的移位操作等可参考下面的文章。
verilog有符号数详解:http://www.cnblogs.com/LJWJL/p/3481995.html,
Verilog-2001新特性及代码实现:http://www.asic-world.com/verilog/verilog2k.html,
逻辑移位与算术移位区别:http://www.cnblogs.com/yuphone/archive/2010/09/21/1832217.html,http://blog.sina.com.cn/s/blog_65311d330100ij9n.html
原来的算法实现针对Verilog-1995中reg和wire没有有符号型,也没有verilog-2001中的算术移位而实现的。根据verilog-2001新特性,引入有符号型reg和算术右移,同样实现了前文的结果。代码如下:
另外,代码中可以适当优化下:1.流水线操作时,定义的中间寄存器在定义是可以选择memory型,且可以单独建立module或者task进行封装迭代过程; 2. 最后对高2位角度寄存时,可以利用for语句选择移位寄存器实现,如下所示。
疑问:有一点比较奇怪的是,转移到第一象限后,x和y不该存在负数的情况,但是现在确实有,这一点比较费解,所以将算术右移改为逻辑右移,在函数极值时存在错误。
联系客服