开发语言：Python，Verilog

开发环境：Pycharm，iVerilog+GTKWave

CORDIC，是一种迭代逼近的方法，仅通过加法、减法和移位运算，就能计算多种超越函数（如三角函数、反三角函数、双曲函数、平方根、对数、指数等）。

它在没有硬件乘法器/除法器的嵌入式系统中尤为重要（如FPGA、单片机）。

一、算法解析

假设在一个直角坐标系中，有一个半径为1，以坐标系原点O为圆心的单位圆；

假设圆上的某一点P(xp，yp)，P点与x轴正方向的夹角为α；

假设圆上的另一个点Q(xq，yq)，向量OQ与向量OP之间的夹角为θ；

由三角函数关系可知，

        xp = cosα，          yp = sinα；

        xq = cos(α+θ)，   yq = sin(α+θ)；

三角函数解算可得，

        xq = cos(α+θ) = cosα cosθ - sinα sinθ = cosθ * xp - sinθ * yp = cosθ(xp - yp tanθ)；

        yq = sin(α+θ) = sinα cosθ + sinθ cosα = cosθ yp + sinθ xp = cosθ(yp + xp tanθ)；

此时就可以得到经由OP向量旋转θ角后，得到的Q点的坐标，且Q点也位于单位圆上；

CORDIC算法的核心就是为了将sin，cos转换为加法，减法，乘法，但是由上式可知，公式中仍然存在着三角函数运算，因此还需要做一下变换。

将需要计算的θ角，分割为多次θi角的旋转，来迭代逐次逼近目标角度，同时为了满足在FPGA硬件层面计算的方便，设tanθi = di *，其中di表示此次旋转的方向，取值为+1 / -1。此外，还需要引入一个变量Z来表示当前旋转后得到的角度与目标角度之间的差值，如果本次旋转后的角度大于目标角度，那么下一次则需要向反方向旋转；如果本次旋转后的角度仍然小于目标角度，那么下一次则继续向原方向旋转；

也就是说，

第0次旋转的角度为：arctan(d0 * 1)， Z1 = Z0 - arctan(d0 * 1)；

                                   若Z1>0，则d1=1；若Z1<0，则d1 = -1；

第1次旋转的角度为：arctan(d1 * 1/2), Z2 = Z1 - arctan(d1 * 1/2);

                                   若Z2>0，则d2=1；若Z2<0，则d2 = -1；

第2次旋转的角度为：arctan(d2 * 1/4), Z3 = Z2 - arctan(d2 * 1/4)；

                                   若Z3>0，则d3=1；若Z3<0，则d3 = -1；

......

(Z0可设置为目标角度值)

以此类推，在经过多次迭代后，旋转角度累积，会逐渐逼近目标角度

根据以上所述旋转计算方法，可以得出旋转后得到新的点的坐标为

x(i+1) = cosθi(x(i) - y(i) tanθi) = cosθi(x(i) - y(i) * di * )

y(i+1) = cosθi(y(i) + x(i) tanθi) = cosθi(y(i) + x(i) * di * )

其中cosθi可以暂时先把它看做一个常数，它只会影响(x(i+1), y(i+1))的模，不会影响其角度。

到了这里基本上已经可以通过Python来实现一个软件版本的Cordic算法了。

 

二、Python实现Cordic算法

# math库
import math


# x0 : 初始点的x坐标
# y0 : 初始点的y坐标
# th : 待计算的角, 注意是角度值， 30°， 45°.....
# num: 迭代次数, 一般迭代到16次即可
def cordic(x0, y0, th, num):
    xi = x0;
    yi = y0;
    di = 1;
    r_sin = 0;
    r_cos = 0;
    # 先根据角度值来计算其对应sin值，和cos值的符号位，并将角度转换到第一象限
    if th >= 0 and th <=90:
        r_sin = 1;
        r_cos = 1;
        Z0 = th/180 * math.pi;
    elif th > 90 and th <= 180:
        r_sin = 1;
        r_cos = -1;
        Z0 = (180-th)/180 * math.pi;
    elif th > 180 and th <= 270:
        r_sin = -1;
        r_cos = -1;
        Z0 = (th - 180)/180 * math.pi;
    elif th > 270 and th <= 360:
        r_sin = -1;
        r_cos = 1;
        Z0 = (360-th)/180*math.pi;


    zi = Z0;
    # 先不管公式中的cosθi, 后续在FPGA实现中，这部分会优化
    for i in range(0, num):
        #print("xi=",xi, "yi=",yi,"zi=",zi, "di=",di);
        xn = xi - yi * di * 2**(-i);
        yn = yi + xi * di * 2**(-i);
        xi = xn;
        yi = yn;
        #print("atan=", math.atan(2**(-i)));
        zi = zi - di * math.atan(2**(-i));
        if zi > 0 :
            di = 1;
        elif zi < 0:
            di = -1;
        else:
            break;

    # 软件暂时将cosθi乘上，FPGA版本再优化
    for i in range(0, num):
        xi = xi * math.cos(math.atan(2**(-i)));
        yi = yi * math.cos(math.atan(2**(-i)));

    xi = xi * r_sin;
    yi = yi * r_cos;
    # xi即是sin
    # yi即是cos
    # yi/xi即是tan
    return xi, yi, yi/xi;



print(cordic(1, 0, 300, 16))

三、下期预告，Cordic算法的FPGA实现，以及仿真