目录
1、 简介
1.1 系统的目的
1.2 系统的背景
2、 需求概括
2.1 系统需求
2.2 当前系统问题
3、 建议的系统
3.1 设计重点
3.2 系统的原理
3.2.1 算法框图
3.3 数据预处理部分
3.4 迭代算法的核心部分
3.5 收敛判断
3.6 输出格式化
3.7 模块接口信号
3.9 仿真测试
4、 系统分析
4.1 精度分析
4.2 资源分析
4.3 时延分析
5、代码
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、 建议的系统
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;
end3.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个，较之前有很大程度改善。文章来源：https://www.toymoban.com/news/detail-819210.html
4.3 时延分析
本设计固定延迟32个cycle，但是在大多数情况下，10个周期内就能够收敛了，所以可以根据实际情况做更改，在时延和计算精度上做折中文章来源地址https://www.toymoban.com/news/detail-819210.html
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
|
|
到了这里，关于基于FPGA的求模运算器的文章就介绍完了。如果您还想了解更多内容，请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章，希望大家以后多多支持TOY模板网！
