FPGA之CORDIC算法：①CORDIC的基本原理介绍

只看该作者 · 2023-3-28 11:12

CORDIC的基本原理介绍

目前，学习与开发FPGA的程序员们大多使用的是Verilog HDL语言（以下简称为Verilog），关于Verilog的诸多优点一休哥就不多介绍了，在此，我们将重点放在Verilog的运算操作上。

我们都知道，在Verilog中，运算一般分为逻辑运算（与或非等）与算术运算（加减乘除等）。而在一开始学习Verilog时，老司机一定会提醒我们，“切记，千万别用‘/’除、‘%’取模（有的也叫取余）和‘**’幂。”这话说的不无道理，因为这三个运算是不可综合的。但，需清楚理解的是，不可综合的具体意思为不能综合为简单的模块，当我们在程序中调用了这些运算时，‘/’除和‘%’取模在Quartus软件中是可以综合的，因此可以正常调用运行，但是会消耗一些逻辑资源，而且会产生延时，即这两个运算的处理时间会很长，可能会大于时序控制时钟的单周期时间。此时呢，我们会建议你调用IP核来实现运算操作，虽然这样也会消耗许多逻辑资源，但产生的延时相对较小满足了你基本的需求。

问题好像迎刃而解了，可是仔细一想，除了这些运算，我们还剩下什么？对呀，三角函数，反三角函数，对数函数，指数函数呢，这些函数我们在高中就学习了的呀，难道在FPGA中就没有用武之地吗？有人会说，查找表呗，首先将某个运算的所有可能的输入与输出对一一罗列出来，然后放进Rom中，然后根据输入查表得到输出。这个方法虽然有效的避免了延时问题，却是一个十分消耗资源的方法，不适合资源紧张的设计。那么，就真的没有办法了吗？

答案就是咱们今天的标题了，CORDIC，而且CORDIC是一个比较全能的算法，通过这一原理，我们可以实现三角函数，反三角函数，对数函数，指数函数等多种运算。接下来，一休哥就带领大家来学习CORDIC的原理吧。（题外话：请相信一休哥，本文不会让你感到太多痛苦~）

本文将分三个小部分来展开介绍：
1、CORDIC的基本原理介绍
2、CORDIC的具体操作流程介绍
3、CORDIC的旋转模式——Verilog仿真

只看该作者 · 2023-3-28 11:13

一、CORDIC的基本原理介绍
CORDIC算法是一个“化繁为简”的算法，将许多复杂的运算转化为一种“仅需要移位和加法”的迭代操作。CORDIC算法有旋转和向量两个模式，分别可以在圆坐标系、线性坐标系和双曲线坐标系使用，从而可以演算出8种运算，而结合这8种运算也可以衍生出其他许多运算。下表展示了8种运算在CORDIC算法中实现的条件。

首先，我们先从旋转模式下的圆坐标系讲起，这也是CORDIC方法最初的用途。

1、CORDIC的几何原理介绍
假设在xy坐标系中有一个点P1（x1，y1），将P1点绕原点旋转θ角后得到点P2（x2，y2）。

于是可以得到P1和P2的关系。

x2= x1cosθ – y1sinθ = cosθ(x1 – y1tanθ)

y2= y1cosθ + x1sinθ = cosθ(y1 +x1tanθ)

以上就是CORDIC的几何原理部分，而我们该如何深入理解这个几何原理的真正含义呢？

从原理中，我们可以知道，当已知一个点P1的坐标，并已知该点P2旋转的角度θ，则可以根据上述公式求得目标点P2的坐标。然后，麻烦来了，我们需要用FPGA去执行上述运算操作，而FPGA的Verilog语言根本不支持三角函数运算。因此，我们需要对上述式子进行简化操作，将复杂的运算操作转换为一种单一的“迭代位移”算法。那么，接下来我们介绍优化算法部分。

只看该作者 · 2023-3-28 11:14

2、CORDIC的优化算法介绍
我们先介绍算法的优化原理：将旋转角θ细化成若干分固定大小的角度θi，并且规定θi满足tanθi = 2-i，因此∑θi的值在[-99.7°，99.7°]范围内，如果旋转角θ超出此范围，则运用简单的三角运算操作即可（加减π）。

然后我们需要修改几何原理部分的假设，假设在xy坐标系中有一个点P0（x0，y0），将P0点绕原点旋转θ角后得到点Pn（xn，yn）。

于是可以得到P0和Pn的关系。

xn= x0cosθ – y0sinθ = cosθ(x0 – y0tanθ)

yn= y0cosθ + x0sinθ = cosθ(y0 + x0tanθ)

