【Day29】Cordic 演算法的实现

假设今天再做某种数位信号处理时，不小心用到了 arctan(y/x) 函数，那麽当然可以用泰勒展开得到多项式，化成一连串的乘法与加法运算，但是在这里其实有另一个更有效率的方法，那就是基於座标旋转式计算机（Coordinate Rotation Digital Computer, CORDIC）算法。

而 Cordic 算法也被运用在众多场合当中，例如 DSP、适应性滤波器、FFT(Fast Fourier Transform)、DCT(Discrete cosine Transform)、解调器及神经网络上。

step1 : 资料的预处理(要将不在 -90 度 ~ 90 度内的座标旋转至此角度内)
- 若 x 为正（座标在右半部）：直接将 x 载入 x 暂存器，y 载入 y 暂存器
- 若 x 为负，y 为正（座标在第二象限）：将座标 -x 载入 y 暂存器，y 载入 x 暂存器（意思是将座标向右旋转 90 度至第一象限）
- 若 x 为负，y 为负（座标在第三象限）将座标 x 载入 y 暂存器，-y 载入 x 暂存器（意思是将座标向左旋转 90 度至第四象限）
step2 :
- 如果 x 轴之上则将此座标向右旋转 45 度（要有变数来记录目前转了多少角度了）
  - x = x + y
  - y = y - x
- 否则向左旋转 45 度
  - x = x - y
  - y = y + x
step3 :
- 如果 x 轴之上则将此座标向右旋转 26 度
  - x = x + y/2
  - y = y - x/2
- 否则向左旋转 90 度。
  - x = x - y/2
  - y = y + x/2
step4 :
- 如果 x 轴之上则将此座标向右旋转 14 度
  - x = x + y/4
  - y = y - x/4
- 否则向左旋转 14 度
  - x = x - y/4
  - y = y + x/4
step4 :
- 依照 step2 ~ step4 的规律做到达到所需精度即可
step5 : 最後结果即是：
- x 为半径也就是开根号的（x2 + y2）（但会是 1.618 倍）
- 而记录角度变数的则是 arctan(y/x) 的最後输出
- y则是误差值

注：角度在下图中

图片出处

这边来写一个简易的 Cordic 电路

宣告状态

/*------parameter------*/
localparam IDLE    = 3'd0;
localparam LOAD    = 3'd1;
localparam FIRST   = 3'd2;
localparam SECOND  = 3'd3;
localparam THIRD   = 3'd4;
localparam FINISH  = 3'd5;

定义输入输出

module Cordic#
(
  parameter width = 7
)
(
  clkSys, 
  rst_n,
  en,
  x_in,
  y_in,
  radius,
  phase,
  error,
  ready
);
/*------ports declaration------*/
input                   clkSys;
input                   rst_n;
input                   en;
input  signed [width:0] x_in;
input  signed [width:0] y_in;
output        [width:0] radius;
output        [width:0] phase;
output        [width:0] error;
output                  ready;
reg                     ready;
reg    signed [width:0] radius;
reg    signed [width:0] phase;
reg    signed [width:0] error;

宣告变数

/*------variables------*/
reg    signed [width:0] x[3:0];
reg    signed [width:0] y[3:0];
reg    signed [width:0] z[3:0];
reg               [2:0] fstate;

状态逻辑

/*------fstate_machine_state------*/
always@(posedge clkSys or negedge rst_n)begin
  if(!rst_n)fstate <= 2'd0;
  else begin
    case(fstate)
      IDLE:begin
        if(en)fstate <= LOAD;
        else  fstate <= IDLE;
      end
      LOAD:    fstate <= FIRST;
      FIRST:   fstate <= SECOND;
      SECOND:  fstate <= THIRD;
      THIRD:   fstate <= FINISH;
      FINISH:  fstate <= IDLE;
      default :fstate <= IDLE;
    endcase
  end
end

输出逻辑

