cordic算法详解Word文档格式.docx

1、1、CORDIC的几何原理介绍假设在xy坐标系中有一个点P1（x1，y1），将P1点绕原点旋转角后得到点P2（x2，y2）。于是可以得到P1和P2的关系。x2 = x1cos y1sin = cos（x1 y1tan） y2 = y1cos + x1sin = cos（y1 +x1tan） 以上就是CORDIC的几何原理部分，而我们该如何深入理解这个几何原理的真正含义呢？从原理中，我们可以知道，当已知一个点P1的坐标，并已知该点P1旋转的角度，则可以根据上述公式求得目标点P2的坐标。然后，麻烦来了，我们需要用FPGA去执行上述运算操作，而FPGA的Verilog语言根本不支持三角函数运算。因此

2、，我们需要对上述式子进行简化操作，将复杂的运算操作转换为一种单一的“迭代位移”算法。那么，接下来我们介绍优化算法部分。2、CORDIC的优化算法介绍我们先介绍算法的优化原理：将旋转角细化成若干分固定大小的角度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 +

3、x0sin = cos（y0 + x0tan） 然后，我们将旋转角细化成i，由于每次的旋转角度i是固定不变的（因为满足tani = 2-i），如果一直朝着一个方向旋转则i一定会超过（如果在-99.7范围内）。因此我们需要对i设定一个方向值di。如果旋转角已经大于，则di为-1，表示下次旋转为顺时针，即向靠近；如果旋转角已经小于，则di为1，表示下次旋转为逆时针，即也向靠近。然后我们可以得到每次旋转的角度值dii，设角度剩余值为zi+1，则有zi+1 = zi - dii，其中z0为。如此随着i的增大，角度剩余值zi+1将会趋近于0，此时运算结束。（注：可以发现，di与zi的符号位相同） 第一次

4、旋转0，d0为旋转方向：x1 = cos0（x0 d0y0tan0） y1 = cos0（y0 + d0x0tan0） 第二次旋转1，d1为旋转方向：x2 = cos1（x1 d1y1tan1） = cos1cos0（x0 d0y0tan0 d1y0tan1 d1d0 x0tan1 tan0） y2 = cos1（y1 + d1x1tan1） = cos1cos0（y0 + d0x0tan0 + d1x0tan1 d1d0y0tan1 tan0） 大家可能已经发现了，在我们旋转的过程中，式子里一直会有tani和cosi，而每次都可以提取出cosi。虽然我们的FPGA无法计算tani，但我们知道

5、tani = 2-i，因此可以执行和tani效果相同的移位操作2-i来取代tani。而对于cosi，我们可以事先全部提取出来，然后等待迭代结束之后（角度剩余值zi+1趋近于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）

6、，于是得到点Pn（xncosi，yncosi），此时的点Pn就近似等于之前假设中的点Pn（xn，yn）了，所以此时的点Pn同样满足之前假设得到的公式：xncosi = 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+t

7、an2i）0.5 = 1/（1+2-2i）0.5，而1/（1+2-2i）0.5的极值为1，因此我们可以得出一个结论：当i的次数很大时，cosi的值趋于一个常数。关于如何计算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） = a

8、tan（a）; cos_value（i） = cos（jiao（i）; if（ i = 1） K（i） = cos_value（i）; K_1（i） = 1/K（i）; else K（i） = K（i-1）*cos_value（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了，cos

9、i的值则约等于0.607253，1/cosi为1.64676。所以当迭代次数等于16时，通过迭代得到的点Pn坐标已经非常接近之前假设中的点Pn。所以，当迭代次数等于16时，这个式子是成立的。此时，已知条件有三个x0、y0和。通过16次迭代，我们可以得到xn和yn。而式中的cosi是个随i变化的值，我们可以预先将其值存入系统中。然后，我们人为设置x0 = cosi，y0 = 0，则根据等式，xn = cos，yn = sin。其中1/cosi的值我们也同样预先存入系统中。如此，我们就实现了正弦和余弦操作了。二、CORDIC的具体操作流程介绍1、CORDIC的旋转模式由于算法较复杂，一休哥再总结一

