定点CORDIC算法求所有三角函数及向量模的原理分析、硬件实现(FPGA)

一、CORDIC算法

  CORDIC(Coordinate Rotation DIgital Computer)是一种通过迭代对多种数学函数求值的方法,它可以对三角函数、双曲函数和平面旋转问题进行求解。

  在CORDIC之前,要对特殊函数求值,最自然的方法便是级数展开,例如利用泰勒展开来逼近目标函数,只要阶数取得足够大,就可以无限逼近目标函数。级数展开在数学上是完美的,但运用到计算机时,我们很快就会发现问题:级数展开本质是用多项式函数来近似目标函数,这其中包括大量复杂浮点运算,对于没有硬件浮点运算单元的平台,只能通过软件浮点实现,效率很低。

  CORDIC的出现解决了这个问题。该算法利用迭代逼近的方法,仅仅通过加/减和移位操作,即可求出特殊函数的值,极大的方便了计算机实现。

本文所做的工作:

  一、从CORDIC算法正向角度分析三角函数(sin, cos, tan)对任意角的求值;

  二、从CORDIC算法逆向角度分析反三角函数(arcsin、arccos、arctan)对任意角的求值、向量模的求值;

  三、分析定点运算的实现及软件模型的建立

  四、通过Verilog HDL设计硬件。

本文代码仓库详见:https://github.com/sci-dev-git/CORDIC-all-in-one-verilog

二、核心思想

  正如该算法的名字所说,CORDIC最初是为一种用来进行坐标轴旋转的专用计算机开发的(其原型硬件于1959年成功应用于实时导航)。既然只是坐标轴旋转,算法最早的出发点也许是解决旋转问题。追根溯源,早期研究者是如何想到这个算法的呢?它与函数求值究竟有什么关系?这正是接下来笔者尝试说明的问题。

  【1】、从算法的正向角度(对于确定的点,已知旋转角,求旋转后的坐标)分析:

  假设现在我们要将一个直角坐标系中的点P(x0, y0)绕原点逆时针旋转z0角度,则变换后的点P1(x1, y1)坐标如下:

定点CORDIC算法求所有三角函数及向量模的原理分析、硬件实现(FPGA)

  我们发现这个转换式中涉及三角函数,但现在假设机器还不能求任意三角函数值,那么能否改进?

  可以考虑查表实现,但查表法要求数据必须是离散的,这样旋转角只能取有限个值。如何对任意的旋转角求解?

  简单,利用迭代逼近,把目标角度分成若干个小角度,每次迭代只旋转一个特定角度并靠近目标,通过若干次迭代后,点就被旋转到了近似的位置上。关键点是,通过特定的分割,使得每次旋转的角度都是特定角,如果将些特定角的三角函数值固化到列表中,就可以利用查表绕开三角函数的计算。对计算机而言虽然需要迭代多次,但每次都是简单运算,总体速度快于复杂的三角求值。一般该算法的结果只是近似值,但通过设置迭代次数就可以控制精度,迭代次数越多,精度越高,并最终趋于稳定。

  

例如:将点P0旋转z = 50°到P‘,求P’坐标。

  首先说明:对所有问题来说,z的取值都是[0,90°],因此假设每次旋转的角度都是90°的n分之一,这样便于制作三角函数表。

  第一步正向旋转45°变换到P1,发现没有达到目标角度;因此又继续将P1正向旋转22.5°变换到P2,发现超出了目标角度;于是继续逆向旋转11.25°变换到P3,也超出了目标角度;再逆向旋转5.625°,到达P4,这时我们发现P4已经很趋近答案了,在该精度内可以用P4近似表示P’坐标。

定点CORDIC算法求所有三角函数及向量模的原理分析、硬件实现(FPGA)

  这个算法基本上解决了旋转变换中三角函数的问题,但并不完美:每次迭代都需要进行4次浮点乘法运算,因为旋转角z的正余弦都是复杂的浮点数,例如cos 22.5° ≈ 0.923879,浮点运算的开销比较明显。

  继续观察旋转变换式,可以变形成:

