目录
1、 简介
1.1 系统的目的
在实际算法中,有很多地方需要用到求模的算法,以及开根的数学计算。本设计在FPGA上实现一款基于Cordic算法的两个向量的求模计算方法,旨在替代传统的计算公式,优化面积和速度。
1.2 系统的背景
在电机控制算法中,有很多地方需要用到开根和平方差开根的运算,然而FPGA内部却没有专门的硬件资源来实现这种算法。目前系统使用的方案是a2+b2然后开根,开根使用的是xilinx的Cordic IP,Xilinx提供的IP保密性太高,使用者无法进行项目FPGA平台的灵活切换。
2、 需求概括
2.1 系统需求
在FPGA开发板上实现一种面积小,性能高,精度满足要求的求模计算模块,要求自主源码实现,后期易于切换平台,易于维护。
其它要求:
数据位宽8-32位可配
主频:100MHz
延迟: 小于32 cycle
支持定点小数,小数位可配
2.2 当前系统问题
同提出的设计相比,目前的系统普遍存在以下问题:
- 计算精度低
- 资源消耗高
- 可移植性差
- 计算延迟高
3.1 设计重点
- 通过提高内部运算数据宽度和舍入方法来提高运算精度。
- 将各参数预留有接口,方面移植到各种设计场合。
- 利用Cordic圆周旋转方式进行结果逼近。
- 输入与输出数据的格式化。
3.2 系统的原理
CORDIC算法是一个“化繁为简”的算法,将许多复杂的运算转化为一种“仅需要移位和加法”的迭代操作。CORDIC算法有旋转和向量两个模式,分别可以在圆坐标系、线性坐标系和双曲线坐标系使用,从而可以演算出8种运算,而结合这8种运算也可以衍生出其他许多运算。在这里我们只介绍求模运算在CORDIC算法中实现的条件。
CORDIC算法有两种工作模式:旋转模式和向量模式,旋转模式是给定旋转角,计算旋转后的坐标,而向量模式相当于旋转模式的变体(倒推得到),根据旋转后的坐标,得到向量的角度和模值。
由于CORDIC算法只有移位和加法运算,所以比较适合于在硬件上实现。
如上图所示,若已知点矢量终点A0(X0,Y0),若将该矢量逆时针旋转θ 可以根据三角运算得到B0(X1,Y1)点坐标:
令θ1 = − θ即顺时针旋转θ角度,则:
联立上述两个式子,引入常数$ d (d=-1,+1)$ ,因此可得:
这个算法的核心在于将一系列已知的t a nθ作为表格键值进行存储,而t a nθ可以约等于1/2-n并且1/2-n在FPGA中可以通过右移进行快速运算。t a nθ各个已知存储值如下:
i | θ | t a nθ | cosθ | accosθ |
0 | 45 | 1 | 0.707106781186548 | 0.707106781186548 |
1 | 25.56505 | 0.50 | 0.894427190999916 | 0.632455532033676 |
2 | 14.03243 | 0.25 | 0.970142500145332 | 0.613571991077897 |
3 | 7.125016 | 0.125000000000000 | 0.992277876713668 | 0.608833912517753 |
4 | 3.576334 | 0.0625000000000000 | 0.998052578482889 | 0.607648256256168 |
5 | 1.789910 | 0.0312500000000000 | 0.999512076087079 | 0.607351770141296 |
6 | 0.895173 | 0.0156250000000000 | 0.999877952034695 | 0.607277644093526 |
7 | 0.447614 | 0.00781250000000000 | 0.999969483818788 | 0.607259112298893 |
8 | 0.223810 | 0.00390625000000000 | 0.999992370692779 | 0.607254479332563 |
9 | 0.111905 | 0.00195312500000000 | 0.999998092656824 | 0.607253321089875 |
10 | 0.055952 | 0.000976562500000000 | 0.999999523163183 | 0.607253031529135 |
11 | 0.027976 | 0.000488281250000000 | 0.999999880790732 | 0.607252959138945 |
12 | 0.013988 | 0.000244140625000000 | 0.999999970197679 | 0.607252941041397 |
13 | 0.006994 | 0.000122070312500000 | 0.999999992549420 | 0.607252936517011 |
14 | 0.003497 | 6.10351562500000e-05 | 0.999999998137355 | 0.607252935385914 |
15 | 0.001748 | 3.05175781250000e-05 | 0.999999999534339 | 0.607252935103140 |
而多次旋转过程中,每次旋转的c o sθ需要连续相乘,而多次相乘极限也趋近与0.607252这一个常数,因此也可做近似处理。那么现在还有最后一个问题,这一系列角度能够通过多次旋转得到任意的角度吗?可以看到每个角度是不断降低减半的,呈递减的分布,从宏观上观察大致是可以进行趋近到某一个常数的。
cordic算法向量模式,已知点坐标(x0,y0),可以求得该向量的角度即arctan(y0/x0)。这种可以理解为需要通过多次旋转,将该向量旋转至x轴上,即y0 = 0,此时旋转过的角度即为向量角度,x最终坐标即为向量的长度。
最终迭代公式如下:
只要y趋向于0,那么x的值就是最终所求的模长。
3.2.1 算法框图
系统的算法框图如图1所示,主要分为数据预处理模块,旋转迭代模块,收敛判断,输出格式化。
3.3 数据预处理部分
数据预处理模块主要负责将输入的数据变换到第一区间,即0-0.5pi的区间。
由于求模运算,都是基于正数去算的,所以需要将负数求绝对值,这样下来输入数据都为正,确保能够落在第一区间。针对特殊情况,x为0或者y为0的,可以直接给出模值,这里就不再详述了。
关键代码:
直接给出模值得情况
If(x==0)
Result <= y;
Else if(y==0)
Result <=x;
Else
Result <=r_result;
变换到第一区间的代码:
r_x <= abs(x);
r_y <= abs(y);
3.4 迭代算法的核心部分
设置初始值让x0=data_a,y0=data_b;
根据以上3.2节迭代公式,每次迭代计算出x和y的值,还有d的值,判断下一次运算是加还是减。
核心代码如下:
Assign s_dx=r_Z[DWIDTH-1]==1’B0? r_data_y: - r_data_y;
Assign s_dy=r_Z[DWIDTH-1]==1’B0? r_data_x: - r_data_x;
Assign s_dz=r_Z[DWIDTH-1]==1’B0? actan: - actan;
If(cvt_end == 1’b1)begin
R_data_x <= r_data_x;
R_data_y <= r_data_y;
R_Z <= R_Z;
End else begin
R_data_x <= r_data_x – (s_dx>>>cycle);
R_data_y <= r_data_y+ (s_dy>>>cycle);;
R_Z <= R_Z-s_dz;
end
3.5 收敛判断
根据大量的测试数据,统计出迭代需要的次数以及精度之间的关系设置一个最大迭代次数,达到这个次数就停止计算,输出结果。
If(start_pulse| r_cvt_end)
Cycle <= 0;
Else
Cycle <= Cycle + 1’b1;
If(start_pulse | r_cvt_end)
Reach_max <= 0;
Else if(cycle == MAX_CYCLE)
Reach_max <= 1;
第二种方案是判断角度值是否小于一个接近于0的数,如果条件满足则停止计算,输出结果。
If(start_pulse|cvt_end)
r_cvt_end <= 0;
else if(cycle >8’d5 && abs(r_Z)<3)
r_cvt_end <= 1;
有一种特殊情况就是迭代一次或者极少的次数就收敛了,经过实际测量,这种情况测出来得数据不准确,误差很大,可能跟算法的最终补偿因子有关系。这里我们跳过迭代次数小于5的情况,让其继续迭代计算。
3.6 输出格式化
到目前为止,输出的数据并不是最终的结果,因为每次旋转的时候都少算了一个乘数cosθ,所以我们算出来的结果需要补偿回来。
3.2.2节中提到过这个补偿值大约为0.607252,所以我们需要将最终的结果乘以这个值再输出。
输出模长信号需要经过用户要求的格式去整理一下才能输出,不然内部运算的数据位宽和小数点位置与输出要求的不一样,后级很容易发生错误。
首先由用户输入的参数来计算出输出数据的位宽和小数位置
例如默认配置下:
Parameter FRACTION_BIT =20
Parameter DIN_WIDTH =32
Parameter DOUT_WIDTH =32
表示输入数据位宽是32位,其中低20为为小数,即Q20
假如系统为了提高运算精度,在迭代的时候将运算数据小数位扩展到50位,整数位不变,则迭代收敛后,所得出来的数据62位宽,50位小数。
在我们输出时候需要根据输出的参数DOUT_WIDTH及FRACTION_BIT整理输出的数据格式,直接对最后的30bit进行四舍五入:
R_result<= r_data_x[30+ DIN_WIDTH-1:30]+ r_data_x[29];
3.7 模块接口信号
参数 | 方向 | 位宽 | 描述 |
DIN_WIDTH | 整型 | 输入位宽 | |
DOUT_WIDTH | 整型 | 输出位宽 | |
FRACTION_BIT | 整型 | 小数位数 |
信号 | 方向 | 位宽 | 描述 |
Clk | In | 1 | 时钟输入 |
Rst | In | 1 | 复位输入 |
Start | In | 1 | 开始计算脉冲指令 |
Data_a | In | X | 输入数据a |
Data_b | In | X | 输入数据a |
Result | Out | X | 输出数据 |
Ovalid | Out | 1 | 输出数据有效信号 |
3.9 仿真测试
基于vivado下仿真出来的波形如下:
计算开始时的仿真波形
计算收敛后的仿真波形
从图上看,给的输入信号是3和4,计算出来的模长是5,非常精确,收敛速度也非常快,在第5个cycle时候已经达到4.996了。
其中start脉冲过来后,busy开始拉高,计算开始,数据收敛后ovalid信号拉高一次,busy信号拉低。
中间数据位宽比较大,62位的,本设计为了提高计算精度,过程中小数位采用50位定点小数。
4、 系统分析
4.1 精度分析
输入1 | 输入2 | 计算结果 | 理论结果 | 精度 |
17.071110725 | 20.328888893 | 26.545932769 | 26.5459327395 | 1.11e-9 |
19.071110725 | 22.328888893 | 29.364715576 | 29.3647159612 | 1.31e-8 |
21.071110725 | 24.328888893 | 32.18519115448 | 32.1851913456 | 5.9e-9 |
23.071110725 | 26.328888893 | 35.006949424 | 35.0069498875 | 1.32e-8 |
25.071110725 | 28.328888893 | 37.829704284668 | 37.8297044516 | 4.4e-9 |
27.071110725 | 30.328888893 | 40.653247833252 | 40.6532475624 | 6.66e-9 |
29.071110725 | 32.328888893 | 43.477425575 | 43.4774255889 | 3.197e-10 |
31.071110725 | 34.328888893 | 46.302122116 | 46.3021223521 | 5.1e-9 |
有以上测试数据可知,精度最低已经达到1.32e-8,已经完全满足大部分场景了。
4.2 资源分析
系统在VIVADO平台下选择xczu5eg系列的FPGA,并将本模块设为顶层,综合后资源报告如下:
因为节省了da*da+db*db的运算,所以节省了两大位宽乘法器和一个加法运算,目前lut消耗了585个 dsp 为4个,较之前有很大程度改善。
4.3 时延分析
本设计固定延迟32个cycle,但是在大多数情况下,10个周期内就能够收敛了,所以可以根据实际情况做更改,在时延和计算精度上做折中
5、代码
Module sincos_cordic #(
Parameter MAX_CYCLE=8’d40,
Parameter FRACTION_BIT=20,
Parameter DIN_WIDTH=12,
Parameter DOUT_WIDTH=22
)(
Input wire clk,
Input wire ce,
Input wire reset,
Input wire start,
Input wire [DIN_WIDTH-1:0] data_a,
Input wire [DIN_WIDTH-1:0] data_b,
output wire [DOUT_WIDTH-1:0] result,
output wire busy,
output wire ovalid
);
Localparam DWIDTH = 30+DIN_WIDTH;
Localparam AMP_WIDTH = DWIDTH- DIN_WIDTH;
Localparam Kn = =32’h26dd3b6a;
Reg r_start_sft;
Reg r_start_pulse;
Reg r_busy;
Reg [DIN_WIDTH-1:0] r_data_a;
Reg [DIN_WIDTH-1:0] r_data_b;
Reg signed [DWIDTH -1:0] r_data_x;
Reg signed [DWIDTH -1:0] r_data_y;
Wire signed [DWIDTH -1:0] s_dx;
Wire signed [DWIDTH -1:0] s_dy;
Reg [7:0] cycle;
Reg [7:0] repeat_point;
Wire cvt_end;
Reg r_busy1;
Reg [DOUT_WIDTH-1:0] r_result;
Reg reach_max=0;
Function [DWIDTH -1:0] abs;
Input [DWIDTH -1:0] din;
Abs= (din[DWIDTH -1:0])==1’b1)?-din:din;
Endfuction
always@(posedge clk) begin
r_data_a <= abs(data_a);
r_data_b <= abs(data_b);
end
always@(*)
r_start_pulse <= start & ~r_start_sft;
always@(posedge clk) begin
if(reset == 1’b1)
cycle <= 0;
else begin
if(r_start_pulse | cvt_end)
cycle <= 0;
else
cycle <= cycle + 1’b1;
end
end
always@(posedge clk) begin
if(r_start_pulse)
reach_max <= 1’b0;
else if(cycle == MAX_CYCLE)
reach_max <= 1’b1;
else
reach_max <= reach_max;
end
assign cvt_end = reach_max;
assign s_dx = (r_data_y[DWIDTH-1]==1’B1)?r_data_y:- r_data_y;
assign s_dy = (r_data_y[DWIDTH-1]==1’B1)?r_data_x:- r_data_x;
always@(posedge clk) begin
if(r_start_pulse)begin
r_data_x <= { r_data_a ,{AMP_WIDTH{1’b0}}};
r_data_x <= { r_data_b ,{AMP_WIDTH{1’b0}}};
end else if(cvt_end == 1’b1)begin
r_data_x <= r_data_x;
r_data_y <= r_data_y;
end else begin
r_data_x <= r_data_x – (s_dx>>>cycle);
r_data_y <= r_data_y + (s_dy>>>cycle);
end
end
always@(posedge clk) begin
if(reset== 1’b1)begin
r_busy <= 1’b0;
r_busy1 <= 1’b0;
end else begin
r_busy1 <= r_busy;
if(r_start_pulse)
r_busy <= 1’b1;
else if(cvt_end == 1;b1)
r_busy <= 1’b0;
else
r_busy <= r_busy;
end
end
wire [DWIDTH-1:0] s_result;
wire [DWIDTH+AMP_WIDTH-1:0] s_result0;
assign s_result0 = r_data_x*kn;
assign s_result = s + result0[DWIDTH + AMP_WIDTH-1: AMP_WIDTH];
always@(posedge clk)
r_result <= s_result [DWIDTH -1: AMP_WIDTH]+ s_result [AMP_WIDTH-1];
assign busy = r_busy;
assign ovalid = r_busy1&~r_busy;
assign result = r_result;
endmodule
|