计算机该怎样求解三角函数？或许你的第一印象是采用泰勒睁开，大概采用多项式进行逼近。对于前者，往返的迭代计算开销本钱很大；对于后者，多项式式逼近在较窄的范围內比力接近，超过肯定范围后，就变得十分不抱负了。例如x–>0时，x~sin(x)
今天，我们将要先容三角函数的另一种求法：CORDIC算法
原理
CORDIC的算法的核心就是通过迭代，利用三角函数的物理性质，不绝累积旋转角度，从而得到所求角度的精确近似。
我们假设圆为单位圆，范围且在第一象限，如下图：
我们假设点（x1,y1）与X轴正半轴夹角为α,那么
{ y 2 = s i n ( α + θ ) x 2 = c o s ( α + θ ) \begin{cases} y_2=sin(α+θ)\\ x_2=cos(α+θ) \end{cases} {y2=sin(α+θ)x2=cos(α+θ)
三角函数睁开有
{ y 2 = s i n ( α ) c o s ( θ ) + c o s ( α ) s i n ( θ ) x 2 = c o s ( α ) c o s ( θ ) − s i n ( α ) s i n ( θ ) \begin{cases} y_2=sin(α)cos(θ)+cos(α)sin(θ)\\ x_2=cos(α)cos(θ)-sin(α)sin(θ) \end{cases} {y2=sin(α)cos(θ)+cos(α)sin(θ)x2=cos(α)cos(θ)−sin(α)sin(θ)
将
{ y 1 = s i n ( α ) x 1 = c o s ( α ) \begin{cases} y_1=sin(α)\\ x_1=cos(α) \end{cases} {y1=sin(α)x1=cos(α)
带入上式，有
{ y 2 = s i n ( α ) c o s ( θ ) + c o s ( α ) s i n ( θ ) = y 1 c o s ( θ ) + x 1 s i n ( θ ) = c o s ( θ ) ( y 1 + x 1 t a n ( θ ) ) x 2 = c o s ( α ) c o s ( θ ) − s i n ( α ) s i n ( θ ) = x 1 c o s ( θ ) − y 1 s i n ( θ ) = c o s ( θ ) ( x 1 − y 1 t a n ( θ ) ) \begin{cases} y_2=sin(α)cos(θ)+cos(α)sin(θ)=y_1cos(θ)+x_1sin(θ)=cos(θ)(y_1+x_1tan(θ))\\ x_2=cos(α)cos(θ)-sin(α)sin(θ)=x_1cos(θ)-y_1sin(θ)=cos(θ)(x_1-y_1tan(θ)) \end{cases} {y2=sin(α)cos(θ)+cos(α)sin(θ)=y1cos(θ)+x1sin(θ)=cos(θ)(y1+x1tan(θ))x2=cos(α)cos(θ)−sin(α)sin(θ)=x1cos(θ)−y1sin(θ)=cos(θ)(x1−y1tan(θ))
默认初始值为（1，0），记为
v 0 = [ 1 0 ] v_0=\begin{bmatrix} 1\\ 0 \end{bmatrix} v0=[10]
以上的等式可以表示为旋转矩阵的情势
[ x 2 y 2 ] = c o s ( α ) [ 1 − t a n ( α ) t a n ( α ) 1 ] [ x 1 y 1 ] \begin{bmatrix} x_2\\ y_2 \end{bmatrix} =cos(α) \begin{bmatrix} 1 & -tan(α)\\ tan(α)&1 \end{bmatrix} \begin{bmatrix} x_1\\ y_1 \end{bmatrix} [x2y2]=cos(α)[1tan(α)−tan(α)1][x1y1]
如果将令角度α，满意tan(α)=2-i, 那么就将tan和乘法运算就变成乘2的负次幂。对应在计算机中，就是移位计算。因而复杂的计算，就变成了简单的加减和移位运算。
以是我们有
[ x n y n ] = c o s ( α n ) [ 1 − 2 − n 2 − n 1 ] [ x n − 1 y n − 1 ] = c o s ( α n ) c o s ( α n − 1 ) . . c o s ( α 0 ) [ 1 − 2 − n 2 − n 1 ] [ 1 − 2 − n + 1 2 − n + 1 1 ] . . [ 1 − 2 − 0 2 − 0 1 ] [ 1 0 ] \begin{bmatrix} x_n\\ y_n \end{bmatrix} =cos(α_n) \begin{bmatrix} 1 & -2^{-n}\\ 2^{-n}&1 \end{bmatrix} \begin{bmatrix} x_{n-1}\\ y_{n-1} \end{bmatrix}= cos(α_n)cos(α_{n-1})..cos(α_0) \begin{bmatrix} 1 & -2^{-n}\\ 2^{-n}&1 \end{bmatrix} \begin{bmatrix} 1 & -2^{-n+1}\\ 2^{-n+1}&1 \end{bmatrix} .. \begin{bmatrix} 1 & -2^{-0}\\ 2^{-0}&1 \end{bmatrix} \begin{bmatrix} 1\\ 0 \end{bmatrix} [xnyn]=cos(αn)[12−n−2−n1][xn−1yn−1]=cos(αn)cos(αn−1)..cos(α0)[12−n−2−n1][12−n+1−2−n+11]..[12−0−2−01][10]
处置惩罚细节
1. 缩放因子
由前面推导，我们可以得到：
[ x n y n ] = c o s ( α n ) [ 1 − 2 − n 2 − n 1 ] [ x n − 1 y n − 1 ] = c o s ( α n ) c o s ( α n − 1 ) . . c o s ( α 0 ) [ 1 − 2 − n 2 − n 1 ] [ 1 − 2 − n + 1 2 − n + 1 1 ] . . [ 1 − 2 − 0 2 − 0 1 ] [ 1 0 ] \begin{bmatrix} x_n\\ y_n \end{bmatrix} =cos(α_n) \begin{bmatrix} 1 & -2^{-n}\\ 2^{-n}&1 \end{bmatrix} \begin{bmatrix} x_{n-1}\\ y_{n-1} \end{bmatrix}= cos(α_n)cos(α_{n-1})..cos(α_0) \begin{bmatrix} 1 & -2^{-n}\\ 2^{-n}&1 \end{bmatrix} \begin{bmatrix} 1 & -2^{-n+1}\\ 2^{-n+1}&1 \end{bmatrix} .. \begin{bmatrix} 1 & -2^{-0}\\ 2^{-0}&1 \end{bmatrix} \begin{bmatrix} 1\\ 0 \end{bmatrix} [xnyn]=cos(αn)[12−n−2−n1][xn−1yn−1]=cos(αn)cos(αn−1)..cos(α0)[12−n−2−n1][12−n+1−2−n+11]..[12−0−2−01][10]
我们可以提前将所有的cos(α)的乘积计算出来，做为一个常量，省去乘法运算，记为K
K = c o s ( α n ) c o s ( α n − 1 ) c o s ( α n − 2 ) . . . c o s ( α 0 ) = 0.60725 K=cos(α_n)cos(α_{n-1})cos(α_{n-2})...cos(α_0)=0.60725 K=cos(αn)cos(αn−1)cos(αn−2)...cos(α0)=0.60725
2. 旋转方向
通常来说CORDIC算法会引入dn ，来判断旋转方向。当前角度大于该次迭代的角度，dn为正，逆时钟旋转，反之为负，顺时针旋转。之以是会采用双旋转，是因为其通常比单向旋转的收敛性更好，结果更精确。
因而我们迭代可以写为
{ y n + 1 = y n + d n ∗ x n ∗ 2 − n x n + 1 = x n − d ∗ y n ∗ 2 − n a n g l e n + 1 = a n g l e n − d ∗ t a b l e o f a n g l e s [ n ] \begin{cases} y_{n+1}=y_n+d_n*x_n*2^{-n}\\ x_{n+1}=x_n-d*y_n*2^{-n}\\ angle_{n+1}=angle_{n}-d*tableofangles[n] \end{cases} ⎩ ⎨ ⎧yn+1=yn+dn∗xn∗2−nxn+1=xn−d∗yn∗2−nanglen+1=anglen−d∗tableofangles[n]
table_of_angles 存储的是θ值， θn=arctan(2-n);
对应的表格如下：
n2^(-n)arctan(2^(-n))010.78539816310.50.46364760920.250.24497866330.1250.12435499540.06250.0624188150.03150.03123983360.0156250.01562372970.00781250.00781234180.003906250.0039062390.0019531250.001953123100.0009765630.000976562110.0004882810.000488281120.0002441410.000244141130.000122070.00012207146.10352E-056.10352E-05153.05176E-053.05176E-05 下面我会手把手领导大家从软件建模到硬件实现CORDIC算法，规定输入和输出都是无符号17位数，1位整数位，16位小数位。
Python 代码
测试代码
- #初始化部分，定义参数
- import math
- from math import floor
- NUM_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.0
- for 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=k
- y=0
- angle_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_next
- y=y_next
- return 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 False
- return 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 begin
- angles[0] = 17'h0C910;
- angles[1] = 17'h076B2;
- angles[2] = 17'h03EB7;
- angles[3] = 17'h01FD6;// i=0~3
- angles[4] = 17'h00FFB;
- angles[5] = 17'h007FF;
- angles[6] = 17'h00400;
- angles[7] = 17'h00200;// i=4~7
- angles[8] = 17'h00100;
- angles[9] = 17'h00080;
- angles[10] = 17'h00040;
- angles[11] = 17'h00020;// i=8~11
- angles[12] = 17'h00010;
- angles[13] = 17'h00008;
- angles[14] = 17'h00004;
- angles[15] = 17'h00002;// i=12~15
- angles[16] = 17'h00001;
- end
- reg 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)begin
- if(!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]};
- end
- x=x_next;
- y=y_next;
- end
- end
- assign sin_out = y[32:16];
- assign cos_out = x[32:16];
- endmodule
- `timescale 1ns / 1ps
- module Cordic_Sin_Test();
-
- reg [16:0] theta;
- wire [16:0] sin_out;
- wire [16:0] cos_out;
- initial begin
-
- theta=17'h0;
- #10;
- //15
- theta=17'h4305;
- #10;
- //30
- theta=17'h860A;
- #10;
- //45
- theta=17'hC90F;
- #10;
- //60
- theta=17'h10C15;
- #10;
- //75
- theta=17'h14F1A;
- #10;
- //90
- theta=17'h1921F;
- #10;
-
- end
- Cordic_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°0655361546553615°1696263303169626330230°3276856756327685675545°4634146341463404634160°5675632768567553276875°6330316962633021696290°65536065536154 除了个别点外的绝对误差比力大外，其余的计算精度相称高，误差很小
免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。 |