10、些具体的操作流程。1、 设置迭代次数为16，则x0 = 0.607253，y0 = 0，并输入待计算的角度，在-99.7范围内。2、 根据三个迭代公式进行迭代，i从0至15：xi+1 = xi d iy i2-i yi+1 = yi + d ix i2-i 注：z0 = ，di与zi同符号。3、 经过16次迭代计算后，得到的x16 和y16分别为cos和sin。至此，关于CORDIC的三角函数cos和sin的计算原理讲解结束。关于CORDIC算法计算三角函数cos和sin的MATLAB代码如下所示：x = zeros（die+1,1）;y = zeros（die+1,1）;z = zeros（

11、die+1,1）;x（1） = 0.607253;%初始设置z（1） = pi/4;%待求角度%迭代操作for i = 1:die if z（i） = 0 d = 1; d = -1; 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）;cosa = vpa（x（17）,10）sina = vpa（y（17）,10）c = vpa（z（17）,10）1234567891011121314151617181920212223242、CORDIC的向量

12、模式讲完了旋转模式后，我们接着讲讲向量模式下的圆坐标系。在这里，我们需从头来过了，假设在xy坐标系中有一个点P0（x0，y0），将P0点绕原点旋转角后得到点Pn（xn，0），在-99.7于是可以得到P0和Pn的关系：yn = y0cos + x0sin = cos（y0 + x0tan） = 0 如何得到Pn（xn，yn）一直是我们的目标。而此时，我们还是列出那几个等式：根据三个迭代公式进行迭代，i从0至15：不过此时我们尝试改变初始条件：设置迭代次数为16，则x0 = x，y0 = y，z0 = 0，di与yi的符号相反。表示，经过n次旋转，使Pn靠近x轴。因此，当迭代结束之后，Pn将近似接

13、近x轴，此时yn = 0，可知旋转了，即zn = = arctan（y/x）。而 因此，可得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。关于反正

14、切函数，一休哥要多啰嗦几句了，由于在-99.7范围内，所以我们输入向量OP0（x，y）时，需要保证其在第一、四象限。关于CORDIC算法计算反三角函数arctan的MATLAB代码如下所示：x（1） = 100;y（1） = 200;k = 0.607253; if y（i） d = vpa（x（17）*k,10）a = vpa（y（17）,10）c = vpa（rad2deg（z（17）,10）12345678910111213141516171819202122232425三、CORDIC的旋转模式Verilog仿真一休哥在编写CORDIC算法时，采用了16级流水线，仿真效果十分明显。以下

15、是顶层文件的代码。为了避免浮点运算，为了满足精度要求，一休哥对每个变量都放大了216倍，并且引入了有符号型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;0 Cos;0 Error;define ro

16、t0 32d2949120 /45度*216define 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

17、 rot9 32d7360 /0.1119度*216define 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 sig

18、ned 31:0 x0=0,y0=0,z0=0;0 x1=0,y1=0,z1=0;0 x2=0,y2=0,z2=0;0 x3=0,y3=0,z3=0;0 x4=0,y4=0,z4=0;0 x5=0,y5=0,z5=0;0 x6=0,y6=0,z6=0;0 x7=0,y7=0,z7=0;0 x8=0,y8=0,z8=0;0 x9=0,y9=0,z9=0;0 x10=0,y10=0,z10=0;0 x11=0,y11=0,z11=0;0 x12=0,y12=0,z12=0;0 x13=0,y13=0,z13=0;0 x14=0,y14=0,z14=0;0 x15=0,y15=0,z15=0;0 x

19、16=0,y16=0,z16=0;reg 1:0 Quadrant Pipeline:0;always （posedge CLK_50M or negedge RST_N）begin if（!RST_N） begin x0 = 1b0; y0 z0 = K;= 32d0;= Phase15:0 16; x1 y1 z1 else if（z031）= x0 + y0;= y0 - x0;= z0 + rot0;= x0 - y0;= y0 + x0;= z0 - rot0; x2 1）; z2 1）; z2 = z1 - rot1; x3 2）; z3 2）; z3 = z2 - rot2; x4 3）; z4 3）; z4 = z3 - rot3; x5 4）; z5 4）; z5 = z4 - rot4; x6 5）; z6 5）; z6 = z5 - rot5; x7 6）; z7 6）; z7 = z6 - rot6; x8 7）; z8 7）; z8 = z7 - rot7; x9 8）; z9 8）; z9 = z8 - rot8;always （posedge CLK_50M or