最近在看cordic算法，由于还不会使用matlab，真是痛苦，一系列的笔算才大概明白了这个算法是怎么回事。于是尝试用verilog来实现。用verilog实现之前先参考软件的程序，于是先看了此博文也不截图了，因为怕图形被其他博客网站检测到后屏蔽图片，造成此博文无法正常阅读。
阅读此博文，需要先阅读上面这个博文的内容。
这是此博文中的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\n",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\n", 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("\n z = %lf \n", 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 endmodulerom的初始化文件为
再看运算单元。
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 endmodulerot_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\n",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\n", 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("\n z = %lf \n", 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次就已经定住了。后面的输出结果白循环了。
本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。
