CORDIC( Coordinate rotation digital algorithm ) It's a computational triangle 、 Numerical algorithms for hyperbolic and other mathematical functions , Each operation produces a result output . This enables us to adjust the accuracy of the algorithm according to the application requirements ; More accurate results can be obtained by increasing the number of operation iterations .CORDIC Is to use only addition 、 Subtraction 、 Shift and look-up table implementation of a simple algorithm , This algorithm is in FPGA High efficiency in , It is often used in hardware algorithm implementation .

The paper [1] This paper introduces a method through CORDIC In the algorithm, CR、CV、HV Three modes are used to split the method of high-performance calculation of complex square root , This article will use Chisel Implement the method , The default here is in the following articles CORDIC Based on the principle .

HLS / Chisel Realization CORDIC The algorithm calculates the inverse trigonometric function （ Vector mode of circular coordinate system ）
HLS / Chisel Realization CORDIC Algorithmic hyperbolic system
HLS / Chisel utilize CORDIC Calculation of square root of hyperbolic system

principle

CORDIC background

CORDIC The algorithm has two modes: rotation and vector , They can be in circular coordinates 、 Linear coordinate system and hyperbolic coordinate system use , So we can work out 8 Operations , And combined with this 8 This operation can also derive many other operations . The following table shows 8 Operation in CORDIC The derivation of the core formula in the algorithm ：

Next, we will use the above CV、HV Three modes realize the square root of complex numbers

The square root of a complex number

First of all, there is the following complex square root formula $\sqrt{x+iy}$ The results are derived , First convert to polar mode $\sqrt{rcos\theta +irsin\theta }$, from De Moivre The formula of complex power can be obtained as follows ：

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

Then find the square root as ：

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

because ：

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

And

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

Then the polar coordinates are transformed into the original coordinate system, and the following results can be obtained ：

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

CORDIC High performance computing principle of complex square root

1. CV The mode is converted to polar form

Reference article HLS / Chisel Realization CORDIC The algorithm calculates the inverse trigonometric function （ Vector mode of circular coordinate system ）

First enter the following values

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

Here, the inverse trigonometric function calculation module can get the following results ：

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

among K stay CV The derivation of the model is introduced in the article , Is a fixed constant .

2. HV_square_root Mode calculation

Reference article HLS / Chisel utilize CORDIC Calculation of square root of hyperbolic system

take 1 Results in r Do the following ：

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

The input HV The value of the module is ：

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

Specifically HV For the calculation inside the module, refer to the above cited articles and HLS / Chisel Realization CORDIC Algorithmic hyperbolic system Realization .

Use two final HV_square_root The module calculates the real part and imaginary part respectively, and the results are as follows ：

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

The specific process modules are as follows ：

HLS practice

#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 Pattern  */
    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;
    
    /*  Intermediate treatment  */
    COS_SIN_TYPE  x_next = (r + p) >> 1;
    COS_SIN_TYPE  y_next = (r - p) >> 1;
    
    /* HV  modular , Call the function implemented in the previous article */
    cordic_hv_square_root(x_next , res.re);
    cordic_hv_square_root(y_next , res.im);
}

CHISEL practice

Fixed point number definition

First define the number of fixed points , Pay attention to chisel.experimental There are fixed-point classes in the package

/* Cordic  Global constant definition of  */
trait Cordic_Const_Data extends HasDataConfig{
    
  /*  Number of iterations configuration  */
  val NUM_ITERATIONS = 20
}

trait HasDataConfig {
    
  /*  Bit width definition of fixed-point number  */ 
  val DataWidth = 32
  val BinaryPoint = 20
}

The source code to achieve

class cordic_complex_square_root extends Module with HasDataConfig with Cordic_Const_Data {
    
  /* * @io.op :  Entered divisor   A complex number represented by a fixed-point number Complex  Input range R+ * @io.res :  Output complex division result   A complex number represented by a fixed-point number Complex  Range R * details:  utilize cordic Realize complex division , The period is 2* The number of iterations 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 Pattern seeking square(x^2+y^2) */
  /*  Initializes the calculated register array , formation NUM_ITERATIONS It's a cascade of water  */
  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)
  /*  Deposit the originally entered x value  */
  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
      /*  The first stage of the pipeline is directly connected to (io.x,io.y) Click to do the operation  */
      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) //  Shift instead of multiplication 
        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 Find the square root of a real number in a pattern  */
  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 :  Entered divisor   A complex number represented by a fixed-point number Complex  Input range R+ * @io.res :  Output complex division result   A complex number represented by a fixed-point number Complex  Range R * details:  utilize cordic Realize complex division , The period is 2* The number of iterations NUM_ITERATIONS **/
    val unit = Module(new cordic_complex_square_root)
    unit.io.op := x
    unit.io.res
  }
}

test

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))
      //  The control group 
// 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))
      //  The control group 
      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}")
    }

  }
}

//  perform  sbt "testOnly Cordic_Complex_Square_Root_Tester"