1、CORDIC算法
CORDIC（Coordinate Rotation DIgital Computer）是一种经过迭代实现快速平面旋转的算法，经过变形扩展，它能够对多种数学函数求值，例如求解三角/反三角函数、双曲函数等，除此以外还能够作特殊的向量运算甚至算术运算，这些内容不在本文讨论的重点之中。git
在CORDIC被发明以前，要对特殊函数求值，最天然的方法即是级数展开，例如利用泰勒展开来逼近目标函数，只要阶数取得足够大，就能够无限逼近目标函数。级数展开在数学上是完美的，但运用到计算机时，咱们很快就会发现问题：级数展开本质是用多项式函数来近似目标函数，这其中包括大量复杂浮点运算，对于没有硬件浮点运算单元的平台，只能经过软件浮点实现。github
CORDIC的出现解决了这个问题。该算法利用迭代逼近的方法，仅仅经过加/减和移位操做便可求解，极大的方便了计算机实现。算法
本文所作的工做：
函数
1、从CORDIC算法正向角度分析三角函数（sin, cos, tan）对任意角的求值；工具
2、从CORDIC算法逆向角度分析反三角函数（arcsin、arccos、arctan）对任意角的求值、向量模的求值；优化
3、分析定点运算的实现及软件模型的创建spa
4、经过Verilog HDL设计硬件。设计
本文代码仓库详见：https://github.com/cassuto/CORDIC-all-in-one-verilog 3d
2、核心思想
正如该算法的名字所说，CORDIC最初是为一种用来进行坐标轴旋转的专用计算机开发的（其原型硬件于1959年成功应用于实时导航）。追根溯源，CORDIC的核心就是坐标轴旋转，算法最先的出发点是解决旋转问题。它与函数求值究竟有什么关系？早期研究者是如何想到这个算法的呢？这正是本文试图说明的问题，固然并不完善，欢迎批评指正。调试
一、正向角度（已知旋转角，求点旋转后的坐标）
假设如今咱们要将一个直角坐标系中的点P(x0, y0)绕原点逆时针旋转z0角度，则变换后的点P1(x1, y1)坐标以下：
咱们发现这个转换式中涉及三角函数，但如今假设机器还不能求任意三角函数值，那么可否改进？
能够考虑查表实现，但查表法要求数据必须是离散的，这样旋转角只能取有限个值。如何对任意的旋转角求解？
简单，利用迭代逼近，把目标角度分红若干个小角度，每次迭代只旋转一个特定角度并靠近目标，经过若干次迭代后，点就被旋转到了近似的位置上。关键点是，经过特定的分割，使得每次旋转的角度都是特定角，若是将些特定角的三角函数值固化到列表中，就能够利用查表绕开三角函数的计算。对计算机而言虽然须要迭代屡次，但每次都是简单运算，整体速度快于复杂的三角求值。通常该算法的结果只是近似值，但经过设置迭代次数就能够控制精度，迭代次数越多，精度越高，并最终趋于稳定。
例如：将点P0旋转z = 50°到P‘，求P’坐标。
为了便于制做三角函数表，首先假设：对全部问题来讲，z的取值都是[0,90°]，所以假设每次旋转的角度都是90°的n分之一。
第一步正向旋转45°变换到P1，发现没有达到目标角度；所以又继续将P1正向旋转22.5°变换到P2，发现超出了目标角度；因而继续逆向旋转11.25°变换到P3，仍然超出目标角度；再逆向旋转5.625°，到达P4，这时若P4已经在设定精度内趋近答案，则该精度内能够用P4近似表示P’坐标。
这个算法基本上解决了旋转变换式中三角函数求值的问题，但这并非咱们想要的：每次迭代都须要进行4次浮点乘法运算，旋转角z的正余弦都涉及浮点数，例如cos 22.5° ≈ 0.923879。
继续观察旋转变换式，能够变造成：
上式仍然涉及两个三角函数，可是cos被放到了一边，因而重点研究tan。将角度z微分为n个小角度，若是像下面这样特殊选取z，使tan z刚好只与2的幂有关，则本来复杂的乘tan zn运算就能够经过右移n位实现（xn，yn都是整数，相关问题在后面讨论），这对二进制计算机是很是天然的。有人指出0.618黄金分割是更优的，但并不适合二进制计算机实现。
结合上面的变换式子，并把cosz隐藏到系数P中，咱们就有了递推关系（其中n表示当前迭代次数，n+1则表示下一次迭代。这里只写出了逆时针旋转的状况，顺时针旋转改变1/2n项的符号便可。）
因为n是离散的，而且z = atan(1/2n)能够提早计算，因此可用查表的方法快速得出z的值，这样系数Pn也能够经过查表求出，但系数Pn还是浮点数。
到这里，原先算法的4次浮点运算，被成功减小到了只有2次，咱们离真正的CORDIC算法已经很近了。
但问题尚未结束，若是迭代n次的话，仍须要进行2n次浮点运算，有没有优化的余地呢？引发2n次浮点运算的缘由是每次都须要与系数Pn相乘，这样作真的有必要吗？
构造辅助三角形，能够得出：
所以Pn的取值只与n有关，而与别的变量没有关系。Pn彻底能够在全部迭代都完成后单独计算，在最后将结果乘P=ΠPn便可。
对于根据迭代较少的状况，能够将P的不一样取值固化到表格中，经过查表快速求出。
而对于迭代次数较多的状况，能够将P看做常数近似处理，这是由于随着n的增长，P逐渐稳定下来：
经过数值统计，咱们能够看出具体多少次迭代之后开始能够将P看做常量。从图中能够看出从n=8开始比较合适。当n<3时P的取值变化较大，但事实上小于3次的迭代实用价值不大。
到这里整个算法只需在最后进行一次浮点运算，而每次迭代都只涉及简单的加/减和位移运算。
虽然有了这个快速旋转变换算法，但咱们的最终目标并不在于此，而是求三角函数的值，如何经过旋转求三角函数？
简单。借助单位圆这个工具，在直角坐标系中将点(0, 1)绕原点旋转至α角上，则旋转后的点的横纵坐标对应的正是cosα、sinα的值，tanα也能够利用比值求出，这正是利用CORDIC求三角函数的核心思路。下图显示了CORDIC算法的迭代过程。
上述算法可用伪代码描述以下：
1 For n=0 to i-1 2 If (Z(n) >= 0) Then 3 X(n + 1) := X(n) – (Yn >> n) 4 Y(n + 1) := Y(n) + (Xn >> n) 5 6 Z(n + 1) := Z(n) - atan(1/2^n) 7 Else 8 X(n + 1) := X(n) + (Yn >> n) 9 Y(n + 1) := Y(n) – (Xn >> n) 10 11 Z(n + 1) := Z(n) + atan(1/2^n) 12 End if 13 End for
二、逆向角度（已知点旋转后的点坐标，求旋转角）
从正向分析或者逆向分析，角度不一样，解决的问题也不一样。
仍然在单位圆内，假设已知正弦或余弦中的一个值，如何求对应的反三角函数？这里以求反正弦函数为例说明，已知反正弦函数自变量为S，咱们将点(1, 0)逆时针旋转，若是转到某个角度该点的纵坐标恰等于S，说明这个角度就是反正弦函数的值；同理反余弦函数也可用用相似的方法算出。
接下来继续分析反正切的状况，已知正切值，则能够得出对应的单位圆上的点P0(x0, y0)。若是将点P0旋转一个角度，使旋转后的点纵坐标变成0，那么这个角度正是反正切函数的值。
以上即是利用CORDIC求反三角函数的全部思路。
如今留下的问题是到底旋转多少角度才合适。CORDIC一样以迭代逼近的方法解决这个问题。
以求反正切函数为例：这里只讨论z在0到90°范围内（即x0, y0都为正）的状况。咱们能够先将原始点旋转atan(-1/2)角度进行试探，旋转完成后发现点纵坐标仍然大于0，因而继续旋转atan(-1/4)、atan(-1/8)到达图中草绿色终点，这个时候发现纵坐标小于0，因而再反方向旋转atan(-1/16)角度到达蓝色终点，若此时点的纵坐标在规定精度内接近0则迭代结束，能够经过累加全部迭代中旋转的角度增量求出α的具体数值，角度α正是atan的值。
经过第一小节的讨论能够看出，递推式中系数P在几何上的效果是缩放了旋转半径（点到坐标原点的距离）。因为求反正切函数过程当中只用到旋转角，而旋转角与旋转半径与无关，因此不须要考虑系数P。
上述算法可用伪代码描述以下：
1 A := 0 2 For n=0 to i-1 3 If (Y(n) >= 0) Then 4 X(n + 1) := X(n) + (Yn >> n) 5 Y(n + 1) := Y(n) - (Xn >> n) 6 7 A := A + atan(1/2^n) 8 Else 9 X(n + 1) := X(n) - (Yn >> n) 10 Y(n + 1) := Y(n) + (Xn >> n) 11 12 A := A - atan(1/2^n) 13 End if 14 End for
理解了求反正切函数的思路后，咱们就能够着手考虑反正弦和反余弦函数的求值。比较可知，旋转过程当中，对反正切函数要求旋转点的纵坐标趋于0，而对反正弦函数要求旋转点纵坐标趋于特定值S（函数自变量的取值）。能够发现，前者所涉及的的旋转是后者的一种特殊状况，这样只需简单修改前者就能够得出适用后者的旋转算法。
1 limP := 0.607253 2 A := 0 3 B := S / limP 4 For n=0 to i-1 5 If (Y(n) >= B) Then 6 X(n + 1) := X(n) + (Yn >> n) 7 Y(n + 1) := Y(n) - (Xn >> n) 8 9 A := A + atan(1/2^n) 10 Else 11 X(n + 1) := X(n) - (Yn >> n) 12 Y(n + 1) := Y(n) + (Xn >> n) 13 14 A := A - atan(1/2^n) 15 End if 16 End for
须要解释的是B := S / limP，这里为何须要将S放大1/limP倍呢？缘由很简单，因为在每次迭代时都忽略了系数Pn<1，这会致使Y(n)比未忽略系数时大出1/Pn倍。为了使Y(n)与S可以进行比较，须要将S同比例放大1/Pn倍，这样每次迭代又将引入浮点运算。直接近似Pn = limP虽然会引入细微的精度损失，但避免了大量浮点运算。
同理，也能够得出旋转点横坐标趋于特定值的算法。
1 limP := 0.607253 2 A := 0 3 B := S / limP 4 For n=0 to i-1 5 If (X(n) >= B) Then 6 X(n + 1) := X(n) - (Yn >> n) 7 Y(n + 1) := Y(n) + (Xn >> n) 8 9 A := A - atan(1/2^n) 10 Else 11 X(n + 1) := X(n) + (Yn >> n) 12 Y(n + 1) := Y(n) - (Xn >> n) 13 14 A := A + atan(1/2^n) 15 End if 16 End for
至此全部关于利用CORDIC求三角和反三角函数的分析就结束了。
借鉴上面求反正切角函数的思路，咱们也能够看出向量求模的运算方法。直角坐标系中，将目标向量V的起点平移到坐标原点，终点用坐标(x0, y0)表示。若是把这个点旋转到x轴（或y轴）上，则旋转后的点对应的横坐标（或纵坐标）正是该向量的模长|V|。
该算法的核心与求反正切函数基本相同，这里再也不赘述。
3、定点算法的软件模型创建
若是能用整数表示浮点数，则能够经过整数运算完成实数运算。在必定精度内，定点数与浮点数能够经过比例放缩互相转换。
首先肯定整数位宽，这里以16位为例，
一、量化角度：在直角范围内，能够将正整数的取值区间均分为coeff = 2^14 / 90 = 182个单位，则一个单位对应角度的数量级为10^-2。任意给定浮点角度Z°，其对应的定点角度为floor(Z * coeff)；
二、量化坐标(x, y)：将原始坐标放大coeff倍并取整便可得出定点坐标。但若坐标取值较大则存在溢出风险，需另选系数。
上述两点措施完成了算法输入的定点量化；而对于算法的输出，只需将定点数缩小coeff倍便可获得最终的浮点结果。
在进行硬件设计以前，有必要先考虑软件模型。由于软件语言的抽象程度高于硬件描述语言，能够更加紧凑地对算法进行描述，同时能快速地进行调试。
软件模型主要涉及两个核心算法：旋转和反旋转。
旋转算法的C语言模型：
https://github.com/cassuto/CORDIC-all-in-one-verilog/blob/master/CORDIC-rotate-fixed-point.c
反旋转算法：
https://github.com/cassuto/CORDIC-all-in-one-verilog/blob/master/CORDIC-anti-rotate-fixed-point.c
4、FPGA实现
这一小节主要解决的问题是硬件实现的相关细节，例如象限转换，流水化处理等。采用流水线的目的在于提升时钟频率，例如在DDS（直接数字频率合成）应用中，CORDIC算法能够代替采样表生成波形数据。为了实现较高的合成频率，流水化是有必要的。
关于Verilog实现有几点须要说明的地方，首先利用verilog标准中规定的generate语句，能够实现任意深度流水线的综合；
其次在此以前对算法的全部讨论都只涉及第Ⅰ象限，如今加入象限的处理。电路直接取相位输入高2位判断象限，并按函数值在4个象限的符号关系对输出进行处理。
完整的Verilog模块以下：该模块经过端口phase_i输入相位。经过sin_o，cos_o输出三角函数值，经过err_o输出迭代引发的相位偏差。从流水线为空开始，需等待PIPE_DEPTH+2个时钟周期才能得出结果；而流水线填满后，对于后续相位输入，模块均可以在2个时钟周期内更新输出。
1 module cordic_dds # ( 2 parameter DW = 16, /* Data width */ 3 parameter PIPE_DEPTH = 14, /* Pipeline depth */ 4 parameter limP = 16'h4dba /* P = 0.607253 * 2^15 */ 5 ) 6 (/*AUTOARG*/ 7 // Outputs 8 sin_o, cos_o, err_o, 9 // Inputs 10 clk, phase_i 11 ); 12 13 input clk; 14 input [DW-1:0] phase_i; /* Phase */ 15 output [DW:0] sin_o, cos_o; /* Function value output */ 16 output [DW:0] err_o; /* Phase Error output */ 17 18 reg [DW:0] cos_r=0, sin_o_r=0; 19 reg [DW:0] x[PIPE_DEPTH:0]; 20 reg [DW:0] y[PIPE_DEPTH:0]; 21 reg [DW:0] z[PIPE_DEPTH:0]; 22 23 reg [DW:0] atan_rom[PIPE_DEPTH:0]; 24 25 reg [1:0] quadrant [PIPE_DEPTH:0]; 26 27 integer i; 28 initial begin 29 for(i=0; i<=PIPE_DEPTH; i=i+1) begin 30 x[i] = 0; y[i] = 0; z[i] = 0; 31 quadrant[i] = 2'b0; 32 end 33 end 34 35 initial begin 36 atan_rom[0] <= 8189; 37 atan_rom[1] <= 4834; 38 atan_rom[2] <= 2554; 39 atan_rom[3] <= 1296; 40 atan_rom[4] <= 650; 41 atan_rom[5] <= 325; 42 atan_rom[6] <= 162; 43 atan_rom[7] <= 81; 44 atan_rom[8] <= 40; 45 atan_rom[9] <= 20; 46 atan_rom[10] <= 10; 47 atan_rom[11] <= 5; 48 atan_rom[12] <= 2; 49 atan_rom[13] <= 1; 50 end 51 52 53 // ================= // 54 // Pipeline stages // 55 // ================= // 56 always @ (posedge clk) begin // stage 0 57 x[0] <= {1'b0, limP}; 58 y[0] <= 0; 59 z[0] <= {3'b0, phase_i[DW-1-2:0]}; // control the phase_i to the range[0-Pi/2] 60 end 61 62 always @ (posedge clk) begin // stage 1 63 x[1] <= x[0] - y[0]; 64 y[1] <= x[0] + y[0]; 65 z[1] <= z[0] - atan_rom[0]; // reversal 45deg 66 end 67 68 generate 69 genvar k; 70 for(k=1; k<PIPE_DEPTH; k=k+1) begin 71 always @ (posedge clk) begin 72 if (z[k][DW]) begin /* the diff is negative on clockwise */ 73 x[k+1] <= x[k] + {{k{y[k][DW]}},y[k][DW:k]}; /* >> k */ 74 y[k+1] <= y[k] - {{k{x[k][DW]}},x[k][DW:k]}; /* >> k */ 75 z[k+1] <= z[k] + atan_rom[k]; 76 end else begin 77 x[k+1] <= x[k] - {{k{y[k][DW]}},y[k][DW:k]}; 78 y[k+1] <= y[k] + {{k{x[k][DW]}},x[k][DW:k]}; 79 z[k+1] <= z[k] - atan_rom[k]; 80 end 81 end 82 end 83 endgenerate 84 85 // ================= // 86 // Count quadrant // 87 // ================= // 88 always @ (posedge clk) begin 89 quadrant[0] <= phase_i[DW-1:DW-2]; 90 end 91 generate 92 genvar j; 93 for(j=0; j<PIPE_DEPTH; j=j+1) begin 94 always @ (posedge clk) begin 95 quadrant[j+1] <= quadrant[j]; 96 end 97 end 98 endgenerate 99 100 // ================= // 101 // Adjust quadrant // 102 // ================= // 103 always @ (posedge clk) 104 case(quadrant[PIPE_DEPTH]) 105 2'b00: begin 106 cos_r <= x[PIPE_DEPTH]; /* cos */ 107 sin_o_r <= y[PIPE_DEPTH]; /* sin */ 108 end 109 2'b01: begin 110 cos_r <= ~(y[PIPE_DEPTH]) + 1'b1; /* -sin */ 111 sin_o_r <= x[PIPE_DEPTH]; /* cos */ 112 end 113 2'b10: begin 114 cos_r <= ~(x[PIPE_DEPTH]) + 1'b1; /* -cos */ 115 sin_o_r <= ~(y[PIPE_DEPTH]) + 1'b1; /* -sin */ 116 end 117 2'b11: begin 118 cos_r <= y[PIPE_DEPTH]; /* sin */ 119 sin_o_r <= ~(x[PIPE_DEPTH]) + 1'b1; /* -cos */ 120 end 121 endcase 122 123 assign cos_o = cos_r; 124 assign sin_o = sin_o_r; 125 assign err_o = z[PIPE_DEPTH]; 126 127 endmodule
经过仿真得出波形以下：
本文仅对CORDIC算法进行了一些浅显的分析，基本覆盖到了算法的核心思路和简单应用，但从工业应用的角度出发，还有不少值得讨论的问题未在文章中分析。因为笔者时间有限，不能将本文作得很是全面，疏漏之处有待往后逐步完善，还望各位读者海涵。