最近在看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\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 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\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次就已經定住了。后面的輸出結果白循環了。
歡迎加入： FPGA廣東交流群：162664354
FPGA開發者聯盟： 485678884
微信公眾號：FPGA攻城獅之家
