最近在看cordic算法，由于还不会使用matlab，真是痛苦，一系列的笔算才大概明白了这个算法是怎么回事。于是尝试用verilog来实现。用verilog实现之前先参考软件的程序，于是先看了此博文http://blog.csdn.net/liyuanbhu/article/details/8458769 也不截图了，因为怕图形被其他博客网站检测到后屏蔽图片，造成此博文无法正常阅读。
阅读此博文，需要先阅读上面这个博文的内容。
这是此博文中的C代码。避免浮点运算，所以angle数组里面的角度值都扩大了256倍。此程序中计算点（10,20）与X轴之间的夹角。最终运行结果是16238。而16238/256=63.43
1 #include <stdio.h> 2 #include <stdlib.h> 3 4 //double my_atan2(double x, double y); 5 6 7 double my_atan2 (int x, int y) 8 { 9 const int angle[] = {11520, 6801, 3593, 1824, 916, 458, 229, 115, 57, 29, 14, 7, 4, 2, 1}; 10 11 int i = 0; 12 int x_new, y_new; 13 int angleSum = 0; 14 15 x *= 1024;// ½« X Y ·Å´óһЩ£¬½á¹û»á¸ü׼ȷ 16 y *= 1024; 17 18 printf ("org_x = %d, org_y=%d ",x,y); 19 20 for(i = 0; i < 15; i++) 21 { 22 if(y > 0) 23 { 24 x_new = x + (y >> i); 25 y_new = y - (x >> i); 26 x = x_new; 27 y = y_new; 28 angleSum += angle[i]; 29 } 30 else 31 { 32 x_new = x - (y >> i); 33 y_new = y + (x >> i); 34 x = x_new; 35 y = y_new; 36 angleSum -= angle[i]; 37 } 38 printf("Debug: i = %d x=%d, y =%d angleSum = %d, angle = %d ", i,x,y ,angleSum, angle[i]); 39 } 40 return angleSum; 41 } 42 43 void main() 44 { 45 double z=0 ; 46 z = my_atan2(10.0, 20.0); 47 printf(" z = %lf ", z); 48 49 }
既然有了C 就很好进行verilog设计。
先谈架构。
1，先将角度数据放到一个rom中进行存储
2，取一个数，运算一个数。直到循环结束。
这也就是所谓的cordic算法的向量模式。用来求角度。先看顶层
1 // 2 // 3 // 4 // 5 6 module cordic_rotation_top ( 7 clock , 8 rst_n , 9 x_crd, 10 y_crd, 11 ena , 12 deg_sum 13 ); 14 input clock ; 15 input rst_n ; 16 input [15:0] x_crd ; 17 input [15:0] y_crd ; 18 input ena ; 19 output [15:0] deg_sum ; 20 21 22 wire [4:0] deg_addr ; 23 wire [15:0] deg_data ; 24 25 alt_ip_rom_cordic u_rom ( 26 .address (deg_addr), 27 .clock (clock), 28 .q (deg_data) 29 ); 30 31 32 33 cordic_rotation u_cord ( 34 .clk (clock), 35 .rst_n (rst_n), 36 .ena (ena), 37 .x_crd (x_crd), 38 .y_crd (y_crd), 39 .deg_data (deg_data), 40 .deg_addr (deg_addr), 41 .deg_sum (deg_sum) 42 ); 43 44 endmodule
rom的初始化文件为
再看运算单元。
1 module cordic_rotation ( 2 clk , 3 rst_n , 4 ena , 5 x_crd, 6 y_crd, 7 deg_data, 8 deg_addr, 9 deg_sum 10 ); 11 12 input clk ; 13 input rst_n ; 14 input ena ; 15 input [15:0] x_crd ; //x coordinate 16 input [15:0] y_crd ; //y coordinate 17 input [15:0] deg_data ; 18 output [4:0] deg_addr ; 19 output reg[15:0] deg_sum ; 20 21 // ------ rotation count 0 - 14 ------- 22 reg [4:0] rot_cnt ; 23 reg [4:0] rot_cnt_r ; 24 wire opr_en ; 25 assign opr_en = ((rot_cnt_r<5'd15)&(ena)) ; 26 27 always @ (posedge clk or negedge rst_n) 28 if (!rst_n) begin 29 rot_cnt <= 5'd0 ; 30 rot_cnt_r<= 5'd0; 31 end 32 else if (opr_en) begin 33 rot_cnt <= rot_cnt + 5'd1 ; 34 rot_cnt_r<= rot_cnt ; 35 end 36 else if (!ena) begin 37 rot_cnt <= 5'd0 ; 38 rot_cnt_r<= rot_cnt ; 39 end 40 assign deg_addr = rot_cnt ; 41 42 //--------------------------------------- 43 reg cal_cnt ; 44 reg signed [16:0] x_d, y_d ; 45 always @ (posedge clk or negedge rst_n) 46 if (!rst_n) begin 47 x_d <= 17'd0 ; 48 y_d <= 17'd0 ; 49 deg_sum <= 16'd0 ; 50 cal_cnt <= 1'd0 ; 51 end 52 else if (opr_en) begin 53 case (cal_cnt) 54 1'd0 : begin 55 x_d <= {1'd0,x_crd}; 56 y_d <= {1'd0,y_crd}; 57 deg_sum <= 16'd0 ; 58 cal_cnt <= 1'd1 ; 59 end 60 1'd1 : begin 61 if ((y_d[16])|(y_d==16'd0)) begin 62 x_d <= x_d - (y_d >>> rot_cnt_r); 63 y_d <= y_d + (x_d >>> rot_cnt_r) ; 64 deg_sum <= deg_sum - deg_data ; 65 end 66 else begin 67 x_d <= x_d + (y_d >>> rot_cnt_r); 68 y_d <= y_d - (x_d >>> rot_cnt_r) ; 69 deg_sum <= deg_sum + deg_data ; 70 end 71 end 72 endcase 73 end 74 else begin 75 cal_cnt <= 1'd0 ; 76 end 77 78 endmodule
rot_cnt作为rom的地址。但是由于rom返回度数值需要一个周期，所以下面的运算单元需要等待一个周期才可以进行运算，于是有了rot_cnt_r这个参数。
于是问题又来了，上面的是向量模式。还有一个旋转模式。于是开始捣鼓。
稍微将上述的C代码进行修改，计算30度的cos30，以及sin30
1 #include <stdio.h> 2 #include <stdlib.h> 3 4 //double my_atan2(double x, double y); 5 6 7 double my_atan2 (int x, int y) 8 { 9 const int angle[] = {11520, 6801, 3593, 1824, 916, 458, 229, 115, 57, 29, 14, 7, 4, 2, 1}; 10 11 int i = 0; 12 int x_new, y_new; 13 int angleSum = 30*256; 14 15 x *= 1024;// ½« X Y ·Å´óһЩ£¬½á¹û»á¸ü׼ȷ 16 y *= 1024; 17 18 printf ("org_x = %d, org_y=%d ",x,y); 19 20 for(i = 0; i < 15; i++) 21 { 22 if(angleSum <= 0) 23 { 24 x_new = x + (y >> i); 25 y_new = y - (x >> i); 26 x = x_new; 27 y = y_new; 28 angleSum += angle[i]; 29 } 30 else 31 { 32 x_new = x - (y >> i); 33 y_new = y + (x >> i); 34 x = x_new; 35 y = y_new; 36 angleSum -= angle[i]; 37 } 38 printf("Debug: i = %d x=%d, y =%d angleSum = %d, angle = %d ", i,x,y ,angleSum, angle[i]); 39 } 40 return angleSum; 41 } 42 43 void main() 44 { 45 double z=0 ; 46 z = my_atan2(1.0, 0.0); 47 printf(" z = %lf ", z); 48 49 }
结果是这样子的
先看已知条件：
1，cordic算法中的Kn极限值=1.6476 。 1/Kn=0.6073.
2，上述软件设置y=0. x=1.并且坐标值扩大了1024倍。从公式上，cosa= x sina= y 角度a=30度
3，运算结果x=1461. y=844
好，开始运算 cosa=1461/(1024*1.6476)=0.8659。 cos30=0.8660.
sina=844/(1024*1.6476) = 0.500 。 sin30 = 0.5
精度由loop次数相关
好，既然验证了这个软件程序是可以的。那么就将它转换成Verilog的程序
1 //cordic 算法的旋转模式 2 3 4 module cordic_rotation ( 5 clk , 6 rst_n , 7 ena , 8 x_crd, 9 y_crd, 10 deg_data, 11 deg_addr, 12 deg_sum 13 ); 14 15 input clk ; 16 input rst_n ; 17 input ena ; 18 input [15:0] x_crd ; //x coordinate 19 input [15:0] y_crd ; //y coordinate 20 input [15:0] deg_data ; 21 output [4:0] deg_addr ; 22 output [15:0] deg_sum ; 23 24 25 parameter DEGREE = 30 ; 26 parameter DEG_PARA = (DEGREE*256) ; 27 // ------ rotation count 0 - 14 ------- 28 reg [4:0] rot_cnt ; 29 reg [4:0] rot_cnt_r ; 30 wire opr_en ; 31 assign opr_en = ((rot_cnt_r<5'd15)&(ena)) ; 32 33 always @ (posedge clk or negedge rst_n) 34 if (!rst_n) begin 35 rot_cnt <= 5'd0 ; 36 rot_cnt_r<= 5'd0; 37 end 38 else if (opr_en) begin 39 rot_cnt <= rot_cnt + 5'd1 ; 40 rot_cnt_r<= rot_cnt ; 41 end 42 else if (!ena) begin 43 rot_cnt <= 5'd0 ; 44 rot_cnt_r<= rot_cnt ; 45 end 46 assign deg_addr = rot_cnt ; 47 48 //--------------------------------------- 49 reg cal_cnt ; 50 reg signed [15:0] deg_sum_r ; 51 reg signed [16:0] x_d, y_d ; 52 always @ (posedge clk or negedge rst_n) 53 if (!rst_n) begin 54 x_d <= 17'd0 ; 55 y_d <= 17'd0 ; 56 deg_sum_r <= 16'd0 ; 57 cal_cnt <= 1'd0 ; 58 end 59 else if (opr_en) begin 60 case (cal_cnt) 61 1'd0 : begin 62 x_d <= {1'd0,x_crd}; 63 y_d <= {1'd0,y_crd}; 64 deg_sum_r <= DEG_PARA ; 65 cal_cnt <= 1'd1 ; 66 end 67 1'd1 : begin 68 if ((deg_sum_r[15])|(deg_sum_r==16'd0)) begin //<=0 69 x_d <= x_d + (y_d >>> rot_cnt_r); 70 y_d <= y_d - (x_d >>> rot_cnt_r) ; 71 deg_sum_r <= deg_sum_r + deg_data ; 72 end 73 else begin 74 x_d <= x_d - (y_d >>> rot_cnt_r); 75 y_d <= y_d + (x_d >>> rot_cnt_r) ; 76 deg_sum_r <= deg_sum_r - deg_data ; 77 end 78 end 79 endcase 80 end 81 else begin 82 cal_cnt <= 1'd0 ; 83 end 84 85 assign deg_sum = deg_sum_r ; 86 87 endmodule
这个程序也是在上一个程序的基础上修改的。所以尽量少改动。而是修改了tb。所以需要看看TB是怎么弄的
1 /// 2 // 3 // 4 // 5 `timescale 1ns/100ps 6 7 8 module cordic_rotation_top_tb ; 9 reg clock ,rst_n ; 10 reg [15:0] x_crd ,y_crd ; 11 reg ena ; 12 13 wire [15:0] deg_sum ; 14 15 16 cordic_rotation_top u_top ( 17 .clock (clock), 18 .rst_n (rst_n), 19 .x_crd (x_crd), 20 .y_crd (y_crd), 21 .ena (ena), 22 .deg_sum (deg_sum) 23 ); 24 25 always #10 clock = ~clock ; 26 initial begin 27 clock = 0 ; rst_n =0 ; 28 x_crd = 1*1024 ; y_crd=0; 29 ena = 1 ; 30 #40 rst_n = 1 ; 31 #350 32 /* 33 ena = 0 ; 34 #40 35 x_crd = 10*1024 ; y_crd=10*1024; 36 ena = 1 ; 37 #350 38 39 ena = 0 ; 40 #40 41 x_crd = 20*1024 ; y_crd=10*1024; 42 ena = 1 ; 43 #350 */ 44 $stop ; 45 end 46 47 endmodule
将第28行替换为 x_crd = 10*1024 ; y_crd=20*1024; 就变成了上一个工程的测试代码。
此程序输出结果为：
和软件程序输出结果相同。
有没有发现问题，在计算这些东西的时候，很多网络博文说精度和loop次数相关。但是这个扩大256倍运算到第11次就已经定住了。后面的输出结果白循环了。
欢迎加入： FPGA广东交流群：162664354
FPGA开发者联盟： 485678884
微信公众号：FPGA攻城狮之家