定点CORDIC算法求所有三角函数及向量模的原理分析、硬件实现(FPGA)

  上式仍然涉及两个三角函数,但是cos被放到了一边,于是重点研究tan。将角度z微分为n个小角度,如果像下面这样特殊选取z,使tan z恰好只与2的幂有关,则原本复杂的乘tan zn运算就可以通过右移n位实现(xn,yn都是整数,相关问题在后面讨论),这对二进制计算机是非常自然的。通过这种方法成功地消除了乘法。

定点CORDIC算法求所有三角函数及向量模的原理分析、硬件实现(FPGA)

  结合上面的变换式子,并把cosz隐藏到系数P中,我们就有了递推关系(其中n表示当前迭代次数,n+1则表示下一次迭代。这里只写出了逆时针旋转的情况,顺时针旋转改变1/2n项的符号即可。)

定点CORDIC算法求所有三角函数及向量模的原理分析、硬件实现(FPGA)

  由于n是离散的,并且z = atan(1/2n)可以提前计算,所以可用查表的方法快速得出z的值,这样系数Pn也可以通过查表求出,但系数Pn仍是浮点数。

  到这里,原先算法的4次浮点运算,被成功减少到了只有2次,我们离真正的CORDIC算法已经很近了。

  但问题还没有结束,如果迭代n次的话,仍需要进行2n次浮点运算,有没有优化的余地呢?引起2n次浮点运算的原因是每次都需要与系数Pn相乘,这样做真的有必要吗?

  构造辅助三角形,可以得出:

定点CORDIC算法求所有三角函数及向量模的原理分析、硬件实现(FPGA)

  因此Pn的取值只与n有关,而与别的变量没有关系。Pn完全可以在所有迭代都完成后单独计算,在最后将结果乘P=ΠPn即可。

  对于根据迭代较少的情况,可以将P的不同取值固化到表格中,通过查表快速求出。

定点CORDIC算法求所有三角函数及向量模的原理分析、硬件实现(FPGA)

  而对于迭代次数较多的情况,可以将P看作常数近似处理,这是因为随着n的增加,P逐渐稳定下来:

定点CORDIC算法求所有三角函数及向量模的原理分析、硬件实现(FPGA)

  通过数值统计,我们可以看出具体多少次迭代以后开始可以将P看作常量。从图中可以看出从n=8开始比较合适。当n<3时P的取值变化较大,但事实上小于3次的迭代基本上是没有意义的(求解特殊角除外)。

  定点CORDIC算法求所有三角函数及向量模的原理分析、硬件实现(FPGA)

  到这里整个算法只需在最后进行一次浮点运算,而每次迭代都只涉及简单的加/减和位移运算了!!!

  虽然有了这个快速旋转变换算法,但我们的最终目标并不在于此,而是求三角函数的值,如何通过旋转求三角函数?

  简单。借助单位圆这个工具,在直角坐标系中将点(0, 1)绕原点旋转至α角上,则旋转后的点的横纵坐标对应的正是cosα、sinα的值,tanα也可以利用比值求出,这正是利用CORDIC求三角函数的核心思路。下图显示了CORDIC算法的迭代过程。

定点CORDIC算法求所有三角函数及向量模的原理分析、硬件实现(FPGA)

  上述算法可用伪代码描述如下:

