1、cordic算法详解cordic算法详解 转载自小一休哥的文章:目前,学习与开发FPGA的程序员们大多使用的是Verilog HDL语言(以下简称为Verilog),关于Verilog的诸多优点一休哥就不多介绍了,在此,我们将重点放在Verilog的运算操作上。 我们都知道,在Verilog中,运算一般分为逻辑运算(与或非等)与算术运算(加减乘除等)。而在一开始学习Verilog时,老司机一定会提醒我们,“切记,千万别用/除、%取模(有的也叫取余)和*幂。”这话说的不无道理,因为这三个运算是不可综合的。但,需清楚理解的是,不可综合的具体意思为不能综合为简单的模块,当我们在程序中调用了这些运算时
2、,/除和%取模在Quartus软件中是可以综合的,因此可以正常调用运行,但是会消耗一些逻辑资源,而且会产生延时,即这两个运算的处理时间会很长,可能会大于时序控制时钟的单周期时间。此时呢,我们会建议你调用IP核来实现运算操作,虽然这样也会消耗许多逻辑资源,但产生的延时相对较小满足了你基本的需求。 问题好像迎刃而解了,可是仔细一想,除了这些运算,我们还剩下什么?对呀,三角函数,反三角函数,对数函数,指数函数呢,这些函数我们在高中就学习了的呀,难道在FPGA中就没有用武之地吗?有人会说,查找表呗,首先将某个运算的所有可能的输入与输出对一一罗列出来,然后放进Rom中,然后根据输入查表得到输出。这个方法
3、虽然有效的避免了延时问题,却是一个十分消耗资源的方法,不适合资源紧张的设计。那么,就真的没有办法了吗? 答案就是咱们今天的标题了,CORDIC,而且CORDIC是一个比较全能的算法,通过这一原理,我们可以实现三角函数,反三角函数,对数函数,指数函数等多种运算。接下来,一休哥就带领大家来学习CORDIC的原理吧。(题外话:请相信一休哥,本文不会让你感到太多痛苦) 本文将分三个小部分来展开介绍: 1、CORDIC的基本原理介绍 2、CORDIC的具体操作流程介绍 3、CORDIC的旋转模式Verilog仿真 本文涉及到的全部资料链接: 链接: 密码:x92u一、CORDIC的基本原理介绍CORDI
4、C算法是一个“化繁为简”的算法,将许多复杂的运算转化为一种“仅需要移位和加法”的迭代操作。CORDIC算法有旋转和向量两个模式,分别可以在圆坐标系、线性坐标系和双曲线坐标系使用,从而可以演算出8种运算,而结合这8种运算也可以衍生出其他许多运算。下表展示了8种运算在CORDIC算法中实现的条件。 这里写图片描述 首先,我们先从旋转模式下的圆坐标系讲起,这也是CORDIC方法最初的用途。1、CORDIC的几何原理介绍假设在xy坐标系中有一个点P1(x1,y1),将P1点绕原点旋转角后得到点P2(x2,y2)。 这里写图片描述 于是可以得到P1和P2的关系。 x2 = x1cos y1sin = c
5、os(x1 y1tan) y2 = y1cos + x1sin = cos(y1 +x1tan) 以上就是CORDIC的几何原理部分,而我们该如何深入理解这个几何原理的真正含义呢? 从原理中,我们可以知道,当已知一个点P1的坐标,并已知该点P1旋转的角度,则可以根据上述公式求得目标点P2的坐标。然后,麻烦来了,我们需要用FPGA去执行上述运算操作,而FPGA的Verilog语言根本不支持三角函数运算。因此,我们需要对上述式子进行简化操作,将复杂的运算操作转换为一种单一的“迭代位移”算法。那么,接下来我们介绍优化算法部分。2、CORDIC的优化算法介绍我们先介绍算法的优化原理:将旋转角细化成若干
6、分固定大小的角度i,并且规定i满足tani = 2-i,因此i的值在-99.7,99.7范围内,如果旋转角超出此范围,则运用简单的三角运算操作即可(加减)。 然后我们需要修改几何原理部分的假设,假设在xy坐标系中有一个点P0(x0,y0),将P0点绕原点旋转角后得到点Pn(xn,yn)。 于是可以得到P0和Pn的关系。 xn = x0cos y0sin = cos(x0 y0tan) yn = y0cos + x0sin = cos(y0 + x0tan) 然后,我们将旋转角细化成i,由于每次的旋转角度i是固定不变的(因为满足tani = 2-i),如果一直朝着一个方向旋转则i一定会超过(如果
7、在-99.7,99.7范围内)。因此我们需要对i设定一个方向值di。如果旋转角已经大于,则di为-1,表示下次旋转为顺时针,即向靠近;如果旋转角已经小于,则di为1,表示下次旋转为逆时针,即也向靠近。然后我们可以得到每次旋转的角度值dii,设角度剩余值为zi+1,则有zi+1 = zi - dii,其中z0为。如此随着i的增大,角度剩余值zi+1将会趋近于0,此时运算结束。(注:可以发现,di与zi的符号位相同) 第一次旋转0,d0为旋转方向: x1 = cos0(x0 d0y0tan0) y1 = cos0(y0 + d0x0tan0) 第二次旋转1,d1为旋转方向: x2 = cos1(x
8、1 d1y1tan1) = cos1cos0(x0 d0y0tan0 d1y0tan1 d1d0 x0tan1 tan0) y2 = cos1(y1 + d1x1tan1) = cos1cos0(y0 + d0x0tan0 + d1x0tan1 d1d0y0tan1 tan0) 大家可能已经发现了,在我们旋转的过程中,式子里一直会有tani和cosi,而每次都可以提取出cosi。虽然我们的FPGA无法计算tani,但我们知道tani = 2-i,因此可以执行和tani效果相同的移位操作2-i来取代tani。而对于cosi,我们可以事先全部提取出来,然后等待迭代结束之后(角度剩余值zi+1趋近于
9、0,一般当系统设置最大迭代次数为16时zi+1已经很小了),然后计算出cosi的值即可。 总结一下: 迭代公式有三: xi+1 = xi d iy i2-i,提取了cosi,2-i等效替换了tani之后 yi+1 = yi + d ix i2-i,提取了cosi,2-i等效替换了tani之后 zi+1 = zi - dii 其中i从0开始迭代,假设当i = n-1时,zn趋近于0,迭代结束。然后对结果乘上cosi(i从0至n-1),于是得到点Pn(xncosi,yncosi),此时的点Pn就近似等于之前假设中的点Pn(xn,yn)了,所以此时的点Pn同样满足之前假设得到的公式: xncosi
10、= x0cos y0sin yncosi = y0cos + x0sin 由于i从0至n-1,所以上式可以转化成下式: xn = 1/cosi(x0cos y0sin),(其中i从0至n-1) yn = 1/cosi(y0cos + x0sin),(其中i从0至n-1) 注意:上式中的xn,yn是经过迭代后的结果,而不是之前假设中的点Pn(xn,yn)。了解这点是十分重要的。 根据高中学的三角函数关系,可以知道cosi = 1/(1+tan2i)0.5 = 1/(1+2-2i)0.5,而1/(1+2-2i)0.5的极值为1,因此我们可以得出一个结论:当i的次数很大时,cosi的值趋于一个常数。
11、 关于如何计算cosi的代码如下所示:close all;clear;clc;% 初始化die = 16;%迭代次数jiao = zeros(die,1);%每次旋转的角度cos_value = zeros(die,1);%每次旋转的角度的余弦值K = zeros(die,1);%余弦值的N元乘积K_1 = zeros(die,1);%余弦值的N元乘积的倒数for i = 1 : die a = 2(-(i-1); jiao(i) = atan(a); cos_value(i) = cos(jiao(i); if( i = 1) K(i) = cos_value(i); K_1(i) = 1/
12、K(i); else K(i) = K(i-1)*cos_value(i); K_1(i) = 1/K(i); endendjiao = vpa(rad2deg(jiao)*256,10) cos_value = vpa(cos_value,10)K = vpa(K,10)K_1 = vpa(K_1,10)12345678910111213141516171819202122232425这里写图片描述 从上表也可以看出,当迭代次数为16,i=15时,cosi的值已经非常趋近于1了,cosi的值则约等于0.607253,1/cosi为1.64676。所以当迭代次数等于16时,通过迭代得到的点Pn
13、坐标已经非常接近之前假设中的点Pn。所以,当迭代次数等于16时,这个式子是成立的。 xn = 1/cosi(x0cos y0sin),(其中i从0至n-1) yn = 1/cosi(y0cos + x0sin),(其中i从0至n-1) 此时,已知条件有三个x0、y0和。通过16次迭代,我们可以得到xn和yn。而式中的cosi是个随i变化的值,我们可以预先将其值存入系统中。 然后,我们人为设置x0 = cosi,y0 = 0,则根据等式,xn = cos,yn = sin。其中1/cosi的值我们也同样预先存入系统中。如此,我们就实现了正弦和余弦操作了。二、CORDIC的具体操作流程介绍1、CO
14、RDIC的旋转模式由于算法较复杂,一休哥再总结一些具体的操作流程。 1、 设置迭代次数为16,则x0 = 0.607253,y0 = 0,并输入待计算的角度,在-99.7,99.7范围内。 2、 根据三个迭代公式进行迭代,i从0至15: xi+1 = xi d iy i2-i yi+1 = yi + d ix i2-i zi+1 = zi - dii 注:z0 = ,di与zi同符号。 3、 经过16次迭代计算后,得到的x16 和y16分别为cos和sin。 至此,关于CORDIC的三角函数cos和sin的计算原理讲解结束。 关于CORDIC算法计算三角函数cos和sin的MATLAB代码如下
15、所示:close all;clear;clc;% 初始化die = 16;%迭代次数x = zeros(die+1,1);y = zeros(die+1,1);z = zeros(die+1,1);x(1) = 0.607253;%初始设置z(1) = pi/4;%待求角度%迭代操作for i = 1:die if z(i) = 0 d = 1; else d = -1; end x(i+1) = x(i) - d*y(i)*(2(-(i-1); y(i+1) = y(i) + d*x(i)*(2(-(i-1); z(i+1) = z(i) - d*atan(2(-(i-1);endcosa
16、= vpa(x(17),10)sina = vpa(y(17),10)c = vpa(z(17),10)1234567891011121314151617181920212223242、CORDIC的向量模式讲完了旋转模式后,我们接着讲讲向量模式下的圆坐标系。 在这里,我们需从头来过了,假设在xy坐标系中有一个点P0(x0,y0),将P0点绕原点旋转角后得到点Pn(xn,0),在-99.7,99.7范围内。 于是可以得到P0和Pn的关系: xn = x0cos y0sin = cos(x0 y0tan) yn = y0cos + x0sin = cos(y0 + x0tan) = 0 如何得到
17、Pn(xn,yn)一直是我们的目标。而此时,我们还是列出那几个等式: 根据三个迭代公式进行迭代,i从0至15: xi+1 = xi d iy i2-i yi+1 = yi + d ix i2-i zi+1 = zi - dii 不过此时我们尝试改变初始条件: 设置迭代次数为16,则x0 = x,y0 = y,z0 = 0,di与yi的符号相反。表示,经过n次旋转,使Pn靠近x轴。 因此,当迭代结束之后,Pn将近似接近x轴,此时yn = 0,可知旋转了,即zn = = arctan(y/x)。 而 xn = 1/cosi(x0cos y0sin),(其中i从0至n-1) yn = 1/cosi(
18、y0cos + x0sin),(其中i从0至n-1) 因此,可得ycos + xsin = 0, xn = 1/cosi(xcos ysin) = 1/cosi (xcos ysin)2(1/2) = 1/cosi x2cos2 + y2sin2 2xysincos(1/2) = 1/cosi x2cos2 + y2sin2 + y2 cos2 + x2sin2(1/2) = 1/cosi x2 + y2(1/2) 由上可以知道,我们通过迭代,就算出了反正切函数zn = = arctan(y/x),以及向量OP0(x,y)的长度 d = xn * cosi。 关于反正切函数,一休哥要多啰嗦几句
19、了,由于在-99.7,99.7范围内,所以我们输入向量OP0(x,y)时,需要保证其在第一、四象限。 关于CORDIC算法计算反三角函数arctan的MATLAB代码如下所示:close all;clear;clc;% 初始化die = 16;%迭代次数x = zeros(die+1,1);y = zeros(die+1,1);z = zeros(die+1,1);x(1) = 100;%初始设置y(1) = 200;%初始设置k = 0.607253;%初始设置%迭代操作for i = 1:die if y(i) = 0 d = -1; else d = 1; end x(i+1) = x(
20、i) - d*y(i)*(2(-(i-1); y(i+1) = y(i) + d*x(i)*(2(-(i-1); z(i+1) = z(i) - d*atan(2(-(i-1);endd = vpa(x(17)*k,10)a = vpa(y(17),10)c = vpa(rad2deg(z(17),10)12345678910111213141516171819202122232425三、CORDIC的旋转模式Verilog仿真一休哥在编写CORDIC算法时,采用了16级流水线,仿真效果十分明显。以下是顶层文件的代码。 为了避免浮点运算,为了满足精度要求,一休哥对每个变量都放大了216倍,并且
21、引入了有符号型reg和算术右移。 关于Verilog代码的编写,一休哥已经不想多说了,因为代码是完全符合我之前所讲的CORDIC的原理与MATLAB仿真代码。相信大家在看完本文的前两个部分之后,对Verilog的理解应该不是难事儿。module Cordic_Test( CLK_50M,RST_N, Phase, Sin,Cos,Error);input CLK_50M;input RST_N;input 31:0 Phase;output 31:0 Sin;output 31:0 Cos;output 31:0 Error;define rot0 32d2949120 /45度*216def
22、ine rot1 32d1740992 /26.5651度*216define rot2 32d919872 /14.0362度*216define rot3 32d466944 /7.1250度*216define rot4 32d234368 /3.5763度*216define rot5 32d117312 /1.7899度*216define rot6 32d58688 /0.8952度*216define rot7 32d29312 /0.4476度*216define rot8 32d14656 /0.2238度*216define rot9 32d7360 /0.1119度*21
23、6define rot10 32d3648 /0.0560度*216define rot11 32d1856 /0.0280度*216define rot12 32d896 /0.0140度*216define rot13 32d448 /0.0070度*216define rot14 32d256 /0.0035度*216define rot15 32d128 /0.0018度*216parameter Pipeline = 16;parameter K = 32h09b74; /K=0.607253*216,32h09b74,reg signed 31:0 Sin;reg signed 3
24、1:0 Cos;reg signed 31:0 Error;reg signed 31:0 x0=0,y0=0,z0=0;reg signed 31:0 x1=0,y1=0,z1=0;reg signed 31:0 x2=0,y2=0,z2=0;reg signed 31:0 x3=0,y3=0,z3=0;reg signed 31:0 x4=0,y4=0,z4=0;reg signed 31:0 x5=0,y5=0,z5=0;reg signed 31:0 x6=0,y6=0,z6=0;reg signed 31:0 x7=0,y7=0,z7=0;reg signed 31:0 x8=0,y
25、8=0,z8=0;reg signed 31:0 x9=0,y9=0,z9=0;reg signed 31:0 x10=0,y10=0,z10=0;reg signed 31:0 x11=0,y11=0,z11=0;reg signed 31:0 x12=0,y12=0,z12=0;reg signed 31:0 x13=0,y13=0,z13=0;reg signed 31:0 x14=0,y14=0,z14=0;reg signed 31:0 x15=0,y15=0,z15=0;reg signed 31:0 x16=0,y16=0,z16=0;reg 1:0 Quadrant Pipel
26、ine:0;always (posedge CLK_50M or negedge RST_N)begin if(!RST_N) begin x0 = 1b0; y0 = 1b0; z0 = 1b0; end else begin x0 = K; y0 = 32d0; z0 = Phase15:0 16; endendalways (posedge CLK_50M or negedge RST_N)begin if(!RST_N) begin x1 = 1b0; y1 = 1b0; z1 = 1b0; end else if(z031) begin x1 = x0 + y0; y1 = y0 -
27、 x0; z1 = z0 + rot0; end else begin x1 = x0 - y0; y1 = y0 + x0; z1 = z0 - rot0; endendalways (posedge CLK_50M or negedge RST_N)begin if(!RST_N) begin x2 1); z2 1); z2 = z1 - rot1; endendalways (posedge CLK_50M or negedge RST_N)begin if(!RST_N) begin x3 2); z3 2); z3 = z2 - rot2; endendalways (posedg
28、e CLK_50M or negedge RST_N)begin if(!RST_N) begin x4 3); z4 3); z4 = z3 - rot3; endendalways (posedge CLK_50M or negedge RST_N)begin if(!RST_N) begin x5 4); z5 4); z5 = z4 - rot4; endendalways (posedge CLK_50M or negedge RST_N)begin if(!RST_N) begin x6 5); z6 5); z6 = z5 - rot5; endendalways (posedg
29、e CLK_50M or negedge RST_N)begin if(!RST_N) begin x7 6); z7 6); z7 = z6 - rot6; endendalways (posedge CLK_50M or negedge RST_N)begin if(!RST_N) begin x8 7); z8 7); z8 = z7 - rot7; endendalways (posedge CLK_50M or negedge RST_N)begin if(!RST_N) begin x9 8); z9 8); z9 = z8 - rot8; endendalways (posedge CLK_50M or
copyright@ 2008-2022 冰豆网网站版权所有
经营许可证编号:鄂ICP备2022015515号-1