CORDIC算法：三角函数的硬件加速革命——从数学原理到FPGA实现的超高效计算方案
计算机该怎样求解三角函数？或许你的第一印象是采用泰勒睁开，大概采用多项式进行逼近。对于前者，往返的迭代计算开销本钱很大；对于后者，多项式式逼近在较窄的范围內比力接近，超过肯定范围后，就变得十分不抱负了。例如x–>0时，x~sin(x)今天，我们将要先容三角函数的另一种求法：CORDIC算法
原理
CORDIC的算法的核心就是通过迭代，利用三角函数的物理性质，不绝累积旋转角度，从而得到所求角度的精确近似。
我们假设圆为单位圆，范围且在第一象限，如下图：
https://i-blog.csdnimg.cn/direct/620ef9c24a3144d7a9e1a7e460f24dfa.png#pic_center
我们假设点（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=
以上的等式可以表示为旋转矩阵的情势
[ 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} =cos(α)
如果将令角度α，满意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} =cos(αn)=cos(αn)cos(αn−1)..cos(α0)..
处置惩罚细节
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} =cos(αn)=cos(αn)cos(αn−1)..cos(α0)..
我们可以提前将所有的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 \end{cases} ⎩ ⎨ ⎧yn+1=yn+dn∗xn∗2−nxn+1=xn−d∗yn∗2−nanglen+1=anglen−d∗tableofangles
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)))
returnfloor(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
else:
x_next=x+y*2**(-i)
y_next=y-x*2**(-i)
angle_temp+=Angles_Table
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-test)>3):
print("Error! Loss of accuracy! ground_truth: %f, test: %f", ground_truth, test)
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()
比力结果
https://i-blog.csdnimg.cn/direct/cff78cc76ff44dc8bba7f8ce1e5c901e.png#pic_center
由此可知，CORDIC算法精度很高
Verilog 代码
模块代码
module Cordic_Sin(
input wire theta, // 输入角度（Q1.16格式，范围0 ~ π/2）
output wire sin_out, // 输出sin值（Q1.16格式）
output wire cos_out
);
// 预计算参数（Q1.16格式）
localparam signed K =17'sh09B75;// 1/1.64676补偿因子;17'h1A592;
//Q1.15
reg signed angles ; //arctan(2^-i)
integer iter;
initial begin
angles= 17'h0C910;
angles= 17'h076B2;
angles= 17'h03EB7;
angles= 17'h01FD6;// i=0~3
angles= 17'h00FFB;
angles= 17'h007FF;
angles= 17'h00400;
angles= 17'h00200;// i=4~7
angles= 17'h00100;
angles= 17'h00080;
angles = 17'h00040;
angles = 17'h00020;// i=8~11
angles = 17'h00010;
angles = 17'h00008;
angles = 17'h00004;
angles = 17'h00002;// i=12~15
angles = 17'h00001;
end
reg signed x,y; //初始化 x=K; y=0
reg signed x_next,y_next;
reg signed 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)begin //
//正向旋转
x_next=x-(y>>>i);//算术移位
y_next=y+(x>>>i);
angle=angle-{1'b0,angles};
end else begin
//负向旋转
x_next=x+(y>>>i);
y_next=y-(x>>>i);
angle=angle+{1'b0,angles};
end
x=x_next;
y=y_next;
end
end
assign sin_out = y;
assign cos_out = x;
endmodule
Test Bench
`timescale 1ns / 1ps
module Cordic_Sin_Test();
reg theta;
wire sin_out;
wire 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
仿真结果
https://i-blog.csdnimg.cn/direct/814cfae8081a498ba4a6dd11e4992884.png#pic_center
以上的数据，输入数据需要将其转换成弧度值，然后转换成S0I1F16定点格式，sin,cos正确值也是一样
输入角度sin正确值cos正确值sin计算值cos计算值0°0655361546553615°1696263303169626330230°3276856756327685675545°4634146341463404634160°5675632768567553276875°6330316962633021696290°65536065536154 除了个别点外的绝对误差比力大外，其余的计算精度相称高，误差很小
免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。
页:
[1]