verilog實現基於Cordic演算法的雙曲函式計算

       Cordic演算法可以用FPGA硬體來實現三角函式，向量旋轉，指數函式以及三角函式等數值計算，它是一種從一般的向量旋轉方程中推導得出。採用用不斷的旋轉求出對應的正弦餘弦值，是一種近似求解法。旋轉的角度很講求，比如求取正餘弦函式值時每次旋轉的角度必須使得 正切值近似等於 1/(2^N)。旋轉的目的是讓Y軸趨近與0。把每次旋轉的角度累加，即得到旋轉的角度和即為正切值。如圖1所示為Cordic基本原理示意圖，圖2位求解三角函式示意圖。

圖1 Cordic基本原理

圖2 Cordic求解三角函式

     顯然Cordic演算法只需要簡單的硬體電路就可以實現。根據旋轉方向判定，可以選擇模式或向量模式進行工作，並通過設定不同的初始條件可以實現特殊函式的數學計算，從而為該演算法開創更廣泛的應用空間。

在雙曲座標下，Cordic演算法的迭代方程為：

• 與其他Cordic計算不同，由於tanh(-1)(2^0)為無窮大，所以迭代序列須從(i=1)開始，可以將各i處的值計算好並儲存在ROM中實現一個小的查詢表（LUT）。並從收斂考慮序列i的取值從第4項開始每隔3i+1項須重複一次。n次迭代輸出方程為：
• 
  