然后，我们将旋转角θ细化成θi，由于每次的旋转角度θi是固定不变的（因为满足tanθi = 2-i），如果一直朝着一个方向旋转则∑θi一定会超过θ（如果θ在[-99.7°，99.7°]范围内）。因此我们需要对θi设定一个方向值di。如果旋转角已经大于θ，则di为-1，表示下次旋转为顺时针，即向θ靠近；如果旋转角已经小于θ，则di为1，表示下次旋转为逆时针，即也向θ靠近。然后我们可以得到每次旋转的角度值diθi，设角度剩余值为zi+1，则有zi+1 = zi - diθi，其中z0为θ。如此随着i的增大，角度剩余值zi+1 将会趋近于0，此时运算结束。（注：可以发现，di与zi的符号位相同）

第一次旋转θ0，d0为旋转方向：

x1 = cosθ0(x0 – d0y0tanθ0)

y1 = cosθ0(y0 + d0x0tanθ0)

第二次旋转θ1，d1为旋转方向：

x2 = cosθ1(x1 – d1y1tanθ1)

= cosθ1cosθ0(x0 – d0y0tanθ0 – d1y0tanθ1 – d1d0 x0tanθ1tanθ0)

y2 = cosθ1(y1 + d1x1tanθ1)

= cosθ1cosθ0(y0 + d0x0tanθ0 + d1x0tanθ1 – d1d0y0tanθ1tanθ0)

大家可能已经发现了，在我们旋转的过程中，式子里一直会有tanθi和cosθi，而每次都可以提取出cosθi。虽然我们的FPGA无法计算tanθi，但我们知道tanθi = 2-i，因此可以执行和tanθi效果相同的移位操作2-i来取代tanθi。而对于cosθi，我们可以事先全部提取出来，然后等待迭代结束之后（角度剩余值zi+1趋近于0，一般当系统设置最大迭代次数为16时zi+1已经很小了），然后计算出∏cosθi的值即可。

总结一下：

迭代公式有三：

xi+1 = xi – diyi2-i，提取了cosθi，2-i等效替换了tanθi之后

yi+1 = yi + dixi2-i，提取了cosθi，2-i等效替换了tanθi之后

zi+1 = zi - diθi

其中i从0开始迭代，假设当i = n-1时，zn趋近于0，迭代结束。然后对结果乘上∏cosθi（i从0至n-1），于是得到点Pn（xn∏cosθi，yn∏cosθi），此时的点Pn就近似等于之前假设中的点Pn（xn，yn）了，所以此时的点Pn同样满足之前假设得到的公式：

xn∏cosθi = x0cosθ – y0sinθ

yn∏cosθi = y0cosθ + x0sinθ

由于i从0至n-1，所以上式可以转化成下式：

xn = 1/∏cosθi (x0cosθ – y0sinθ)，（其中i从0至n-1）

yn = 1/∏cosθi (y0cosθ + x0sinθ)，（其中i从0至n-1）

注意：上式中的xn，yn是经过迭代后的结果，而不是之前假设中的点Pn（xn，yn）。了解这点是十分重要的。

根据高中学的三角函数关系，可以知道cosθi = 1/[(1+tan2θi)^0.5] = 1/[(1+2-2i)^0.5]，而1/[(1+2-2i)^0.5]的极值为1，因此我们可以得出一个结论：当i的次数很大时，∏cosθi 的值趋于一个常数。

关于如何计算∏cosθi 的代码如下所示：

复制

close all;



clear;



clc;



% 初始化



die=16;%迭代次数



jiao=zeros(die,1);%每次旋转的角度



cos_value=zeros(die,1);%每次旋转的角度的余弦值



K=zeros(die,1);%余弦值的N元乘积



K_1=zeros(die,1);%余弦值的N元乘积的倒数



fori=1:die



    a=2^(-(i-1));



    jiao(i)=atan(a);



    cos_value(i)=cos(jiao(i));



    if(i==1)



        K(i)=cos_value(i);



        K_1(i)=1/K(i);



    else



        K(i)=K(i-1)*cos_value(i);



        K_1(i)=1/K(i);



    end



end



jiao=vpa(rad2deg(jiao)*256,10)



cos_value=vpa(cos_value,10)



K=vpa(K,10)



K_1=vpa(K_1,10)

从上表也可以看出，当迭代次数为16，i=15时，cosθi 的值已经非常趋近于1了，∏cosθi 的值则约等于0.607253，1/∏cosθi为1.64676。所以当迭代次数等于16时，通过迭代得到的点Pn坐标已经非常接近之前假设中的点Pn。所以，当迭代次数等于16时，这个式子是成立的。

xn = 1/∏cosθi (x0cosθ – y0sinθ)，（其中i从0至n-1）

yn = 1/∏cosθi (y0cosθ + x0sinθ)，（其中i从0至n-1）

此时，已知条件有三个x0、y0和θ。通过16次迭代，我们可以得到xn和yn。而式中的∏cosθi 是个随i变化的值，我们可以预先将其值存入系统中。

然后，我们人为设置x0 = ∏cosθi ，y0 = 0，则根据等式，xn = cosθ，yn = sinθ。其中1/∏cosθi 的值我们也同样预先存入系统中。如此，我们就实现了正弦和余弦操作了。