HLS / Chisel 实现CORDIC算法高性能复数除法

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三种模式实现复数除法

复数除法高性能计算

首先有如下的复数除法公式，我们一般会做如下的拆分：

$$\mathrm{Z}=\frac{a+ib}{p+iq}=\frac{ap+bq}{p^2+q^2}+i\frac{bp-aq}{p^2+q^2}$$

1. CV模式转化为极坐标形式

首先输入如下的数值

$${\mathrm{X}}_0=p,Y_0=q,Z_0=0$$

在此反三角函数计算模块可以得到如下的结果：

$$X=K\sqrt{p^2+q^2}=r,\mathrm{Y}=0,Z=\tan (q/p{)}^{-1}=\theta$$

其中K在CV模式的推导文章中有介绍，是一个固定的常数。接下来的步骤我们需要极坐标系的结果$(r,\theta )$。

2. LV模式计算实数除法

将1中的结果 $X=r$作为输入此模块，输入初始化为如下：

$${\mathrm{X}}_0=K\sqrt{p^2+q^2},Y_0=a,Z_0=0$$

最后通过实数除法得到如下的结果：

$$X=K\sqrt{p^2+q^2},\mathrm{Y}=0,Z=\frac{a}{K\sqrt{p^2+q^2}}$$

同理再输入

$${\mathrm{X}}_0=K\sqrt{p^2+q^2},Y_0=b,Z_0=0$$

获得如下结果：

$$X=K\sqrt{p^2+q^2},\mathrm{Y}=0,Z=\frac{b}{K\sqrt{p^2+q^2}}$$

3. CR模式计算旋转推导出复数除法

将1的结果$Z=\tan (q/p{)}^{-1}=\theta$，2中的结果$Z=\frac{a}{K\sqrt{p^2+q^2}}，Z=\frac{b}{K\sqrt{p^2+q^2}}$输入CR模块，输入如下：

$${\mathrm{X}}_0=\frac{a}{K\sqrt{p^2+q^2}},Y_0=\frac{b}{K\sqrt{p^2+q^2}},Z_0=-\theta$$

已知CR模式的推导公式如下：

$$\begin{array}{rl}& K(x\cos z-y\sin z)\\ & K(y\cos z+x\sin z)\end{array}$$

已知

$$\tan (q/p{)}^{-1}=\theta$$

那么可以得到：

$$\begin{array}{rl}& \sin \theta =\frac{q}{\sqrt{q^2+p^2}}\\ & \cos \theta =\frac{p}{\sqrt{q^2+p^2}}\end{array}$$

将x,y,z带入CR模式的推导模式则可以得到：

$$X=\frac{ap+bq}{p^2+q^2}Y=\frac{bp-aq}{p^2+q^2}$$

即复数的实部和虚部计算完成。

输入输出范围限定

整个CORDIC计算复数除法的输入输出如上图。由于CR、CV、LV模块有输入输出的范围限制，则对于整个复数除法计算大模块的输入输出也有范围限制，首先我们单独讨论子模块：

• CR模块执行旋转功能，将输入的点（x,y）∈R旋转$\theta \in [-90°，90°]$角度得到点（X,Y）
• CV模块执行极坐标转化（arctan计算），将输入的点（x>0,y∈R）即一四象限的点转化为极坐标（r∈R,$\theta \in [-90°，90°]$）
• LV模块执行除法运算，在我们设计的模块中，我们使用二进制移位处理使得可以实现实数除法，所以 输入（x,y）∈R，输出z∈R

基于如上分析，我们分别对复数除法模块的a、b、q、p进行取值范围分析：

1. 对于a、b

a、b参与的运算从LV模块开始，然后LV模块输出影响CR模块的输入。由于LV和CR的输入都是实数域R，所以输入a、b∈R

2. 对于p、q

p、q参与了CV、CR、LV三个模块的运算，有以下分析。

a. 对于q<0，q>0

此时在首先CV模块的输入，(q,p)是在第二象限，由于CV模块的输入范围限制，我们可以将$(q,p)$作原点对称到$(q‘,p‘)=(-p,-q)$,则计算出来的结果$\theta$为(q`,p`) 对应的结果，但极坐标的$K\sqrt{p^2+q^2}=r$是没有任何影响的

然后在LV模块中计算由于r是不变化的，所以LV模块没有任何影响。

最后在CR模块中将$-\theta$带入，得到的应该是

$$X=\frac{ap‘+bq‘}{p{‘}^2+q{‘}^2}Y=\frac{bp‘-aq‘}{p{‘}^2+q{‘}^2}$$

带入$(q‘,p‘)=(-p,-q)$即可以将负号加上就可以得到复数除法结果

$$X=\frac{-ap-bq}{p^2+q^2}Y=\frac{-bp+aq}{p^2+q^2}$$

$$\mathrm{Z}=\frac{a+ib}{p+iq}=\frac{ap+bq}{p^2+q^2}+i\frac{bp-aq}{p^2+q^2}=-\frac{ap‘+bq‘}{p{‘}^2+q{‘}^2}-i\frac{bp‘-aq‘}{p{‘}^2+q{‘}^2}$$

b.对于q<0，q<0

此时在首先CV模块的输入，(q,p)是在第三象限，同理如上小节方法也可以转换到第一象限，最后在第三步CR模块得到的处理也一样，即如下:

$$X=\frac{-ap-bq}{p^2+q^2}Y=\frac{-bp+aq}{p^2+q^2}$$

c.对于q>0，q>0 或 q>0，q<0

此时首先CV模块的输入是在允许的一四象限处理范围内，可以直接计算得到正确的结果

$$X=\frac{ap+bq}{p^2+q^2}Y=\frac{bp-aq}{p^2+q^2}$$

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.