定点CORDIC算法求所有三角函数及向量模的原理分析、硬件实现(FPGA)

     ) Then
         X(n + ) := X(n) – (Yn >> n)
         Y(n + ) := Y(n) + (Xn >> n)

         Z(n + ) := Z(n) - atan(/^n)
     Else
         X(n + ) := X(n) + (Yn >> n)
         Y(n + ) := Y(n) – (Xn >> n)

         Z(n + ) := Z(n) + atan(/^n)
     End if
 End for

  【2】、从算法的逆向角度(对于确定的点,已知旋转后的点坐标,求旋转角)说起。

  从正向分析或者逆向分析,角度不同,解决的问题也不同。

  仍然在单位圆内,假设已知正弦或余弦中的一个值,如何求对应的反三角函数?这里以求反正弦函数为例说明,已知反正弦函数自变量为S,我们将点(1, 0)逆时针旋转,如果转到某个角度该点的纵坐标恰等于S,说明这个角度就是反正弦函数的值;同理反余弦函数也可用用类似的方法算出。

  接下来继续分析反正切的情况,已知正切值,则可以得出对应的单位圆上的点P0(x0, y0)。如果将点P0旋转一个角度,使旋转后的点纵坐标变成0,那么这个角度正是反正切函数的值

  以上便是利用CORDIC求反三角函数的所有思路。

  现在留下的问题是到底旋转多少角度才合适。CORDIC同样以迭代逼近的方法解决这个问题。

  以求反正切函数为例:这里只讨论z在0到90°范围内(即x0, y0都为正)的情况。我们可以先将原始点旋转atan(-1/2)角度进行试探,旋转完成后发现点纵坐标仍然大于0,于是继续旋转atan(-1/4)、atan(-1/8)到达图中草绿色终点,这个时候发现纵坐标小于0,说明旋转过头了,于是再反方向旋转atan(-1/16)角度到达蓝色终点,若此时点的纵坐标在规定精度内接近0则迭代结束,可以通过累加所有迭代中旋转的角度增量求出α的具体数值,这个总角度α正是atan的值。

定点CORDIC算法求所有三角函数及向量模的原理分析、硬件实现(FPGA)

  通过第一小节的讨论可以看出,递推式中系数P在几何上的效果是缩放了旋转半径(点到坐标原点的距离)。由于求反正切函数过程中只用到旋转角,而旋转角与旋转半径与无关,所以不需要考虑系数P

  上述算法可用伪代码描述如下:

定点CORDIC算法求所有三角函数及向量模的原理分析、硬件实现(FPGA)

 A := 

     ) Then
         X(n + ) := X(n) + (Yn >> n)
         Y(n + ) := Y(n) - (Xn >> n)

         A := A + atan(/^n)
     Else
         X(n + ) := X(n) - (Yn >> n)
         Y(n + ) := Y(n) + (Xn >> n)

         A := A - atan(/^n)
     End if
 End for

  理解了求反正切函数的思路后,我们就可以着手考虑反正弦和反余弦函数的求值。比较可知,旋转过程中,对反正切函数要求旋转点的纵坐标趋于0,而对反正弦函数要求旋转点纵坐标趋于特定值S(函数自变量的取值)。可以发现,前者所涉及的的旋转是后者的一种特殊情况,这样只需简单修改前者就可以得出适用后者的旋转算法。

定点CORDIC算法求所有三角函数及向量模的原理分析、硬件实现(FPGA)

 limP := 0.607253
 A :=
 B := S / limP

     If (Y(n) >= B) Then
         X(n + ) := X(n) + (Yn >> n)
         Y(n + ) := Y(n) - (Xn >> n)

         A := A + atan(/^n)
     Else
         X(n + ) := X(n) - (Yn >> n)
         Y(n + ) := Y(n) + (Xn >> n)

         A := A - atan(/^n)
     End if
 End for

  需要解释的是B := S / limP,这里为什么需要将S放大1/limP倍呢?原因很简单,由于在每次迭代时都忽略了系数Pn<1,这会导致Y(n)比未忽略系数时大出1/Pn倍。为了使Y(n)与S能够进行比较,需要将S同比例放大1/Pn倍,这样每次迭代又将引入浮点运算。直接近似Pn = limP虽然会引入细微的精度损失,但避免了大量浮点运算。

  同理,也可以得出旋转点横坐标趋于特定值的算法。

 limP := 0.607253
 A :=
 B := S / limP

     If (X(n) >= B) Then
         X(n + ) := X(n) - (Yn >> n)
         Y(n + ) := Y(n) + (Xn >> n)

         A := A - atan(/^n)
     Else
         X(n + ) := X(n) + (Yn >> n)
         Y(n + ) := Y(n) - (Xn >> n)

         A := A + atan(/^n)
     End if
 End for

  

  至此所有关于利用CORDIC求三角和反三角函数的分析就结束了。

  借鉴上面求反正切角函数的思路,我们也可以看出向量求模的运算方法。直角坐标系中,将目标向量V的起点平移到坐标原点,终点用坐标(x0, y0)表示。如果把这个点旋转到x轴(或y轴)上,则旋转后的点对应的横坐标(或纵坐标)正是该向量的模长|V|。

  该算法的核心与求反正切函数基本相同,这里不再赘述。

 

  除三角反三角函数和向量模外,CORDIC可以完成的运算还有很多,本文不再重点讨论。

