CORDIC(坐标旋转数字算法)是一种计算三角、双曲和其他数学函数的数字算法,每次运算均产生一次结果输出。这使我们能够根据应用需求调整算法精度;增加运算迭代次数可以得到更精确的结果。CORDIC 是只使用加法、减法、移位和查找表实现的简单算法,这种算法在FPGA中实现效率高,在硬件算法实现中经常用到。
论文[1]介绍了一种通过CORDIC算法中CR、CV、LV三种模式拆分高性能计算复数除法的方法,本文将使用Chisel 实现该方法,这里默认在以下几篇文章的CORDIC原理基础之上完成。
HLS / Chisel 实现CORDIC算法计算三角函数(圆周系统之旋转模式)HLS / Chisel 实现CORDIC算法计算反三角函数(圆坐标系向量模式)HLS / Chisel 实现CORDIC算法计算除法(线性坐标系向量模式LV)
原理
CORDIC背景
CORDIC算法有旋转和向量两个模式,分别可以在圆坐标系、线性坐标系和双曲线坐标系使用,从而可以演算出8种运算,而结合这8种运算也可以衍生出其他许多运算。下表展示了8种运算在CORDIC算法中计算的核心公式推导:
接下来我们将利用如上的CR、CV、LV三种模式实现复数除法
复数除法高性能计算
首先有如下的复数除法公式,我们一般会做如下的拆分:
- CV模式转化为极坐标形式
首先输入如下的数值
在此反三角函数计算模块可以得到如下的结果:
其中K在CV模式的推导文章中有介绍,是一个固定的常数。接下来的步骤我们需要极坐标系的结果。
- LV模式计算实数除法
将1中的结果 作为输入此模块,输入初始化为如下:
最后通过实数除法得到如下的结果:
同理再输入
获得如下结果:
- CR模式计算旋转推导出复数除法
将1的结果,2中的结果输入CR模块,输入如下:
已知CR模式的推导公式如下:
已知
那么可以得到:
将x,y,z带入CR模式的推导模式则可以得到:
即复数的实部和虚部计算完成。
输入输出范围限定
整个CORDIC计算复数除法的输入输出如上图。由于CR、CV、LV模块有输入输出的范围限制,则对于整个复数除法计算大模块的输入输出也有范围限制,首先我们单独讨论子模块:
- CR模块执行旋转功能,将输入的点(x,y)∈R旋转角度得到点(X,Y)
- CV模块执行极坐标转化(arctan计算),将输入的点(x>0,y∈R)即一四象限的点转化为极坐标(r∈R,)
- LV模块执行除法运算,在我们设计的模块中,我们使用二进制移位处理使得可以实现实数除法,所以 输入(x,y)∈R,输出z∈R
基于如上分析,我们分别对复数除法模块的a、b、q、p进行取值范围分析:
对于a、b
a、b参与的运算从LV模块开始,然后LV模块输出影响CR模块的输入。由于LV和CR的输入都是实数域R,所以输入a、b∈R
对于p、q
p、q参与了CV、CR、LV三个模块的运算,有以下分析。
a. 对于q<0,q>0
此时在首先CV模块的输入,(q,p)是在第二象限,由于CV模块的输入范围限制,我们可以将作原点对称到,则计算出来的结果为(q`,p`) 对应的结果,但极坐标的是没有任何影响的
然后在LV模块中计算由于r是不变化的,所以LV模块没有任何影响。
最后在CR模块中将带入,得到的应该是
带入即可以将负号加上就可以得到复数除法结果
b.对于q<0,q<0
此时在首先CV模块的输入,(q,p)是在第三象限,同理如上小节方法也可以转换到第一象限,最后在第三步CR模块得到的处理也一样,即如下:
c.对于q>0,q>0 或 q>0,q<0
此时首先CV模块的输入是在允许的一四象限处理范围内,可以直接计算得到正确的结果
HLS实现
- 定义复数
- 对输入的数做处理转移到输入合适的区间
#define COS_SIN_TYPE float
#define NUM_ITERATIONS 50
typedef struct { COS_SIN_TYPE re, im; } Complex;
void cordic_complex_divide(Complex dividend, Complex divisor, Complex &z)
{
COS_SIN_TYPE a,b,p,q;
a = dividend.re;
b = dividend.im;
p = divisor.re;
q = divisor.im;
bool flag = 1; // 默认不改变最后输出实部虚部的符号
if((p < 0 and q > 0) or (p < 0 and q < 0)){
// 将2,3象限的点旋转到1,4象限,方便CV模块操作
p = -p;
q = -q;
flag = 0;
}
COS_SIN_TYPE theta, r, arctan; // arctan为舍弃的
cordic_cv(p, q, r, theta, arctan);
COS_SIN_TYPE temp_divide_a, temp_divide_b;
cordic_lv(r, a, temp_divide_a);
cordic_lv(r, b, temp_divide_b);
cordic_cr(x = temp_divide_a, y = temp_divide_b, theta = -theta, z.re, z.im);
if(flag) {
z.re = -z.re;
z.im = -z.im;
}
}
Chisel实现
定义数据类型
import chisel3._
import chisel3.experimental._
trait HasDataConfig {
/* 定点数的位宽定义 */
val DataWidth = 32
val BinaryPoint = 20
}
class Complex extends Bundle with HasDataConfig{
/* 复数的 实部,虚部*/
val re: FixedPoint = FixedPoint(DataWidth.W, BinaryPoint.BP)
val im: FixedPoint = FixedPoint(DataWidth.W, BinaryPoint.BP)
}
class ComplexOperationIO extends Bundle with HasDataConfig{
val op1: Complex = Input(new Complex())
val op2: Complex = Input(new Complex())
val res: Complex = Output(new Complex())
}
源码
import chisel3._
import chisel3.experimental._
class cordic_complex_divide extends Module with HasDataConfig with Cordic_Const_Data {
/*
* @io.op1 : 输入的被除数 定点数表示的复数Complex 输入范围R
* @io.op2 : 输入的除数 定点数表示的复数Complex 输入范围R
* @io.res : 输出的复数除法结果 定点数表示的复数Complex 范围R
* details: 利用cordic实现复数除法,周期为3*迭代次数NUM_ITERATIONS
**/
val io = IO(new ComplexOperationIO)
val a = WireInit(io.op1.re)
val b = WireInit(io.op1.im)
val p = Wire(FixedPoint(DataWidth.W,BinaryPoint.BP))
val q = Wire(FixedPoint(DataWidth.W,BinaryPoint.BP))
/* (p,q)所在象限必须调整到一四象限才能在CV阶段输入,最后通过flag再回调即可,flag 需要存储NUM_ITERATIONS个流水层级到最后输出时使用 */
val flag_reg: Vec[Bool] = RegInit(VecInit(Seq.fill(NUM_ITERATIONS*3)(1.B)))
when(io.op2.re < FixedPoint.fromDouble(0,DataWidth.W,BinaryPoint.BP)){
p := -io.op2.re
q := -io.op2.im
flag_reg(0) := 0.B
}.otherwise{
p := io.op2.re
q := io.op2.im
flag_reg(0) := 1.B
}
for(i <- 1 until NUM_ITERATIONS*3){
flag_reg(i) := flag_reg(i-1)
}
/* a,b为LV模块的输入,但是需要存储保存NUM_ITERATIONS个流水层级到LV模块使用 */
val a_reg: Vec[FixedPoint] = RegInit(VecInit(Seq.fill(NUM_ITERATIONS)(FixedPoint.fromDouble(0.0, DataWidth.W, BinaryPoint.BP))))
val b_reg: Vec[FixedPoint] = RegInit(VecInit(Seq.fill(NUM_ITERATIONS)(FixedPoint.fromDouble(0.0, DataWidth.W, BinaryPoint.BP))))
for(i <- 0 until NUM_ITERATIONS){
if(i == 0){
a_reg(0) := a
b_reg(0) := b
}else{
a_reg(i) := a_reg(i-1)
b_reg(i) := b_reg(i-1)
}
}
/* CV模块得到极坐标坐标系结果 */
val (theta, r, arctan) = cordic_cv(p,q)
/* theta为CV模块的输出,但是需要存储保存NUM_ITERATIONS个流水层级到CR模块使用 */
val theta_reg: Vec[FixedPoint] = RegInit(VecInit(Seq.fill(NUM_ITERATIONS)(FixedPoint.fromDouble(0.0, DataWidth.W, BinaryPoint.BP))))
for(i <- 0 until NUM_ITERATIONS){
if(i == 0){
theta_reg(0) := theta
}else{
theta_reg(i) := theta_reg(i-1)
}
}
/* LV模块并行计算得到两个除法结果 */
val z_a = cordic_divide(r, a_reg(NUM_ITERATIONS-1))
val z_b = cordic_divide(r, b_reg(NUM_ITERATIONS-1))
/* CR旋转点计算复数除法 */
val (re, im) = cordic_cr(z_a, z_b, -theta_reg(NUM_ITERATIONS-1))
when(flag_reg(NUM_ITERATIONS*3-1)){
io.res.re := re
io.res.im := im
}.otherwise{
io.res.re := -re
io.res.im := -im
}
}
object cordic_complex_divide_APP extends App {
println(getVerilogString(new cordic_complex_divide))
(new chisel3.stage.ChiselStage).emitVerilog(new cordic_complex_divide)
}
定义伴生对象工厂方法
最后我们定义两个个伴生对象工厂方法来方便直接调用函数,无需多加连线
/*
* 定义这个类的伴生对象,并定义一个工厂方法来简化模块的例化和连线。
**/
object cordic_complex_divide {
def apply(divided: Complex, divisor: Complex): Complex = {
/*
* @divided : 输入的被除数 定点数表示的复数Complex 输入范围R
* @divisor : 输入的除数 定点数表示的复数Complex 输入范围R
* @res : 输出的复数除法结果 定点数表示的复数Complex 范围R
* details: 利用cordic实现复数除法,周期为3*迭代次数NUM_ITERATIONS
**/
val unit = Module(new cordic_complex_divide)
unit.io.op1 := divided
unit.io.op2 := divisor
unit.io.res
}
}
参考论文
[1] Yang B , Dong W , Liu L . Complex division and square-root using CORDIC. IEEE, 2012.