/*------fstate_machine_output------*/
always@(posedge clkSys or negedge rst_n)begin
  if(!rst_n)begin
    x[0]<=0;x[1]<=0;x[2]<=0;x[3]<=0;
    y[0]<=0;y[1]<=0;y[2]<=0;y[3]<=0;
    z[0]<=0;z[1]<=0;z[2]<=0;z[3]<=0;
    ready  <= 1'b0;
    radius <= 0;
    phase  <= 0;
    error  <= 0;
  end
  else begin
    case(fstate)
      IDLE:begin
        x[0]<=0;x[1]<=0;x[2]<=0;x[3]<=0;
        y[0]<=0;y[1]<=0;y[2]<=0;y[3]<=0;
        z[0]<=0;z[1]<=0;z[2]<=0;z[3]<=0;
        ready <= 1'b0;
      end
      LOAD:begin
        if(x_in >= 0)begin
          x[0] <= x_in;
          y[0] <= y_in;
          z[0] <= 8'd0;
        end
        else if(y_in >= 0)begin
          x[0] <= y_in;
          y[0] <= -x_in;
          z[0] <= 8'd90;
        end
        else begin
          x[0] <= -y_in;
          y[0] <= x_in;
          z[0] <= -8'd90;
        end
      end
      FIRST:begin
        if(y[0] >= 0)begin
          x[1] <= x[0] + y[0];
          y[1] <= y[0] - x[0];
          z[1] <= z[0] + 8'd45;
        end
        else begin
          x[1] <= x[0] - y[0];
          y[1] <= y[0] + x[0];
          z[1] <= z[0] - 8'd45;
        end
      end
      SECOND:begin
        if(y[1] >= 0)begin
          x[2] <= x[1] + ({y[1][width], y[1][width:1]});//y[1]>>>1
          y[2] <= y[1] - ({x[1][width], x[1][width:1]});//x[1]>>>1
          z[2] <= z[1] + 8'd26;
        end
        else begin
          x[2] <= x[1] - ({y[1][width], y[1][width:1]});
          y[2] <= y[1] + ({x[1][width], x[1][width:1]});
          z[2] <= z[1] - 8'd26;
        end
      end
      THIRD:begin
        if(y[2] >= 0)begin
          x[3] <= x[2] + ({{2{y[2][width]}}, y[2][width:2]});//y[2]>>>2
          y[3] <= y[2] - ({{2{x[2][width]}}, x[2][width:2]});//x[2]>>>2
          z[3] <= z[2] + 8'd14;
        end
        else begin
          x[3] <= x[2] - ({{2{y[2][width]}}, y[2][width:2]});
          y[3] <= y[2] + ({{2{x[2][width]}}, x[2][width:2]});
          z[3] <= z[2] - 8'd14;
        end
      end
      FINISH:begin
        radius <= x[3];
        phase  <= z[3];
        error  <= y[3];
        ready  <= 1'b1;
      end
      default:begin
        x[0]<=0;x[1]<=0;x[2]<=0;x[3]<=0;
        y[0]<=0;y[1]<=0;y[2]<=0;y[3]<=0;
        z[0]<=0;z[1]<=0;z[2]<=0;z[3]<=0;
        ready  <= 1'b0;
        radius <= 0;
        phase  <= 0;
        error  <= 0;
      end
    endcase
  end
end

endmodule

TestBench


`timescale 1ns/1ns
module tb_Cordic();
localparam width = 7;
reg                    clkSys;
reg                    rst_n;
reg                    en;
reg   signed [width:0] x_in;
reg   signed [width:0] y_in;
wire  signed [width:0] radius;
wire  signed [width:0] phase;
wire  signed [width:0] error;
wire                   ready;
Cordic UUT (
  .clkSys(clkSys), 
  .rst_n(rst_n),
  .en(en),
  .x_in(x_in),
  .y_in(y_in),
  .radius(radius),
  .phase(phase),
  .error(error),
  .ready(ready)
);

always #5 clkSys = ~clkSys;

initial begin
  clkSys = 0;
  rst_n = 0;
  x_in = 0;
  y_in = 0;
  en = 0;
  repeat(5)@(posedge clkSys)rst_n = 0;
  rst_n = 1;
  test(-41,55);
  test(30, 40);
  #1000 $stop;
end
/*----------task---------*/
task test;
input signed [width:0] in1;
input signed [width:0] in2;
  begin
    #30;
    x_in = in1;
    y_in = in2;
    #10 en = 1;
    #10 en = 0;
    wait (ready == 1);
  end
endtask

endmodule

Wave