HLS / Chisel 实践CORDIC高性能计算复数平方根

CORDIC(坐标旋转数字算法)是一种计算三角、双曲和其他数学函数的数字算法，每次运算均产生一次结果输出。这使我们能够根据应用需求调整算法精度；增加运算迭代次数可以得到更精确的结果。CORDIC 是只使用加法、减法、移位和查找表实现的简单算法，这种算法在FPGA中实现效率高，在硬件算法实现中经常用到。

论文[1]介绍了一种通过CORDIC算法中CR、CV、HV三种模式拆分高性能计算复数平方根的方法，本文将使用Chisel 实现该方法，这里默认在以下几篇文章的CORDIC原理基础之上完成。

HLS / Chisel 实现CORDIC算法计算反三角函数（圆坐标系向量模式）
HLS / Chisel 实现CORDIC算法双曲系统
HLS / Chisel 利用CORDIC双曲系统实现平方根计算

原理

CORDIC背景

CORDIC算法有旋转和向量两个模式，分别可以在圆坐标系、线性坐标系和双曲线坐标系使用，从而可以演算出8种运算，而结合这8种运算也可以衍生出其他许多运算。下表展示了8种运算在CORDIC算法中计算的核心公式推导：

接下来我们将利用如上的CV、HV三种模式实现复数平方根

复数平方根

首先有如下的复数平方根公式$\sqrt{x+iy}$的结果推导，首先转换为极坐标模式$\sqrt{rcos\theta +irsin\theta }$，由De Moivre公式可以得到复数乘方的公式如下：

