計算機該如何求解三角函數？或許你的第一印象是采用泰勒展開，或者采用多項式進行逼近。對于前者，來回的迭代計算開銷成本很大；對于后者，多項式式逼近在較窄的范圍內比較接近，超過一定范圍后，就變得十分不理想了。例如x–>0時，x~sin(x)

今天，我們將要介紹三角函數的另一種求法：CORDIC算法

原理

CORDIC的算法的核心就是通過迭代，利用三角函數的物理性質，不斷累積旋轉角度，從而得到所求角度的精確近似。

我們假設圓為單位圓，范圍且在第一象限，如下圖：

我們假設點（x1,y1）與X軸正半軸夾角為α,那么
$$\{\begin{array}{l}{y}_2=sin(\alpha +\theta )\\ x_2=cos(\alpha +\theta )\end{array}$$
三角函數展開有
$$\{\begin{array}{l}{y}_2=sin(\alpha )cos(\theta )+cos(\alpha )sin(\theta )\\ x_2=cos(\alpha )cos(\theta )?sin(\alpha )sin(\theta )\end{array}$$

將
$$\{\begin{array}{l}{y}_1=sin(\alpha )\\ x_1=cos(\alpha )\end{array}$$
帶入上式，有
$$\{\begin{array}{l}{y}_2=sin(\alpha )cos(\theta )+cos(\alpha )sin(\theta )=y_{1}cos(\theta )+x_{1}sin(\theta )=cos(\theta )(y_1+x_{1}tan(\theta ))\\ x_2=cos(\alpha )cos(\theta )?sin(\alpha )sin(\theta )=x_{1}cos(\theta )?y_{1}sin(\theta )=cos(\theta )(x_1?y_{1}tan(\theta ))\end{array}$$

默認初始值為（1，0），記為
$$v_0=\left[\begin{array}{c}1\\ 0\end{array}\right]$$
以上的等式可以表示為旋轉矩陣的形式
$$\left[\begin{array}{c}{x}_2\\ y_2\end{array}\right]=cos(\alpha )\left[\begin{array}{cc}1& ?tan(\alpha )\\ tan(\alpha )& 1\end{array}\right]\left[\begin{array}{c}{x}_1\\ y_1\end{array}\right]$$
如果將令角度α，滿足tan(α)=2^-i, 那么就將tan和乘法運算就變成乘2的負次冪。對應在計算機中，就是移位計算。因而復雜的計算，就變成了簡單的加減和移位運算。

所以我們有
$$\left[\begin{array}{c}{x}_n\\ y_n\end{array}\right]=cos({\alpha }_n)\left[\begin{array}{cc}1& ?2^{?n}\\ 2^{?n}& 1\end{array}\right]\left[\begin{array}{c}{x}_{n?1}\\ y_{n?1}\end{array}\right]=cos({\alpha }_n)cos({\alpha }_{n?1})..cos({\alpha }_0)\left[\begin{array}{cc}1& ?2^{?n}\\ 2^{?n}& 1\end{array}\right]\left[\begin{array}{cc}1& ?2^{?n+1}\\ 2^{?n+1}& 1\end{array}\right]..\left[\begin{array}{cc}1& ?2^{?0}\\ 2^{?0}& 1\end{array}\right]\left[\begin{array}{c}1\\ 0\end{array}\right]$$

處理細節

1. 縮放因子

由前面推導，我們可以得到：
$$\left[\begin{array}{c}{x}_n\\ y_n\end{array}\right]=cos({\alpha }_n)\left[\begin{array}{cc}1& ?2^{?n}\\ 2^{?n}& 1\end{array}\right]\left[\begin{array}{c}{x}_{n?1}\\ y_{n?1}\end{array}\right]=cos({\alpha }_n)cos({\alpha }_{n?1})..cos({\alpha }_0)\left[\begin{array}{cc}1& ?2^{?n}\\ 2^{?n}& 1\end{array}\right]\left[\begin{array}{cc}1& ?2^{?n+1}\\ 2^{?n+1}& 1\end{array}\right]..\left[\begin{array}{cc}1& ?2^{?0}\\ 2^{?0}& 1\end{array}\right]\left[\begin{array}{c}1\\ 0\end{array}\right]$$

我們可以提前將所有的cos(α)的乘積計算出來，做為一個常量，省去乘法運算，記為K
$$K=cos({\alpha }_n)cos({\alpha }_{n?1})cos({\alpha }_{n?2})...cos({\alpha }_0)=0.60725$$
2. 旋轉方向

通常來說CORDIC算法會引入d_n ，來判斷旋轉方向。當前角度大于該次迭代的角度，dn_{為正，逆時鐘旋轉，反之為負，順時針旋轉。之所以會采用雙旋轉，是因為其通常比單向旋轉的收斂性更好，結果更精確。}

因而我們迭代可以寫為
$$?\begin{array}{l}{y}_{n+1}=y_n+d_n?x_n?2^{?n}\\ x_{n+1}=x_n?d?y_n?2^{?n}\\ angl{e}_{n+1}=angl{e}_n?d?tableofangles[n]\end{array}$$
table_of_angles 存儲的是θ值， θ_n=arctan(2^-n);

對應的表格如下：

n	2^(-n)	arctan(2^(-n))
0	1	0.785398163
1	0.5	0.463647609
2	0.25	0.244978663
3	0.125	0.124354995
4	0.0625	0.06241881
5	0.0315	0.031239833
6	0.015625	0.015623729
7	0.0078125	0.007812341
8	0.00390625	0.00390623
9	0.001953125	0.001953123
10	0.000976563	0.000976562
11	0.000488281	0.000488281
12	0.000244141	0.000244141
13	0.00012207	0.00012207
14	6.10352E-05	6.10352E-05
15	3.05176E-05	3.05176E-05

下面我會手把手帶領大家從軟件建模到硬件實現CORDIC算法，規定輸入和輸出都是無符號17位數，1位整數位，16位小數位。

Python 代碼

測試代碼

#初始化部分，定義參數
import math
from math import floorNUM_ITER = 16
Frac_Bits=16
Data_Scale=2**Frac_Bits
Angles_Table=[]#創建對應的對應的角度表
def create_angel_table():   for i in range(NUM_ITER):angles=math.atan(2**(-i))angles=floor(angles*Data_Scale+0.5)/Data_Scale#print(angles)#print(angles)#angles=angles*(1<<Frac_Bits)+0.5#angles=floor(angles)#print(angles)#print(hex(angles))Angles_Table.append(angles)#計算出縮放因子
def compute_k():k=1.0for i in range(NUM_ITER):angles=math.atan(2**(-i))k=k*math.cos(angles)#print(K)#print(hex(floor(K*(1<<Frac_Bits)+0.5)))return  floor(k*Data_Scale+0.5)/Data_Scale  # cordic 算法迭代
def cordic(theta,k):x=ky=0angle_temp= floor(math.radians(theta)*Data_Scale+0.5)/Data_Scale  for i in range(NUM_ITER):if(angle_temp>=0):x_next=x-y*2**(-i)y_next=y+x*2**(-i)angle_temp-=Angles_Table[i]else:   x_next=x+y*2**(-i)y_next=y-x*2**(-i)angle_temp+=Angles_Table[i]x=x_nexty=y_nextreturn x,y#cordic 算法算出的結果，與真實結果進行比較
def compare(ground_truth, test):for i in range(len(ground_truth)): # 如果誤差超過 3*2^（-16）次，那么退出比較if( abs(ground_truth[i]-test[i])>3):print("Error! Loss of accuracy! ground_truth: %f, test: %f", ground_truth[i], test[i])return Falsereturn True
#得到cordic算法結果，經行比較
def main():create_angel_table()k=compute_k()cos_truth=[]sin_truth=[]cos_test=[]sin_test=[]for i in range(90):cos_truth.append(floor(math.cos(i*math.pi/180)*Data_Scale+0.5))sin_truth.append(floor(math.sin(i*math.pi/180)*Data_Scale+0.5))cos_temp,sin_temp=cordic(i,k)cos_test.append(floor(cos_temp*Data_Scale+0.5))sin_test.append(floor(sin_temp*Data_Scale+0.5))if (compare(cos_truth,cos_test) and compare(sin_truth,sin_test)):print("Test Pass")else:print("Test Fail")if __name__ == "__main__":main()

比較結果

由此可知，CORDIC算法精度很高

Verilog 代碼

模塊代碼

module Cordic_Sin(input wire [16:0] theta,   // 輸入角度（Q1.16格式，范圍0 ~ π/2）output wire [16:0] sin_out, // 輸出sin值（Q1.16格式）output wire [16:0] cos_out
);// 預計算參數（Q1.16格式）
localparam signed [16:0] K =17'sh09B75;  // 1/1.64676補償因子;  17'h1A592;         //Q1.15
reg signed [16:0] angles [0:16];    //arctan(2^-i)
integer iter;
initial beginangles[0]  = 17'h0C910;angles[1]  = 17'h076B2;angles[2]  = 17'h03EB7;angles[3]  = 17'h01FD6;// i=0~3angles[4]  = 17'h00FFB;angles[5]  = 17'h007FF;angles[6]  = 17'h00400;angles[7]  = 17'h00200;// i=4~7angles[8]  = 17'h00100;angles[9]  = 17'h00080;angles[10] = 17'h00040;angles[11] = 17'h00020;// i=8~11angles[12] = 17'h00010;angles[13] = 17'h00008;angles[14] = 17'h00004;angles[15] = 17'h00002;// i=12~15angles[16] = 17'h00001;
endreg signed [32:0]x,y; //初始化 x=K; y=0
reg signed [32:0]x_next,y_next;
reg signed [17:0] angle; // 初始化角度等于輸出角度
integer i;always@(*)begin//初始化x={K,16'b0};y=33'h0;angle={1'b0,theta};//迭代計算for(i=0;i<16;i=i+1)beginif(!angle[17])begin   ////正向旋轉x_next=x-(y>>>i);  //算術移位y_next=y+(x>>>i);angle=angle-{1'b0,angles[i]};end else begin//負向旋轉x_next=x+(y>>>i);y_next=y-(x>>>i);angle=angle+{1'b0,angles[i]};endx=x_next;y=y_next;end
endassign sin_out = y[32:16];
assign cos_out = x[32:16];endmodule

Test Bench

`timescale 1ns / 1psmodule Cordic_Sin_Test();reg  [16:0] theta; wire [16:0] sin_out;wire [16:0] cos_out;initial begintheta=17'h0;#10;//15 theta=17'h4305;#10;//30theta=17'h860A;#10;//45theta=17'hC90F;#10;//60theta=17'h10C15;#10;//75theta=17'h14F1A;#10;//90theta=17'h1921F;#10;endCordic_Sin uut(.theta(theta),   // 輸入角度（Q1.16格式，范圍0 ~ π/2）.sin_out(sin_out),// 輸出sin值（Q1.16格式）.cos_out(cos_out)
);endmodule

仿真結果

以上的數據，輸入數據需要將其轉換成弧度值，然后轉換成S0I1F16定點格式，sin,cos準確值也是一樣

輸入角度	sin準確值	cos準確值	sin計算值	cos計算值
0°	0	65536	154	65536
15°	16962	63303	16962	63302
30°	32768	56756	32768	56755
45°	46341	46341	46340	46341
60°	56756	32768	56755	32768
75°	63303	16962	63302	16962
90°	65536	0	65536	154

除了個別點外的絕對誤差比較大外，其余的計算精度相當高，誤差很小