三、定点运算实现及软件模型建立

  如果能用整数表示浮点数,则可以通过整数运算完成实数运算。在一定精度内,定点数与浮点数可以通过比例放缩互相转换。

  首先确定整数位宽,这里以16位为例,

  1、量化角度:在直角范围内,可以将正整数的取值区间均分为coeff = 2^14 / 90 = 182个单位,则一个单位对应角度的数量级为10^-2。任意给定浮点角度Z°,其对应的定点角度为floor(Z * coeff);

  2、量化坐标(x, y):将原始坐标放大coeff倍并取整即可得出定点坐标。但若坐标取值较大则存在溢出风险,需另选系数。

  上述两点措施完成了算法输入的定点量化;而对于算法的输出,只需将定点数缩小coeff倍即可得到最终的浮点结果。

  现在终于可以开始设计工作了。不过在硬件设计之前,先考虑软件模型是有必要的,因为软件语言的抽象程度高于硬件描述语言,可以更加紧凑地对算法进行描述,同时能快速地进行调试。

  软件模型主要涉及两个核心算法:旋转和反旋转。

旋转算法的C语言模型:

  https://github.com/sci-dev-git/CORDIC-all-in-one-verilog/blob/master/CORDIC-rotate-fixed-point.c

反旋转算法:

  https://github.com/sci-dev-git/CORDIC-all-in-one-verilog/blob/master/CORDIC-anti-rotate-fixed-point.c

