目录
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;
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个，较之前有很大程度改善。
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
|
|
384
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