$$Z^n=[r(cos\theta +isin\theta {]}^n=r^n[cos(n\theta )+isin(n\theta )]$$

则求平方根为：

$$Z^{\frac{1}{2}}=\sqrt{r}(cos\frac{\theta }{2}+sin\frac{\theta }{2})$$

由于：

$$cos\frac{\theta }{2}=\pm \sqrt{\frac{1+cos\theta }{2}},sin\frac{\theta }{2}=\pm \sqrt{\frac{1-cos\theta }{2}}$$

且

$$x>0,\theta \in [-\frac{\pi }{2},\frac{\pi }{2}]$$

再将极坐标转化为原坐标系可以得到如下的推到结果：

$$\sqrt{x+iy}=\sqrt{\frac{\sqrt{x^2+y^2}+x}{2}}\pm i\sqrt{\frac{\sqrt{x^2+y^2}-x}{2}}$$

CORDIC复数平方根高性能计算原理

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

参考文章HLS / Chisel 实现CORDIC算法计算反三角函数（圆坐标系向量模式）

首先输入如下的数值

$${\mathrm{X}}_0=x,Y_0=y,Z_0=0$$

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

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

其中K在CV模式的推导文章中有介绍，是一个固定的常数。

2. HV_square_root模式计算

参考文章HLS / Chisel 利用CORDIC双曲系统实现平方根计算

将1中的结果 r 作如下处理：

$$a_{re}=(r+x)/2，a_{im}=(r-x)/2$$

即输入HV模块的值是：

$${\mathrm{X}}_0=a+1,Y_0=a-1,Z_0=0$$

具体HV模块内部的计算参考以上引用文章和HLS / Chisel 实现CORDIC算法双曲系统 实现.

用两个最终的HV_square_root模块分别计算实部虚部结果得到如下：

$$X=\sqrt{\frac{\sqrt{x^2+y^2}+x}{2}}Y=\sqrt{\frac{\sqrt{x^2+y^2}-x}{2}}$$

具体的流程模块如下：

HLS实践

#define COS_SIN_TYPE float
#define NUM_ITERATIONS 50

typedef struct { COS_SIN_TYPE re, im; } Complex; 

THETA_TYPE cordic_K = 0.6072529350088814

void cordic_complex_square_root(Complex x, Complex &res)
{
    COS_SIN_TYPE p,q;
    p = divisor.re;
    q = divisor.im;
 
    /* CV模式 */
    for(int j = 0;j < NUM_ITERATIONS; j++){
        COS_SIN_TYPE temp_cos = p;
        if(y < 0){
            p = p - q >> j;
            q = q + temp_cos >> j;
        } else{
            x = x + y >> j;
            q = q - temp_cos >> j;
        }
    }
    
    // Set ther final sine and cosine values
    r = p * cordic_K_n;
    
    /* 中间处理 */
    COS_SIN_TYPE  x_next = (r + p) >> 1;
    COS_SIN_TYPE  y_next = (r - p) >> 1;
    
    /* HV 模块，调用之前文章实现了的函数*/
    cordic_hv_square_root(x_next , res.re);
    cordic_hv_square_root(y_next , res.im);
}

CHISEL实践

定点数定义

首先定义定点数，注意在chisel.experimental包里有定点数的类

/* Cordic 的全局常数定义 */
trait Cordic_Const_Data extends HasDataConfig{
  /* 迭代次数配置 */
  val NUM_ITERATIONS = 20
}

trait HasDataConfig {
  /* 定点数的位宽定义 */ 
  val DataWidth = 32
  val BinaryPoint = 20
}

源码实现

class cordic_complex_square_root extends Module with HasDataConfig with Cordic_Const_Data {
  /*
   * @io.op : 输入的被除数 定点数表示的复数Complex 输入范围R+
   * @io.res : 输出的复数除法结果 定点数表示的复数Complex 范围R
   * details: 利用cordic实现复数除法，周期为2*迭代次数NUM_ITERATIONS
  **/
  val io = IO(new Bundle {
    val op: Complex = Input(new Complex())
    val res: Complex = Output(new Complex())
//    val debug_cv_out: FixedPoint = Output(FixedPoint(DataWidth.W, BinaryPoint.BP))
  })

  val x = WireInit(io.op.re)
  val y = WireInit(io.op.im)

  /* CV模式求square(x^2+y^2) */
  /* 初始化计算的寄存器数组，形成NUM_ITERATIONS级流水 */
  val current_x: Vec[FixedPoint] = RegInit(VecInit(Seq.fill(NUM_ITERATIONS)(FixedPoint.fromDouble(0.0, DataWidth.W, BinaryPoint.BP)))) // cos
  val current_y: Vec[FixedPoint] = RegInit(VecInit(Seq.fill(NUM_ITERATIONS)(FixedPoint.fromDouble(0.0, DataWidth.W, BinaryPoint.BP)))) // sin
  val cordic_K: FixedPoint = FixedPoint.fromDouble(0.6072529350088814, DataWidth.W, BinaryPoint.BP)
  /* 寄存最初输入的x值 */
  val reg_x: Vec[FixedPoint] = RegInit(VecInit(Seq.fill(NUM_ITERATIONS)(FixedPoint.fromDouble(0.0, DataWidth.W, BinaryPoint.BP)))) // cos
  for (i <- 0 until NUM_ITERATIONS) {
    /*
     * x[i] = K(x[i-1] - sigma * 2^(-i) * y[i-1])
     * y[i] = K(y[i-1] + sigma * 2^(-i) * x[i-1])
    **/
    if (i == 0) {
      reg_x(0) := x
      /* 流水线第一级直接对(io.x,io.y)点做运算 */
      when(y < FixedPoint.fromDouble(0.0, DataWidth.W, BinaryPoint.BP)) {
        current_x(i) := x - (y >> i)
        current_y(i) := y + (x >> i)
      }.otherwise {
        current_x(i) := x + (y >> i)
        current_y(i) := y - (x >> i)
      }
    } else {
      reg_x(i) := reg_x(i - 1)
      when(current_y(i - 1) < FixedPoint.fromDouble(0.0, DataWidth.W, BinaryPoint.BP)) {
        current_x(i) := current_x(i - 1) - (current_y(i - 1) >> i) // 移位替代乘法
        current_y(i) := current_y(i - 1) + (current_x(i - 1) >> i)
      }.otherwise {
        current_x(i) := current_x(i - 1) + (current_y(i - 1) >> i)
        current_y(i) := current_y(i - 1) - (current_x(i - 1) >> i)
      }
    }
  }

//  io.debug_cv_out := current_x(NUM_ITERATIONS - 1)*cordic_K
  val temp_cv_out = current_x(NUM_ITERATIONS - 1)*cordic_K

  /* HV模式求实数平方根 */
  val next_x = (temp_cv_out + reg_x(NUM_ITERATIONS - 1)) >> 1
  val next_y = (temp_cv_out - reg_x(NUM_ITERATIONS - 1)) >> 1
  io.res.re := cordic_hv_square_root(next_x)
  io.res.im := cordic_hv_square_root(next_y)
}

object cordic_complex_square_root{
  def apply(x: Complex): Complex = {
    /*
     * @io.op : 输入的被除数 定点数表示的复数Complex 输入范围R+
     * @io.res : 输出的复数除法结果 定点数表示的复数Complex 范围R
     * details: 利用cordic实现复数除法，周期为2*迭代次数NUM_ITERATIONS
    **/
    val unit = Module(new cordic_complex_square_root)
    unit.io.op := x
    unit.io.res
  }
}

测试

import org.scalatest._
import chisel3._
import chiseltest._
import chisel3.experimental._
import scala.math._

class Cordic_Complex_Square_Root_Tester extends FlatSpec with ChiselScalatestTester with Matchers {
  behavior of "mytest2"
  it should "do something" in {
    val NUM_ITERATIONS = 20
    test(new cordic_complex_square_root) { c =>
      val opre = 0.6
      val opim = 0.6
      c.io.op.re.poke(FixedPoint.fromDouble(opre, 32.W, 20.BP))
      c.io.op.im.poke(FixedPoint.fromDouble(opim, 32.W, 20.BP))
      // 对比组
//      val cv_out = pow(opre * opre + opim * opim, 0.5)
      val opposite_resre = pow((pow(opre * opre + opim * opim, 0.5) + opre) / 2, 0.5)
      val opposite_resim = pow((pow(opre * opre + opim * opim, 0.5) - opre) / 2, 0.5)

      for (i <- 0 to NUM_ITERATIONS * 2) {
        println(s"*****************${i}******************")
        println(s"res.re ${c.io.res.re.peek}")
        println(s"res.im ${c.io.res.im.peek}\n")
//        println(s"debug_cv_out ${c.io.debug_cv_out.peek}\n")
        c.clock.step(1)
      }

//      println(s"cv_out ${cv_out}")
      println(s"opposite_res.re ${opposite_resre}")
      println(s"opposite_res.im ${opposite_resim}")
    }

    test(new cordic_complex_square_root) { c =>
      val opre = 3.0
      val opim = -3.0
      c.io.op.re.poke(FixedPoint.fromDouble(opre, 32.W, 20.BP))
      c.io.op.im.poke(FixedPoint.fromDouble(opim, 32.W, 20.BP))
      // 对比组
      val opposite_resre = pow((pow(opre * opre + opim * opim, 0.5) + opre) / 2, 0.5)
      val opposite_resim = pow((pow(opre * opre + opim * opim, 0.5) - opre) / 2, 0.5)
      c.clock.step(41)

      println(s"res.re ${c.io.res.re.peek}")
      println(s"res.im ${c.io.res.im.peek}\n")

      println(s"opposite_res.re ${opposite_resre}")
      println(s"opposite_res.im ${opposite_resim}")
    }

  }
}

// 执行 sbt "testOnly Cordic_Complex_Square_Root_Tester"