四、FPGA实现

  本文到这里已经接近尾声了,通过铺垫,FPGA的实现自然水到渠成。这一小节主要解决的问题是硬件实现的相关细节,例如象限转换,流水化处理等。采用流水线的目的在于提高时钟频率,例如在DDS(直接数字频率合成)应用中,CORDIC算法可以代替采样表生成波形数据。为了实现较高的合成频率,流水化是有必要的。

  关于Verilog实现有几点需要说明的地方,首先利用verilog标准中规定的generate语句,可以实现任意深度流水线的综合;

  其次在此之前对算法的所有讨论都只涉及第Ⅰ象限,现在加入象限的处理。电路直接取相位输入高2位判断象限,并按函数值在4个象限的符号关系对输出进行处理。

  完整的Verilog模块如下:该模块通过端口phase_i输入相位。通过sin_o,cos_o输出三角函数值,通过err_o输出迭代引起的相位误差。从流水线为空开始,需等待PIPE_DEPTH+2个时钟周期才能得出结果;而流水线填满后,对于后续相位输入,模块都可以在2个时钟周期内更新输出。

 module cordic_dds # (
    ,               /* Data width */
    ,       /* Pipeline depth */
    'h4dba        /* P = 0.607253 * 2^15 */
 )
 (/*AUTOARG*/
    // Outputs
    sin_o, cos_o, err_o,
    // Inputs
    clk, phase_i
 );

    input             clk;
    :]    phase_i;       /* Phase */
    ]     sin_o, cos_o;  /* Function value output */
    ]     err_o;         /* Phase Error output */

    ] cos_r=, sin_o_r=;
    ] x[PIPE_DEPTH:];
    ] y[PIPE_DEPTH:];
    ] z[PIPE_DEPTH:];

    ] atan_rom[PIPE_DEPTH:];

    :] quadrant [PIPE_DEPTH:];

    integer i;
    initial begin
       ; i<=PIPE_DEPTH; i=i+) begin
          x[i] = ; y[i] = ; z[i] = ;
          quadrant[i] = 'b0;
       end
    end

    initial begin
       atan_rom[] <= ;
       atan_rom[] <= ;
       atan_rom[] <= ;
       atan_rom[] <= ;
       atan_rom[] <= ;
       atan_rom[] <= ;
       atan_rom[] <= ;
       atan_rom[] <= ;
       atan_rom[] <= ;
       atan_rom[] <= ;
       atan_rom[] <= ;
       atan_rom[] <= ;
       atan_rom[] <= ;
       atan_rom[] <= ;
    end

    // ================= //
    // Pipeline stages   //
    // ================= //
    always @ (posedge clk) begin // stage 0
       x[] <= {'b0, limP};
       y[] <= ;
       z[] <= {'b0, phase_i[DW-1-2:0]}; // control the phase_i to the range[0-Pi/2]
    end

    always @ (posedge clk) begin // stage 1
       x[] <= x[] - y[];
       y[] <= x[] + y[];
       z[] <= z[] - atan_rom[]; // reversal 45deg
    end

    generate
       genvar k;
       ; k<PIPE_DEPTH; k=k+) begin
          always @ (posedge clk) begin
             if (z[k][DW]) begin /* the diff is negative on clockwise */
                x[k+] <= x[k] + {{k{y[k][DW]}},y[k][DW:k]}; /* >> k */
                y[k+] <= y[k] - {{k{x[k][DW]}},x[k][DW:k]}; /* >> k */
                z[k+] <= z[k] + atan_rom[k];
             end else begin
                x[k+] <= x[k] - {{k{y[k][DW]}},y[k][DW:k]};
                y[k+] <= y[k] + {{k{x[k][DW]}},x[k][DW:k]};
                z[k+] <= z[k] - atan_rom[k];
             end
          end
       end
    endgenerate

    // ================= //
    // Count quadrant    //
    // ================= //
    always @ (posedge clk) begin
       quadrant[] <= phase_i[DW-:DW-];
    end
    generate
       genvar j;
       ; j<PIPE_DEPTH; j=j+) begin
          always @ (posedge clk) begin
             quadrant[j+] <= quadrant[j];
          end
       end
    endgenerate

    // ================= //
    // Adjust quadrant   //
    // ================= //
    always @ (posedge clk)
       case(quadrant[PIPE_DEPTH])
          'b00: begin
             cos_r <= x[PIPE_DEPTH]; /* cos */
             sin_o_r <= y[PIPE_DEPTH]; /* sin */
          end
          'b01: begin
             cos_r <= ~(y[PIPE_DEPTH]) + 'b1; /* -sin */
             sin_o_r <= x[PIPE_DEPTH]; /* cos */
          end
          'b10: begin
             cos_r <= ~(x[PIPE_DEPTH]) + 'b1; /* -cos */
             sin_o_r <= ~(y[PIPE_DEPTH]) + 'b1; /* -sin */
          end
          'b11: begin
             cos_r <= y[PIPE_DEPTH]; /* sin */
             sin_o_r <= ~(x[PIPE_DEPTH]) + 'b1; /* -cos */
          end
       endcase

    assign cos_o = cos_r;
    assign sin_o = sin_o_r;
    assign err_o = z[PIPE_DEPTH];

 endmodule

  通过仿真得出波形如下:

定点CORDIC算法求所有三角函数及向量模的原理分析、硬件实现(FPGA)

  本文仅对CORDIC算法进行了一些浅显的分析,基本覆盖到了算法的核心思路和简单应用,但从工业应用的角度出发,还有很多值得讨论的问题未在文章中分析。由于笔者时间有限,不能将本文做得非常全面,疏漏之处有待日后逐步完善,还望各位读者海涵。

上一篇:Android 利用二次贝塞尔曲线模仿购物车加入物品抛物线动画


下一篇:java调用matlab生成exe文件