• 可以通過選擇適當的初始值及多種操作模式組合完成tanh,exp的計算。下面是初步簡單實現的Verilog程式碼,程式碼中同時有用到ISE自帶的Cordic演算法的IP核作了模擬對比：
• [cpp] view plain copy

  1. module tanh_test(  
  2.         clk,  
  3.     rst_n,  
  4.     z0,//input  
  5.     en,//input  
  6.     cosh_out,  
  7.     sinh_out,  
  8.     e_out,//output  
  9.     busy_end//output  
  10.       
  11.     );  
  12.   
  13. input clk,rst_n;  
  14. input [31:0] z0;//16  
  15. input en;  
  16. output [31:0] sinh_out,cosh_out;//16  
  17. output [31:0] e_out;//16  
  18. output busy_end;  
  19.   
  20. reg en_buff;  
  21.   
  22. wire sub_busy;  
  23.   
  24. reg [31:0] cnt,h_cnt;  
  25. wire [31:0] sub_result;  
  26. reg busy_e;//  
  27.   
  28.   
  29. always @(posedge clk or negedge rst_n)  
  30. begin  
  31.          if(!rst_n)  
  32.              begin  
  33.              en_buff <= 0;  
  34.              busy_e <= 1'b1;//  
  35.              cnt <= 0;  
  36.              h_cnt <= 32'b0;  
  37.              end  
  38.          else begin  
  39.                   en_buff <= en;  
  40.               cnt <= cnt + 1;  
  41.                  end           
  42. end  
  43.   
  44. reg[31:0] angel [0:12];  
  45.   
  46. always @(posedge clk or negedge rst_n)  
  47. begin  
  48.       if(!rst_n)  
  49.       begin  
  50.        angel[0]<=32'b00000000000000001000110010011111;// 0.549306,1/2,tanh  
  51.            angel[1]<=32'b00000000000000000100000101100010;//0.255413,1/4  
  52.            angel[2]<=32'b00000000000000000010000000101011;//0.125657,1/8  
  53.            angel[3]<=32'b00000000000000000001000000000101;//0.062582,1/16  
  54.            angel[4]<=32'b00000000000000000001000000000101;//重複迭代一次  
  55.            angel[5]<=32'b00000000000000000000100000000000;//0.031260，1/32  
  56.            angel[6]<=32'b00000000000000000000010000000000;//0.015626，1/64  
  57.            angel[7]<=32'b00000000000000000000001000000000;//0.007813，1/128  
  58.        angel[8]<=32'b00000000000000000000000100000000;//0.003906250，1/256  
  59.            angel[9]<=32'b00000000000000000000000010000000;//0.001953125，1/128  
  60.        angel[10]<=32'b00000000000000000000000001000000;//0.0009765625，1/128  
  61.        angel[11]<=32'b00000000000000000000000000100000;//0.00048828125，1/128  
  62.        angel[12]<=32'b00000000000000000000000000010000;//0.000244140625，1/128  
  63.        //angel[8]<=17'b00000000000110011;//0.2,1/256,17'b00000000000110011  
  64.            //angel[9]<=17'b00000000000011001;//0.1,1/512,angel[9]<=17'b00000000000011001  
  65.      end  
  66. end  
  67. reg[31:0] reg_z[0:13];//1符號，15整數，16小數  
  68. reg[31:0] reg_x[0:13];  
  69. reg[31:0] reg_y[0:13];  
  70. reg[4:0] i;  
  71.   
  72. always @(posedge clk or negedge rst_n)  
  73. begin  
  74.          if(!rst_n)  
  75.          begin  
  76.              h_cnt <= 0;  
  77.              end  
  78.          else if(en && !en_buff)//&& (busy_e == 1'b1)  
  79.          begin  
  80.               reg_x[0] <= 32'b00000000000000010011010100011000;//（1/0.8282=1.2074）  
  81.               reg_y[0] <= 0;  
  82.               reg_z[0] <= z0;//初始值為v=1  
  83.              end  
  84.          else begin  
  85.              h_cnt <= h_cnt+1'b1;/////////////////////////  
  86.              for(i = 1;i <= 13;i = i+1'b1)  
  87.              begin  
  88.                    if(reg_z[i-1][31])  
  89.                    begin  
  90.                       reg_x[i] <= reg_x[i-1]-(reg_y[i-1]>>i);//z<0,d=-1,否則d=1  
  91.                       reg_y[i] <= reg_y[i-1]-(reg_x[i-1]>>i);  
  92.                       reg_z[i] <= reg_z[i-1]+angel[i-1];  
  93.                        end  
  94.                    else begin  
  95.                       reg_x[i] <= reg_x[i-1]+(reg_y[i-1]>>i);  
  96.                       reg_y[i] <= reg_y[i-1]+(reg_x[i-1]>>i);  
  97.                       reg_z[i] <= reg_z[i-1]-angel[i-1];  
  98.                    end  
  99.              end  
  100.                     
  101.              if(h_cnt >= 20)  
  102.              begin  
  103.                   busy_e <= 1'b0;  
  104.              end  
  105.                     
  106.          end  
  107. end  
  108.   
  109. assign sinh_out = reg_y[11];//<<4  
  110. assign cosh_out = reg_x[11];  
  111. assign e_out = sub_result;  
  112. assign busy_end =(sub_busy | busy_e);  
  113.   
  114. reg [15:0] in_phase;  
  115. wire [15:0] out_cosh;  
  116. wire [15:0] out_sinh;  
  117. initial  
  118. begin  
  119.      in_phase <= 16'b0010000000000000;//（1，第16位符號位，第15，14位整數位，13位小數位）  
  120. end  
  121.   
  122.   cordic_e your_instance_name (  
  123.            .phase_in(in_phase), // input [15 : 0] phase_in  
  124.            .x_out(out_cosh), // output [15 : 0] x_out  
  125.            .y_out(out_sinh), // output [15 : 0] y_out  
  126.            .clk(clk) // input clk  
  127.   );  
  128.   e_sub cosh_sub_sinh(  
  129.        
  130.         .clk(clk),  
  131.         .rst_n(rst_n),  
  132.         .sub_a(cosh_out),  
  133.         .sub_b(sinh_out),  
  134.         .sub_result(sub_result),  
  135.         .busy_e(busy_e),  
  136.         .sub_busy(sub_busy)  
  137.            
  138.    );  
  139.   
  140.   
  141. endmodule  

        模擬結果圖，計算角度為1時的值：

（32'h00018aae=1.5417,sinh(1)真實值1.5431;32'h00012bfd=1.1718,cosh(1)=1.1752;32'h00005eb1=0.3699,e^(-1)=0.3679。使用IP核計算結果：16'h62d1=1.5440,16'h4b4b=1.1765）