CORDIC 算法实现 _FPGA

news2024/12/20 20:39:19

注：本文为 “CORDIC 算法” 相关文章合辑。

未整理去重。

如有内容异常，请看原文。

Cordic 算法的原理介绍

乐富道 2014-01-28 23:05

Cordic 算法知道正弦和余弦值，求反正切，即角度。

采用用不断的旋转求出对应的正弦余弦值，是一种近似求解发。

旋转的角度很讲求，每次旋转的角度必须使得 正切值近似等于 1 2 N \frac{1}{2^N} 2N1​。旋转的目的是让 Y Y Y 轴趋近与 0 0 0。把每次旋转的角度累加，即得到旋转的角度和即为正切值。

比如 Y Y Y 轴旋转 4 5 ∘ 45^{\circ} 45∘，则值减小 1 2 \frac{1}{2} 21​;

再旋转 26.5650 5 ∘ 26.56505^{\circ} 26.56505∘，再减少 1 4 \frac{1}{4} 41​;

再旋转角度 14.0362 4 ∘ 14.03624^{\circ} 14.03624∘，再减少 1 8 \frac{1}{8} 81​; 依次减少 1 16 \frac{1}{16} 161​， 1 32 ⋯ ⋯ \frac{1}{32}\cdots\cdots 321​⋯⋯，最后 Y Y Y 轴的值无限小，趋近于 0 0 0 。

比如 X = 1 X = 1 X=1， Y = 1 Y = 1 Y=1，的角度，角度是 4 5 ∘ 45^{\circ} 45∘。经过一次旋转，要使得 Y = 0 Y = 0 Y=0，这个角度必须是 4 5 ∘ 45^{\circ} 45∘。

img

如上图

img

如图中，直角坐标系中点（ X 0 X_0 X0​， Y 0 Y_0 Y0​）逆时钟旋转角度 θ \theta θ，变换成坐标（ X 1 X_1 X1​， Y 1 Y_1 Y1​）,那么用 X 0 X_0 X0​， Y 0 Y_0 Y0​,以及 θ \theta θ 的三角函数，如果表示 X 1 X_1 X1​， Y 1 Y_1 Y1​ 呢？

img

请想象，如果坐标也旋转角度 θ \theta θ，那么 X 1 X_1 X1​， Y 1 Y_1 Y1​ 的坐标依然是( X 0 X_0 X0​， Y 0 Y_0 Y0​)。接着往下看：

img

看完以上这副图，就该明白这个等式了：

x 1 = x 0 cos ⁡ θ − y 0 sin ⁡ θ x_1 = x_0\cos\theta - y_0\sin\theta x1​=x0​cosθ−y0​sinθ

y 1 = x 0 sin ⁡ θ + y 0 cos ⁡ θ y_1 = x_0\sin\theta + y_0\cos\theta y1​=x0​sinθ+y0​cosθ

再把这个式子化成正切函数。

img

Cordic 算法的思想是通过迭代的方法，不断的旋转特定的角度(这个特定的角度就是使得 Y Y Y 为上次的 1 2 \frac{1}{2} 21​)，使得累计旋转的角度的和无限接近某一设定的角度，

每次旋转的角度的 θ = arctan ⁡ ( 1 2 n ) \theta = \arctan(\frac{1}{2^n}) θ=arctan(2n1​)；

img

具体迭代如下表： Z 0 = 3 0 ∘ Z_0 = 30^{\circ} Z0​=30∘， Y 0 = 0 Y_0 = 0 Y0​=0， X 0 = 0.6073 X_0 = 0.6073 X0​=0.6073

输入 3 0 ∘ 30^{\circ} 30∘，经过 9 9 9 次迭代后， Z 0 = 0 Z_0 = 0 Z0​=0， Y 0 = 0.5006 Y_0 = 0.5006 Y0​=0.5006， X 0 = 0.8657 X_0 = 0.8657 X0​=0.8657

img

x ( i + 1 ) ′ = ( x i ′ − y i ′ ( σ i ) 2 − i ) x'_{(i + 1)} = (x'_i - y'_i(\sigma_i)2^{-i}) x(i+1)′​=(xi′​−yi′​(σi​)2−i)

y ( i + 1 ) ′ = ( x i ′ ( σ i ) 2 − i + y i ′ ) y'_{(i + 1)} = (x'_i(\sigma_i)2^{-i} + y'_i) y(i+1)′​=(xi′​(σi​)2−i+yi′​)

(当 i = 0 i = 0 i=0)

x 1 ′ = 0.607 − 0 ⋅ ( + 1 ) ⋅ 1 = 0.607 x'_1 = 0.607 - 0 \cdot (+1 ) \cdot 1 = 0.607 x1′​=0.607−0⋅(+1)⋅1=0.607

y 1 ′ = 0.607 ⋅ ( + 1 ) ⋅ 1 + 0 = 0.607 y'_1 = 0.607 \cdot (+1 ) \cdot 1 + 0 = 0.607 y1′​=0.607⋅(+1)⋅1+0=0.607

通过 Cordic 算法后，得到

y 9 = 0.5006 ( = sin ⁡ ( 3 0 ∘ ) ) y_9 = 0.5006 (=\sin(30^{\circ})) y9​=0.5006(=sin(30∘))

x 9 = 0.8657 ( = cos ⁡ ( 3 0 ∘ ) ) x_9 = 0.8657 (=\cos(30^{\circ})) x9​=0.8657(=cos(30∘))

所以也可以用 cordic 算法求出正切值的。

或者求反正切值：

img

计算公式：

img

posted @ 2014-01-28 23:05 乐富道

基于 FPGA 的 CORDIC 算法实现 —— Verilog 版

善良的一休君于 2017-08-21 20:07:34 发布

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

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

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

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

本文将分三个小部分来展开介绍：

1、CORDIC 的基本原理介绍

2、CORDIC 的具体操作流程介绍

3、CORDIC 的旋转模式 ——Verilog 仿真

本文涉及到的全部资料链接：

链接：http://pan.baidu.com/s/1gfrJzMj
密码：x92u

一、CORDIC 的基本原理介绍

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

这里写图片描述

这里写图片描述

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

1、CORDIC 的几何原理介绍

假设在 x y xy xy 坐标系中有一个点 P 1 ( x 1 , y 1 ) P1(x_1,y_1) P1(x1​,y1​)，将 P 1 P1 P1 点绕原点旋转 θ \theta θ 角后得到点 P 2 ( x 2 , y 2 ) P2(x_2,y_2) P2(x2​,y2​)。

这里写图片描述

这里写图片描述

于是可以得到 P 1 P1 P1 和 P 2 P2 P2 的关系：

x 2 = x 1 cos ⁡ θ − y 1 sin ⁡ θ = cos ⁡ θ ( x 1 − y 1 tan ⁡ θ ) y 2 = y 1 cos ⁡ θ + x 1 sin ⁡ θ = cos ⁡ θ ( y 1 + x 1 tan ⁡ θ ) \begin{align*} x_2 &= x_1\cos\theta - y_1\sin\theta = \cos\theta(x_1 - y_1\tan\theta)\\ y_2 &= y_1\cos\theta + x_1\sin\theta = \cos\theta(y_1 + x_1\tan\theta) \end{align*} x2​y2​​=x1​cosθ−y1​sinθ=cosθ(x1​−y1​tanθ)=y1​cosθ+x1​sinθ=cosθ(y1​+x1​tanθ)​

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

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

2、CORDIC的优化算法介绍

我们先介绍算法的优化原理：将旋转角 θ \theta θ 细化成若干分固定大小的角度 θ i \theta_i θi​，并且规定 θ i \theta_i θi​ 满足 tan ⁡ θ i = 2 − i \tan\theta_i = 2^{-i} tanθi​=2−i，因此 ∑ θ i \sum\theta_i ∑θi​ 的值在 [ − 99. 7 ∘ , 99. 7 ∘ ] [-99.7^{\circ},99.7^{\circ}] [−99.7∘,99.7∘] 范围内，如果旋转角 θ \theta θ 超出此范围，则运用简单的三角运算操作即可（加减 π \pi π）。

然后我们需要修改几何原理部分的假设，假设在 x y xy xy 坐标系中有一个点 P 0 ( x 0 , y 0 ) P0(x_0,y_0) P0(x0​,y0​)，将 P 0 P0 P0 点绕原点旋转 θ \theta θ 角后得到点 P n ( x n , y n ) Pn(x_n,y_n) Pn(xn​,yn​)。

于是可以得到 P 0 P0 P0 和 P n Pn Pn 的关系：

x n = x 0 cos ⁡ θ − y 0 sin ⁡ θ = cos ⁡ θ ( x 0 − y 0 tan ⁡ θ ) y n = y 0 cos ⁡ θ + x 0 sin ⁡ θ = cos ⁡ θ ( y 0 + x 0 tan ⁡ θ ) \begin{align*} x_n &= x_0\cos\theta - y_0\sin\theta = \cos\theta(x_0 - y_0\tan\theta)\\ y_n &= y_0\cos\theta + x_0\sin\theta = \cos\theta(y_0 + x_0\tan\theta) \end{align*} xn​yn​​=x0​cosθ−y0​sinθ=cosθ(x0​−y0​tanθ)=y0​cosθ+x0​sinθ=cosθ(y0​+x0​tanθ)​

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

第一次旋转 θ 0 \theta_0 θ0​， d 0 d_0 d0​ 为旋转方向：

x 1 = cos ⁡ θ 0 ( x 0 − d 0 y 0 tan ⁡ θ 0 ) y 1 = cos ⁡ θ 0 ( y 0 + d 0 x 0 tan ⁡ θ 0 ) \begin{align*} x_1 &= \cos\theta_0(x_0 - d_0y_0\tan\theta_0)\\ y_1 &= \cos\theta_0(y_0 + d_0x_0\tan\theta_0) \end{align*} x1​y1​​=cosθ0​(x0​−d0​y0​tanθ0​)=cosθ0​(y0​+d0​x0​tanθ0​)​

第二次旋转 θ 1 \theta_1 θ1​， d 1 d_1 d1​ 为旋转方向：

x 2 = cos ⁡ θ 1 ( x 1 − d 1 y 1 tan ⁡ θ 1 ) = cos ⁡ θ 1 cos ⁡ θ 0 ( x 0 − d 0 y 0 tan ⁡ θ 0 − d 1 y 0 tan ⁡ θ 1 − d 1 d 0 x 0 tan ⁡ θ 1 tan ⁡ θ 0 ) y 2 = cos ⁡ θ 1 ( y 1 + d 1 x 1 tan ⁡ θ 1 ) = cos ⁡ θ 1 cos ⁡ θ 0 ( y 0 + d 0 x 0 tan ⁡ θ 0 + d 1 x 0 tan ⁡ θ 1 − d 1 d 0 y 0 tan ⁡ θ 1 tan ⁡ θ 0 ) \begin{align*} x_2 &= \cos\theta_1(x_1 - d_1y_1\tan\theta_1)\\ &= \cos\theta_1\cos\theta_0(x_0 - d_0y_0\tan\theta_0 - d_1y_0\tan\theta_1 - d_1d_0x_0\tan\theta_1\tan\theta_0)\\ y_2 &= \cos\theta_1(y_1 + d_1x_1\tan\theta_1)\\ &= \cos\theta_1\cos\theta_0(y_0 + d_0x_0\tan\theta_0 + d_1x_0\tan\theta_1 - d_1d_0y_0\tan\theta_1\tan\theta_0) \end{align*} x2​y2​​=cosθ1​(x1​−d1​y1​tanθ1​)=cosθ1​cosθ0​(x0​−d0​y0​tanθ0​−d1​y0​tanθ1​−d1​d0​x0​tanθ1​tanθ0​)=cosθ1​(y1​+d1​x1​tanθ1​)=cosθ1​cosθ0​(y0​+d0​x0​tanθ0​+d1​x0​tanθ1​−d1​d0​y0​tanθ1​tanθ0​)​

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

总结一下：

迭代公式有三：

x i + 1 = x i − d i y i 2 − i (提取了 cos ⁡ θ i , 2 − i 等效替换了 tan ⁡ θ i 之后) y i + 1 = y i + d i x i 2 − i (提取了 cos ⁡ θ i , 2 − i 等效替换了 tan ⁡ θ i 之后) z i + 1 = z i − d i θ i \begin{align*} x_{i + 1} &= x_i - d_iy_i2^{-i} \quad \text{(提取了}\cos\theta_i, 2^{-i}\text{等效替换了}\tan\theta_i\text{之后)}\\ y_{i + 1} &= y_i + d_ix_i2^{-i} \quad \text{(提取了}\cos\theta_i, 2^{-i}\text{等效替换了}\tan\theta_i\text{之后)}\\ z_{i + 1} &= z_i - d_i\theta_i \end{align*} xi+1​yi+1​zi+1​​=xi​−di​yi​2−i(提取了cosθi​,2−i等效替换了tanθi​之后)=yi​+di​xi​2−i(提取了cosθi​,2−i等效替换了tanθi​之后)=zi​−di​θi​​

其中 i i i 从 0 0 0 开始迭代，假设当 i = n − 1 i = n - 1 i=n−1 时， z n z_n zn​ 趋近于 0 0 0，迭代结束。然后对结果乘上 ∏ i = 0 n − 1 cos ⁡ θ i \prod_{i = 0}^{n - 1}\cos\theta_i ∏i=0n−1​cosθi​，于是得到点 P n ( x n ∏ i = 0 n − 1 cos ⁡ θ i , y n ∏ i = 0 n − 1 cos ⁡ θ i ) Pn(x_n\prod_{i = 0}^{n - 1}\cos\theta_i,y_n\prod_{i = 0}^{n - 1}\cos\theta_i) Pn(xn​∏i=0n−1​cosθi​,yn​∏i=0n−1​cosθi​)，此时的点 P n Pn Pn 就近似等于之前假设中的点 P n ( x n , y n ) Pn(x_n,y_n) Pn(xn​,yn​) 了，所以此时的点 P n Pn Pn 同样满足之前假设得到的公式：

x n ∏ i = 0 n − 1 cos ⁡ θ i = x 0 cos ⁡ θ − y 0 sin ⁡ θ y n ∏ i = 0 n − 1 cos ⁡ θ i = y 0 cos ⁡ θ + x 0 sin ⁡ θ \begin{align*} x_n\prod_{i = 0}^{n - 1}\cos\theta_i &= x_0\cos\theta - y_0\sin\theta\\ y_n\prod_{i = 0}^{n - 1}\cos\theta_i &= y_0\cos\theta + x_0\sin\theta \end{align*} xn​i=0∏n−1​cosθi​yn​i=0∏n−1​cosθi​​=x0​cosθ−y0​sinθ=y0​cosθ+x0​sinθ​

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

x n = 1 ∏ i = 0 n − 1 cos ⁡ θ i ( x 0 cos ⁡ θ − y 0 sin ⁡ θ ) (其中 i 从 0 至 n − 1 ) y n = 1 ∏ i = 0 n − 1 cos ⁡ θ i ( y 0 cos ⁡ θ + x 0 sin ⁡ θ ) (其中 i 从 0 至 n − 1 ) \begin{align*} x_n &= \frac{1}{\prod_{i = 0}^{n - 1}\cos\theta_i}(x_0\cos\theta - y_0\sin\theta) \quad \text{(其中}i\text{从}0\text{至}n - 1\text{)}\\ y_n &= \frac{1}{\prod_{i = 0}^{n - 1}\cos\theta_i}(y_0\cos\theta + x_0\sin\theta) \quad \text{(其中}i\text{从}0\text{至}n - 1\text{)} \end{align*} xn​yn​​=∏i=0n−1​cosθi​1​(x0​cosθ−y0​sinθ)(其中i从0至n−1)=∏i=0n−1​cosθi​1​(y0​cosθ+x0​sinθ)(其中i从0至n−1)​

注意：上式中的 x n x_n xn​， y n y_n yn​ 是经过迭代后的结果，而不是之前假设中的点 P n ( x n , y n ) Pn(x_n,y_n) Pn(xn​,yn​)。了解这点是十分重要的。

根据高中学的三角函数关系，可以知道 cos ⁡ θ i = 1 ( 1 + tan ⁡ 2 θ i ) 1 2 = 1 ( 1 + 2 − 2 i ) 1 2 \cos\theta_i = \frac{1}{(1 + \tan^{2}\theta_i)^{\frac{1}{2}}} = \frac{1}{(1 + 2^{-2i})^{\frac{1}{2}}} cosθi​=(1+tan2θi​)21​1​=(1+2−2i)21​1​，而 1 ( 1 + 2 − 2 i ) 1 2 \frac{1}{(1 + 2^{-2i})^{\frac{1}{2}}} (1+2−2i)21​1​ 的极值为 1 1 1，因此我们可以得出一个结论：当 i i i 的次数很大时， ∏ i = 0 n − 1 cos ⁡ θ i \prod_{i = 0}^{n - 1}\cos\theta_i ∏i=0n−1​cosθi​ 的值趋于一个常数。

关于如何计算 ∏ cos ⁡ θ i \prod\cos\theta_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元乘积的倒数
for i = 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 16 16， i = 15 i = 15 i=15 时， cos ⁡ θ i \cos\theta_i cosθi​ 的值已经非常趋近于 1 1 1 了， ∏ i = 0 n − 1 cos ⁡ θ i \prod_{i = 0}^{n - 1}\cos\theta_i ∏i=0n−1​cosθi​ 的值则约等于 0.607253 0.607253 0.607253， 1 ∏ i = 0 n − 1 cos ⁡ θ i \frac{1}{\prod_{i = 0}^{n - 1}\cos\theta_i} ∏i=0n−1​cosθi​1​ 为 1.64676 1.64676 1.64676。所以当迭代次数等于 16 16 16 时，通过迭代得到的点 P n Pn Pn 坐标已经非常接近之前假设中的点 P n Pn Pn。所以，当迭代次数等于 16 16 16 时，这个式子是成立的。

x n = 1 ∏ i = 0 n − 1 cos ⁡ θ i ( x 0 cos ⁡ θ − y 0 sin ⁡ θ ) ( 其中  i  从  0  至  n − 1 ) y n = 1 ∏ i = 0 n − 1 cos ⁡ θ i ( y 0 cos ⁡ θ + x 0 sin ⁡ θ ) ( 其中  i  从  0  至  n − 1 ) \begin{align*} x_n &= \frac{1}{\prod_{i = 0}^{n - 1}\cos\theta_i}(x_0\cos\theta - y_0\sin\theta) \quad (\text{其中 }i \text{ 从 }0 \text{ 至 }n - 1)\\ y_n &= \frac{1}{\prod_{i = 0}^{n - 1}\cos\theta_i}(y_0\cos\theta + x_0\sin\theta) \quad (\text{其中 }i \text{ 从 }0 \text{ 至 }n - 1) \end{align*} xn​yn​​=∏i=0n−1​cosθi​1​(x0​cosθ−y0​sinθ)(其中 i 从 0 至 n−1)=∏i=0n−1​cosθi​1​(y0​cosθ+x0​sinθ)(其中 i 从 0 至 n−1)​

此时，已知条件有三个 x 0 x_0 x0​、 y 0 y_0 y0​ 和 θ \theta θ。通过 16 16 16 次迭代，我们可以得到 x n x_n xn​ 和 y n y_n yn​。而式中的 ∏ i = 0 n − 1 cos ⁡ θ i \prod_{i = 0}^{n - 1}\cos\theta_i ∏i=0n−1​cosθi​ 是个随 i i i 变化的值，我们可以预先将其值存入系统中。

然后，我们人为设置 x 0 = ∏ i = 0 n − 1 cos ⁡ θ i x_0 = \prod_{i = 0}^{n - 1}\cos\theta_i x0​=∏i=0n−1​cosθi​， y 0 = 0 y_0 = 0 y0​=0，则根据等式， x n = cos ⁡ θ x_n = \cos\theta xn​=cosθ， y n = sin ⁡ θ y_n = \sin\theta yn​=sinθ。其中 1 ∏ i = 0 n − 1 cos ⁡ θ i \frac{1}{\prod_{i = 0}^{n - 1}\cos\theta_i} ∏i=0n−1​cosθi​1​ 的值我们也同样预先存入系统中。如此，我们就实现了正弦和余弦操作了。

二、CORDIC 的具体操作流程介绍

1、CORDIC 的旋转模式

由于算法较复杂，一休哥再总结一些具体的操作流程。

1、 设置迭代次数为 16 16 16，则 x 0 = 0.607253 x_0 = 0.607253 x0​=0.607253， y 0 = 0 y_0 = 0 y0​=0，并输入待计算的角度 θ \theta θ， θ \theta θ 在 [ − 99. 7 ∘ , 99. 7 ∘ ] [-99.7^{\circ}, 99.7^{\circ}] [−99.7∘,99.7∘] 范围内。

2、 根据三个迭代公式进行迭代， i i i 从 0 0 0 至 15 15 15：

x i + 1 = x i − d i y i 2 − i y i + 1 = y i + d i x i 2 − i z i + 1 = z i − d i θ i \begin{align*} x_{i + 1} &= x_i - d_iy_i2^{-i}\\ y_{i + 1} &= y_i + d_ix_i2^{-i}\\ z_{i + 1} &= z_i - d_i\theta_i \end{align*} xi+1​yi+1​zi+1​​=xi​−di​yi​2−i=yi​+di​xi​2−i=zi​−di​θi​​

注： z 0 = θ z_0 = \theta z0​=θ， d i d_i di​ 与 z i z_i zi​ 同符号。

3、 经过 16 16 16 次迭代计算后，得到的 x 16 x_{16} x16​ 和 y 16 y_{16} y16​ 分别为 cos ⁡ θ \cos\theta cosθ 和 sin ⁡ θ \sin\theta sinθ。

至此，关于 CORDIC 的三角函数 cos ⁡ θ \cos\theta cosθ 和 sin ⁡ θ \sin\theta sinθ 的计算原理讲解结束。

关于 CORDIC 算法计算三角函数 cos ⁡ θ \cos\theta cosθ 和 sin ⁡ θ \sin\theta sinθ 的 MATLAB 代码 如下所示：

close all;
clear;
clc;
% 初始化
die = 16;%迭代次数
x = zeros(die+1,1);
y = zeros(die+1,1);
z = zeros(die+1,1);
x(1) = 0.607253;%初始设置
z(1) = pi/4;%待求角度θ
%迭代操作
for i = 1:die
    if z(i) >= 0
        d = 1;
    else
        d = -1;
    end
    x(i+1) = x(i) - d*y(i)*(2^(-(i-1)));
    y(i+1) = y(i) + d*x(i)*(2^(-(i-1)));
    z(i+1) = z(i) - d*atan(2^(-(i-1)));
end
cosa = vpa(x(17),10)
sina = vpa(y(17),10)
c = vpa(z(17),10)

2、CORDIC的向量模式

讲完了旋转模式后，我们接着讲讲向量模式下的圆坐标系。

在这里，我们需从头来过了，假设在 x y xy xy 坐标系中有一个点 P 0 ( x 0 , y 0 ) P0(x_0,y_0) P0(x0​,y0​) ，将 P 0 P0 P0 点绕原点旋转 θ \theta θ 角后得到点 P n ( x n , 0 ) Pn(x_n,0) Pn(xn​,0) ， θ \theta θ 在 [ − 99. 7 ∘ , 99. 7 ∘ ] [-99.7^{\circ},99.7^{\circ}] [−99.7∘,99.7∘] 范围内。

于是可以得到 P 0 P0 P0 和 P n Pn Pn 的关系：

x n = x 0 cos ⁡ θ − y 0 sin ⁡ θ = cos ⁡ θ ( x 0 − y 0 tan ⁡ θ ) y n = y 0 cos ⁡ θ + x 0 sin ⁡ θ = cos ⁡ θ ( y 0 + x 0 tan ⁡ θ ) = 0 \begin{align*} x_n &= x_0\cos\theta - y_0\sin\theta = \cos\theta(x_0 - y_0\tan\theta)\\ y_n &= y_0\cos\theta + x_0\sin\theta = \cos\theta(y_0 + x_0\tan\theta) = 0 \end{align*} xn​yn​​=x0​cosθ−y0​sinθ=cosθ(x0​−y0​tanθ)=y0​cosθ+x0​sinθ=cosθ(y0​+x0​tanθ)=0​

如何得到 P n ( x n , y n ) Pn(x_n,y_n) Pn(xn​,yn​) 一直是我们的目标。而此时，我们还是列出那几个等式：

根据三个迭代公式进行迭代， i i i 从 0 0 0 至 15 15 15 ：

x i + 1 = x i − d i y i 2 − i y i + 1 = y i + d i x i 2 − i z i + 1 = z i − d i θ i \begin{align*} x_{i + 1} &= x_i - d_iy_i2^{-i}\\ y_{i + 1} &= y_i + d_ix_i2^{-i}\\ z_{i + 1} &= z_i - d_i\theta_i \end{align*} xi+1​yi+1​zi+1​​=xi​−di​yi​2−i=yi​+di​xi​2−i=zi​−di​θi​​

不过此时我们尝试改变初始条件：

设置迭代次数为 16 16 16 ，则 x 0 = x x_0 = x x0​=x ， y 0 = y y_0 = y y0​=y ， z 0 = 0 z_0 = 0 z0​=0 ， d i d_i di​ 与 y i y_i yi​ 的符号相反。表示，经过 n n n 次旋转，使 P n Pn Pn 靠近 x x x 轴。

因此，当迭代结束之后， P n Pn Pn 将近似接近 x x x 轴，此时 y n = 0 y_n = 0 yn​=0 ，可知旋转了 θ \theta θ ，即 z n = θ = arctan ⁡ ( y / x ) z_n = \theta = \arctan(y/x) zn​=θ=arctan(y/x) 。

而

x n = 1 ∏ i = 0 n − 1 cos ⁡ θ i ( x 0 cos ⁡ θ − y 0 sin ⁡ θ ) ( 其中  i  从  0  至  n − 1 ) y n = 1 ∏ i = 0 n − 1 cos ⁡ θ i ( y 0 cos ⁡ θ + x 0 sin ⁡ θ ) ( 其中  i  从  0  至  n − 1 ) \begin{align*} x_n &= \frac{1}{\prod_{i = 0}^{n - 1}\cos\theta_i}(x_0\cos\theta - y_0\sin\theta) \quad (\text{其中 }i \text{ 从 }0 \text{ 至 }n - 1)\\ y_n &= \frac{1}{\prod_{i = 0}^{n - 1}\cos\theta_i}(y_0\cos\theta + x_0\sin\theta) \quad (\text{其中 }i \text{ 从 }0 \text{ 至 }n - 1) \end{align*} xn​yn​​=∏i=0n−1​cosθi​1​(x0​cosθ−y0​sinθ)(其中 i 从 0 至 n−1)=∏i=0n−1​cosθi​1​(y0​cosθ+x0​sinθ)(其中 i 从 0 至 n−1)​

因此，可得 y cos ⁡ θ + x sin ⁡ θ = 0 y\cos\theta + x\sin\theta = 0 ycosθ+xsinθ=0 ，

x n = 1 ∏ i = 0 n − 1 cos ⁡ θ i ( x cos ⁡ θ − y sin ⁡ θ ) = 1 ∏ i = 0 n − 1 cos ⁡ θ i { [ ( x cos ⁡ θ − y sin ⁡ θ ) 2 ] 1 2 } = 1 ∏ i = 0 n − 1 cos ⁡ θ i { [ x 2 cos ⁡ 2 θ + y 2 sin ⁡ 2 θ − 2 x y sin ⁡ θ cos ⁡ θ ] 1 2 } = 1 ∏ i = 0 n − 1 cos ⁡ θ i { [ x 2 cos ⁡ 2 θ + y 2 sin ⁡ 2 θ + y 2 cos ⁡ 2 θ + x 2 sin ⁡ 2 θ ] 1 2 } = 1 ∏ i = 0 n − 1 cos ⁡ θ i { [ x 2 + y 2 ] 1 2 } \begin{align*} x_n &= \frac{1}{\prod_{i = 0}^{n - 1}\cos\theta_i}(x\cos\theta - y\sin\theta)\\ &= \frac{1}{\prod_{i = 0}^{n - 1}\cos\theta_i}\{ [ (x\cos\theta - y\sin\theta)^2]^{\frac{1}{2}}\}\\ &= \frac{1}{\prod_{i = 0}^{n - 1}\cos\theta_i}\{ [ x^2\cos^2\theta + y^2\sin^2\theta - 2xy\sin\theta\cos\theta]^{\frac{1}{2}}\}\\ &= \frac{1}{\prod_{i = 0}^{n - 1}\cos\theta_i}\{ [ x^2\cos^2\theta + y^2\sin^2\theta + y^2\cos^2\theta + x^2\sin^2\theta]^{\frac{1}{2}}\}\\ &= \frac{1}{\prod_{i = 0}^{n - 1}\cos\theta_i}\{ [ x^2 + y^2]^{\frac{1}{2}}\} \end{align*} xn​​=∏i=0n−1​cosθi​1​(xcosθ−ysinθ)=∏i=0n−1​cosθi​1​{[(xcosθ−ysinθ)2]21​}=∏i=0n−1​cosθi​1​{[x2cos2θ+y2sin2θ−2xysinθcosθ]21​}=∏i=0n−1​cosθi​1​{[x2cos2θ+y2sin2θ+y2cos2θ+x2sin2θ]21​}=∏i=0n−1​cosθi​1​{[x2+y2]21​}​

由上可以知道，我们通过迭代，就算出了反正切函数 z n = θ = arctan ⁡ ( y / x ) z_n = \theta = \arctan(y/x) zn​=θ=arctan(y/x) ，以及向量 O P 0 → ( x , y ) \overrightarrow{OP_0}(x,y) OP0​ ​(x,y) 的长度 d = x n ⋅ ∏ i = 0 n − 1 cos ⁡ θ i d = x_n \cdot \prod_{i = 0}^{n - 1}\cos\theta_i d=xn​⋅∏i=0n−1​cosθi​ 。

关于反正切函数，一休哥要多啰嗦几句了，由于 θ \theta θ 在 [ − 99. 7 ∘ , 99. 7 ∘ ] [-99.7^{\circ},99.7^{\circ}] [−99.7∘,99.7∘] 范围内，所以我们输入向量 O P 0 → ( x , y ) \overrightarrow{OP_0}(x,y) OP0​ ​(x,y) 时，需要保证其在第一、四象限。

关于 CORDIC 算法计算反三角函数 arctan ⁡ θ \arctan\theta arctanθ 的 MATLAB 代码如下所示：

close all;
clear;
clc;
% 初始化
die = 16;%迭代次数
x = zeros(die+1,1);
y = zeros(die+1,1);
z = zeros(die+1,1);
x(1) = 100;%初始设置
y(1) = 200;%初始设置
k = 0.607253;%初始设置
%迭代操作
for i = 1:die
    if y(i) >= 0
        d = -1;
    else
        d = 1;
    end
    x(i+1) = x(i) - d*y(i)*(2^(-(i-1)));
    y(i+1) = y(i) + d*x(i)*(2^(-(i-1)));
    z(i+1) = z(i) - d*atan(2^(-(i-1)));
end
d = vpa(x(17)*k,10)
a = vpa(y(17),10)
c = vpa(rad2deg(z(17)),10)

三、CORDIC 的旋转模式——Verilog 仿真

一休哥在编写CORDIC算法时，采用了 16 级流水线，仿真效果十分明显。以下是顶层文件的代码。

为了避免浮点运算，为了满足精度要求，一休哥对每个变量都放大了 2 16 2^{16} 216 倍，并且引入了有符号型 reg \text{reg} reg 和算术右移。

关于 Verilog 代码的编写，一休哥已经不想多说了，因为代码是完全符合我之前所讲的 CORDIC 的原理与 MATLAB 仿真代码。相信大家在看完本文的前两个部分之后，对 Verilog 的理解应该不是难事儿。

module Cordic_Test
(
    CLK_50M,RST_N,
    Phase,
    Sin,Cos,Error
);

input                       CLK_50M;
input                       RST_N;
input       [31:0]          Phase;
output      [31:0]          Sin;
output      [31:0]          Cos;
output      [31:0]          Error;

`define rot0  32'd2949120       //45度*2^16
`define rot1  32'd1740992       //26.5651度*2^16
`define rot2  32'd919872        //14.0362度*2^16
`define rot3  32'd466944        //7.1250度*2^16
`define rot4  32'd234368        //3.5763度*2^16
`define rot5  32'd117312        //1.7899度*2^16
`define rot6  32'd58688         //0.8952度*2^16
`define rot7  32'd29312         //0.4476度*2^16
`define rot8  32'd14656         //0.2238度*2^16
`define rot9  32'd7360          //0.1119度*2^16
`define rot10 32'd3648          //0.0560度*2^16
`define rot11 32'd1856          //0.0280度*2^16
`define rot12 32'd896           //0.0140度*2^16
`define rot13 32'd448           //0.0070度*2^16
`define rot14 32'd256           //0.0035度*2^16
`define rot15 32'd128           //0.0018度*2^16

parameter Pipeline = 16;
parameter K = 32'h09b74;    //K=0.607253*2^16,32'h09b74,

reg signed  [31:0]      Sin;
reg signed  [31:0]      Cos;
reg signed  [31:0]      Error;
reg signed  [31:0]      x0=0,y0=0,z0=0;
reg signed  [31:0]      x1=0,y1=0,z1=0;
reg signed  [31:0]      x2=0,y2=0,z2=0;
reg signed  [31:0]      x3=0,y3=0,z3=0;
reg signed  [31:0]      x4=0,y4=0,z4=0;
reg signed  [31:0]      x5=0,y5=0,z5=0;
reg signed  [31:0]      x6=0,y6=0,z6=0;
reg signed  [31:0]      x7=0,y7=0,z7=0;
reg signed  [31:0]      x8=0,y8=0,z8=0;
reg signed  [31:0]      x9=0,y9=0,z9=0;
reg signed  [31:0]      x10=0,y10=0,z10=0;
reg signed  [31:0]      x11=0,y11=0,z11=0;
reg signed  [31:0]      x12=0,y12=0,z12=0;
reg signed  [31:0]      x13=0,y13=0,z13=0;
reg signed  [31:0]      x14=0,y14=0,z14=0;
reg signed  [31:0]      x15=0,y15=0,z15=0;
reg signed  [31:0]      x16=0,y16=0,z16=0;
reg         [ 1:0]      Quadrant [Pipeline:0];

always @ (posedge CLK_50M or negedge RST_N)
begin
    if(!RST_N)
    begin
        x0 <= 1'b0;                         
        y0 <= 1'b0;
        z0 <= 1'b0;
    end
    else
    begin
        x0 <= K;
        y0 <= 32'd0;
        z0 <= Phase[15:0] << 16;
    end
end

always @ (posedge CLK_50M or negedge RST_N)
begin
    if(!RST_N)
    begin
        x1 <= 1'b0;                         
        y1 <= 1'b0;
        z1 <= 1'b0;
    end
    else if(z0[31])
    begin
      x1 <= x0 + y0;
      y1 <= y0 - x0;
      z1 <= z0 + `rot0;
    end
    else
    begin
      x1 <= x0 - y0;
      y1 <= y0 + x0;
      z1 <= z0 - `rot0;
    end
end

always @ (posedge CLK_50M or negedge RST_N)
begin
    if(!RST_N)
    begin
        x2 <= 1'b0;                         
        y2 <= 1'b0;
        z2 <= 1'b0;
    end
    else if(z1[31])
   begin
        x2 <= x1 + (y1 >>> 1);
        y2 <= y1 - (x1 >>> 1);
        z2 <= z1 + `rot1;
   end
   else
   begin
       x2 <= x1 - (y1 >>> 1);
       y2 <= y1 + (x1 >>> 1);
       z2 <= z1 - `rot1;
   end
end

always @ (posedge CLK_50M or negedge RST_N)
begin
    if(!RST_N)
    begin
        x3 <= 1'b0;                         
        y3 <= 1'b0;
        z3 <= 1'b0;
    end
    else if(z2[31])
   begin
       x3 <= x2 + (y2 >>> 2);
       y3 <= y2 - (x2 >>> 2);
       z3 <= z2 + `rot2;
   end
   else
   begin
       x3 <= x2 - (y2 >>> 2);
       y3 <= y2 + (x2 >>> 2);
       z3 <= z2 - `rot2;
   end
end

always @ (posedge CLK_50M or negedge RST_N)
begin
    if(!RST_N)
    begin
        x4 <= 1'b0;                         
        y4 <= 1'b0;
        z4 <= 1'b0;
    end
    else if(z3[31])
   begin
       x4 <= x3 + (y3 >>> 3);
       y4 <= y3 - (x3 >>> 3);
       z4 <= z3 + `rot3;
   end
   else
   begin
       x4 <= x3 - (y3 >>> 3);
       y4 <= y3 + (x3 >>> 3);
       z4 <= z3 - `rot3;
   end
end

always @ (posedge CLK_50M or negedge RST_N)
begin
    if(!RST_N)
    begin
        x5 <= 1'b0;                         
        y5 <= 1'b0;
        z5 <= 1'b0;
    end
    else if(z4[31])
   begin
       x5 <= x4 + (y4 >>> 4);
       y5 <= y4 - (x4 >>> 4);
       z5 <= z4 + `rot4;
   end
   else
   begin
       x5 <= x4 - (y4 >>> 4);
       y5 <= y4 + (x4 >>> 4);
       z5 <= z4 - `rot4;
   end
end

always @ (posedge CLK_50M or negedge RST_N)
begin
    if(!RST_N)
    begin
        x6 <= 1'b0;                         
        y6 <= 1'b0;
        z6 <= 1'b0;
    end
    else if(z5[31])
   begin
       x6 <= x5 + (y5 >>> 5);
       y6 <= y5 - (x5 >>> 5);
       z6 <= z5 + `rot5;
   end
   else
   begin
       x6 <= x5 - (y5 >>> 5);
       y6 <= y5 + (x5 >>> 5);
       z6 <= z5 - `rot5;
   end
end

always @ (posedge CLK_50M or negedge RST_N)
begin
    if(!RST_N)
    begin
        x7 <= 1'b0;                         
        y7 <= 1'b0;
        z7 <= 1'b0;
    end
    else if(z6[31])
   begin
       x7 <= x6 + (y6 >>> 6);
       y7 <= y6 - (x6 >>> 6);
       z7 <= z6 + `rot6;
   end
   else
   begin
       x7 <= x6 - (y6 >>> 6);
       y7 <= y6 + (x6 >>> 6);
       z7 <= z6 - `rot6;
   end
end

always @ (posedge CLK_50M or negedge RST_N)
begin
    if(!RST_N)
    begin
        x8 <= 1'b0;                         
        y8 <= 1'b0;
        z8 <= 1'b0;
    end
    else if(z7[31])
   begin
       x8 <= x7 + (y7 >>> 7);
       y8 <= y7 - (x7 >>> 7);
       z8 <= z7 + `rot7;
   end
   else
   begin
       x8 <= x7 - (y7 >>> 7);
       y8 <= y7 + (x7 >>> 7);
       z8 <= z7 - `rot7;
   end
end

always @ (posedge CLK_50M or negedge RST_N)
begin
    if(!RST_N)
    begin
        x9 <= 1'b0;                         
        y9 <= 1'b0;
        z9 <= 1'b0;
    end
    else if(z8[31])
   begin
       x9 <= x8 + (y8 >>> 8);
       y9 <= y8 - (x8 >>> 8);
       z9 <= z8 + `rot8;
   end
   else
   begin
       x9 <= x8 - (y8 >>> 8);
       y9 <= y8 + (x8 >>> 8);
       z9 <= z8 - `rot8;
   end
end

always @ (posedge CLK_50M or negedge RST_N)
begin
    if(!RST_N)
    begin
        x10 <= 1'b0;                        
        y10 <= 1'b0;
        z10 <= 1'b0;
    end
    else if(z9[31])
   begin
       x10 <= x9 + (y9 >>> 9);
       y10 <= y9 - (x9 >>> 9);
       z10 <= z9 + `rot9;
   end
   else
   begin
       x10 <= x9 - (y9 >>> 9);
       y10 <= y9 + (x9 >>> 9);
       z10 <= z9 - `rot9;
   end
end

always @ (posedge CLK_50M or negedge RST_N)
begin
    if(!RST_N)
    begin
        x11 <= 1'b0;                        
        y11 <= 1'b0;
        z11 <= 1'b0;
    end
    else if(z10[31])
   begin
       x11 <= x10 + (y10 >>> 10);
       y11 <= y10 - (x10 >>> 10);
       z11 <= z10 + `rot10;
   end
   else
   begin
       x11 <= x10 - (y10 >>> 10);
       y11 <= y10 + (x10 >>> 10);
       z11 <= z10 - `rot10;
   end
end

always @ (posedge CLK_50M or negedge RST_N)
begin
    if(!RST_N)
    begin
        x12 <= 1'b0;                        
        y12 <= 1'b0;
        z12 <= 1'b0;
    end
    else if(z11[31])
   begin
       x12 <= x11 + (y11 >>> 11);
       y12 <= y11 - (x11 >>> 11);
       z12 <= z11 + `rot11;
   end
   else
   begin
       x12 <= x11 - (y11 >>> 11);
       y12 <= y11 + (x11 >>> 11);
       z12 <= z11 - `rot11;
   end
end

always @ (posedge CLK_50M or negedge RST_N)
begin
    if(!RST_N)
    begin
        x13 <= 1'b0;                        
        y13 <= 1'b0;
        z13 <= 1'b0;
    end
    else if(z12[31])
   begin
       x13 <= x12 + (y12 >>> 12);
       y13 <= y12 - (x12 >>> 12);
       z13 <= z12 + `rot12;
   end
   else
   begin
       x13 <= x12 - (y12 >>> 12);
       y13 <= y12 + (x12 >>> 12);
       z13 <= z12 - `rot12;
   end
end

always @ (posedge CLK_50M or negedge RST_N)
begin
    if(!RST_N)
    begin
        x14 <= 1'b0;                        
        y14 <= 1'b0;
        z14 <= 1'b0;
    end
    else if(z13[31])
   begin
       x14 <= x13 + (y13 >>> 13);
       y14 <= y13 - (x13 >>> 13);
       z14 <= z13 + `rot13;
   end
   else
   begin
       x14 <= x13 - (y13 >>> 13);
       y14 <= y13 + (x13 >>> 13);
       z14 <= z13 - `rot13;
   end
end

always @ (posedge CLK_50M or negedge RST_N)
begin
    if(!RST_N)
    begin
        x15 <= 1'b0;                        
        y15 <= 1'b0;
        z15 <= 1'b0;
    end
    else if(z14[31])
   begin
       x15 <= x14 + (y14 >>> 14);
       y15 <= y14 - (x14 >>> 14);
       z15 <= z14 + `rot14;
   end
   else
   begin
       x15 <= x14 - (y14 >>> 14);
       y15 <= y14 + (x14 >>> 14);
       z15 <= z14 - `rot14;
   end
end

always @ (posedge CLK_50M or negedge RST_N)
begin
    if(!RST_N)
    begin
        x16 <= 1'b0;                        
        y16 <= 1'b0;
        z16 <= 1'b0;
    end
    else if(z15[31])
   begin
       x16 <= x15 + (y15 >>> 15);
       y16 <= y15 - (x15 >>> 15);
       z16 <= z15 + `rot15;
   end
   else
   begin
       x16 <= x15 - (y15 >>> 15);
       y16 <= y15 + (x15 >>> 15);
       z16 <= z15 - `rot15;
   end
end

always @ (posedge CLK_50M or negedge RST_N)
begin
    if(!RST_N)
    begin
        Quadrant[0] <= 1'b0;
        Quadrant[1] <= 1'b0;
        Quadrant[2] <= 1'b0;
        Quadrant[3] <= 1'b0;
        Quadrant[4] <= 1'b0;
        Quadrant[5] <= 1'b0;
        Quadrant[6] <= 1'b0;
        Quadrant[7] <= 1'b0;
        Quadrant[8] <= 1'b0;
        Quadrant[9] <= 1'b0;
        Quadrant[10] <= 1'b0;
        Quadrant[11] <= 1'b0;
        Quadrant[12] <= 1'b0;
        Quadrant[13] <= 1'b0;
        Quadrant[14] <= 1'b0;
        Quadrant[15] <= 1'b0;
        Quadrant[16] <= 1'b0;
    end
    else
    begin
        Quadrant[0] <= Phase[17:16];
        Quadrant[1] <= Quadrant[0];
        Quadrant[2] <= Quadrant[1];
        Quadrant[3] <= Quadrant[2];
        Quadrant[4] <= Quadrant[3];
        Quadrant[5] <= Quadrant[4];
        Quadrant[6] <= Quadrant[5];
        Quadrant[7] <= Quadrant[6];
        Quadrant[8] <= Quadrant[7];
        Quadrant[9] <= Quadrant[8];
        Quadrant[10] <= Quadrant[9];
        Quadrant[11] <= Quadrant[10];
        Quadrant[12] <= Quadrant[11];
        Quadrant[13] <= Quadrant[12];
        Quadrant[14] <= Quadrant[13];
        Quadrant[15] <= Quadrant[14];
        Quadrant[16] <= Quadrant[15];
    end
end

always @ (posedge CLK_50M or negedge RST_N)
begin
    if(!RST_N)
    begin
        Cos <= 1'b0;
        Sin <= 1'b0;
        Error <= 1'b0;
    end
    else
    begin
        Error <= z16;
        case(Quadrant[16])
            2'b00: //if the Phase is in first Quadrant,the Sin(X)=Sin(A),Cos(X)=Cos(A)
                begin
                    Cos <= x16;
                    Sin <= y16;
                end
            2'b01: //if the Phase is in second Quadrant,the Sin(X)=Sin(A+90)=CosA,Cos(X)=Cos(A+90)=-SinA
                begin
                    Cos <= ~(y16) + 1'b1;//-Sin
                    Sin <= x16;//Cos
                end
            2'b10: //if the Phase is in third Quadrant,the Sin(X)=Sin(A+180)=-SinA,Cos(X)=Cos(A+180)=-CosA
                begin
                    Cos <= ~(x16) + 1'b1;//-Cos
                    Sin <= ~(y16) + 1'b1;//-Sin
                end
            2'b11: //if the Phase is in forth Quadrant,the Sin(X)=Sin(A+270)=-CosA,Cos(X)=Cos(A+270)=SinA
                begin
                    Cos <= y16;//Sin
                    Sin <= ~(x16) + 1'b1;//-Cos
                end
        endcase
    end
end

endmodule

以下是testbench文件代码

`timescale 1 ps/ 1 ps
 
module Cordic_Test_tb;
 
// Inputs
reg                         CLK_50M;
reg                         RST_N;
reg             [15:0]      cnt;
reg             [15:0]      cnt_n;
reg             [31:0]      Phase;
reg             [31:0]      Phase_n;
wire            [31:0]      Sin;
wire            [31:0]      Cos;
wire            [31:0]      Error;
 
// Instantiate the Unit Under Test (UUT)
Cordic_Test                 uut 
(
    .CLK_50M                (CLK_50M    ),
    .RST_N                  (RST_N      ),
    .Phase                  (Phase      ),
    .Sin                    (Sin        ),
    .Cos                    (Cos        ),
    .Error                  (Error      )
);
 
initial
begin
    #0 CLK_50M = 1'b0;
    #10000 RST_N = 1'b0;
    #10000 RST_N = 1'b1;
    #10000000 $stop;
end 
always #10000 
begin
    CLK_50M = ~CLK_50M;
end
always @ (posedge CLK_50M or negedge RST_N)
begin
    if(!RST_N)
        cnt <= 1'b0;
    else
        cnt <= cnt_n;
end
 
always @ (*)
begin
    if(cnt == 16'd359)
        cnt_n = 1'b0;
    else
        cnt_n = cnt + 1'b1;
end
//生成相位0-359度,Phase[17:16]为相位的象限,Phase[15:10]为相位的值
always @ (posedge CLK_50M or negedge RST_N)
begin
    if(!RST_N)
        Phase <= 1'b0;
    else
        Phase <= Phase_n;
end
 
always @ (*)
begin
    if(cnt <= 16'd90)
        Phase_n = cnt;
    else if(cnt > 16'd90 && cnt <= 16'd180)
        Phase_n = {2'd01,cnt - 16'd90};
    else if(cnt > 16'd180 && cnt <= 16'd270)
        Phase_n = {2'd10,cnt - 16'd180};
    else if(cnt > 16'd270)
        Phase_n = {2'd11,cnt - 16'd270};
end
 
endmodule

最后来一张效果图，可以发现，我们的 16 级流水线已经正常的运行起来了，由于我们仿真输入的相位值为 0-359 度循环，因此 sin 和 cos 也循环了～～～

这里写图片描述

CORDIC 算法详解及 FPGA 实现

ngany 于 2021-05-30 18:52:04 发布

CORDIC 算法详解

1 平面坐标系旋转

CORDIC 算法的思想是通过迭代的方法，使得累计旋转过的角度的和无限接近目标角度。它是一种数值计算逼近的方法，运算只有移位和加减。

通过圆坐标系可以了解 CORDIC 算法的基本思想，如图 1 所示，初始向量 ( x 1 , y 1 ) \left (x_{1},y_{1} \right) (x1​,y1​)旋转 θ \theta θ角度之后得到向量 ( x 2 , y 2 ) \left ( x_{2},y_{2} \right) (x2​,y2​)，两者之间满足（公式 1）关系。

img

图 1 CORDIC 算法原理示意图
$$\left{ x2=x1cosθn−y1sinθny2=y1cosθn−x1sinθn\begin {matrix} x_{2} = x_{1} \cos\theta_{n} - y_{1} \sin\theta_{n} \ y_{2} = y_{1} \cos\theta_{n} - x_{1}\sin\theta_{n} \ \end {matrix} \right.$$

通过提取因数 cos ⁡ θ \cos\theta cosθ，方程式可以改写成下面的形式：

x 2 = x 1 cos ⁡ θ − y 1 sin ⁡ θ = cos ⁡ θ ( x 1 − y 1 tan ⁡ θ ) y 2 = y 1 cos ⁡ θ + x 1 sin ⁡ θ = cos ⁡ θ ( y 1 + x 1 tan ⁡ θ ) \begin{align*} x_{2} &= x_{1} \cos\theta - y_{1} \sin\theta = \cos\theta(x_{1} - y_{1}\tan\theta) \\ y_{2} &= y_{1} \cos\theta + x_{1} \sin\theta = \cos\theta(y_{1} + x_{1}\tan\theta) \end{align*} x2​y2​​=x1​cosθ−y1​sinθ=cosθ(x1​−y1​tanθ)=y1​cosθ+x1​sinθ=cosθ(y1​+x1​tanθ)​

2 伪旋转

如果去掉 cos ⁡ θ \cos\theta cosθ，我们可以的到伪旋转方程式（公式 3），即旋转的角度是正确的，但是 x x x， y y y 的值增加了 cos ⁡ − 1 θ \cos^{- 1}\theta cos−1θ 倍，由于 cos ⁡ − 1 θ > 1 \cos^{- 1}\theta > 1 cos−1θ>1，所以模值变大。这里并不能通过数学方法去除 cos ⁡ θ \cos\theta cosθ项，但是可以化简坐标平面旋转的计算。

x ^ 2 = x 1 − y 1 tan ⁡ θ y ^ 2 = y 1 + x 1 tan ⁡ θ \begin{align*} \widehat{x}_{2} &= x_{1} - y_{1}\tan\theta \\ \widehat{y}_{2} &= y_{1} + x_{1}\tan\theta \end{align*} x 2​y ​2​​=x1​−y1​tanθ=y1​+x1​tanθ​

img

图 2 去除 cos ⁡ θ \cos\theta cosθ后伪坐标系

3 CORDIC 方法

这里为了便于硬件的计算，采用迭代的思想，旋转角度 θ \theta θ 可以通过若干步实现，每一步旋转一个角度 θ i \theta^{i} θi 。并且约定每一次旋转的角度的正切值为 2 的倍数，即 tan ⁡ θ i = 2 − i {\tan\theta^{i} = 2}^{- i} tanθi=2−i 。下面表格是用于 CORDIC 算法中每个迭代 i i i 的旋转角度。这样，乘以正切值就可以变成移位操作。

表 1 迭代角度 θ i \theta^{i} θi正切值

i i i θ i \theta^{i} θi（degree） tan ⁡ θ i = 2 − i {\tan\theta^{i} = 2}^{- i} tanθi=2−i
045.01
126.5550511770.5
214.0362434670.25
37.1250163480.125
43.5763343740.0625

4 角度累加器

引入控制算子 d i d_{i} di​，用于确定旋转的方向。其中第 i i i 步顺时针旋转时 d i = − 1 d_{i} = - 1 di​=−1，逆时针旋转 d i = 1 d_{i} = 1 di​=1。前面的伪旋转坐标方程现在可以表示为：

x i + 1 = x i − d i 2 − i y i y i + 1 = y i + d i 2 − i x i \begin{align*} x_{i + 1} &= x_{i} - d_{i}2^{- i}y_{i} \\ y_{i + 1} &= y_{i} + d_{i}2^{- i}x_{i} \end{align*} xi+1​yi+1​​=xi​−di​2−iyi​=yi​+di​2−ixi​​

这里，我们引入第三个方程，被称为角度累加器，用来在每次迭代过程中追踪累加的旋转角度，表示第 n n n 次旋转后剩下未旋转的角度，定义为

z i + 1 = z i − d i θ i z_{i + 1} = z_{i} - d_{i} \theta^{i} zi+1​=zi​−di​θi，其中 d i = ± 1 d_{i} = \pm 1 di​=±1

5 移位 - 加法算法

因此，原始的算法现在已经被简化为使用向量的伪旋转方程式（公式 6）来表示迭代（移位 - 加法）算法。

  • 2 次移位；
    （ 2 − i 2^{- i} 2−i用移位实现，每右移 n 位就把原来数值乘以 2 的 - i 次方了）

  • 1 次查找表；（每一次迭代都会有一个固定角度 θ i \theta^{i} θi的累加，这个角度是 2 − i 2^{- i} 2−i对应的角度值，使用查表实现。 θ i = arc tan ⁡ 2 − i \theta^{i} = {\text {arc}\tan} 2^{- i} θi=arctan2−i）

  • 3 次加法；（x，y，z 的累加，共三次）

x i + 1 = x i − d i 2 − i y i y i + 1 = y i + d i 2 − i x i z i + 1 = z i − d i θ i \begin{align*} x_{i + 1} &= x_{i} - d_{i}2^{- i}y_{i} \\ y_{i + 1} &= y_{i} + d_{i}2^{- i}x_{i} \\ z_{i + 1} &= z_{i} - d_{i}\theta_{i} \end{align*} xi+1​yi+1​zi+1​​=xi​−di​2−iyi​=yi​+di​2−ixi​=zi​−di​θi​​

6 伸缩因子

当简化算法以伪旋转时， cos ⁡ θ \cos\theta cosθ项被忽略，这样 ( x n , y n ) \left ( x_{n},y_{n} \right) (xn​,yn​) 就被伸缩了 K n K_{n} Kn​ 倍。如果迭代次数已知，可以预先计算伸缩因子 K n K_{n} Kn​。同样 1 / K n 1/K_{n} 1/Kn​也可被预先计算以得到 ( x n , y n ) \left ( x_{n},y_{n} \right) (xn​,yn​) 的真值。

K n = ∏ n 1 cos ⁡ θ i = ∏ n ( 1 + 2 − 2 i ) K_{n} = \prod_{n}^{}\frac {1}{\cos\theta^{i}} = \prod_{n}^{}\left ( \sqrt {1 + 2^{- 2i}} \right) Kn​=∏n​cosθi1​=∏n​(1+2−2i ​)

这里 K n − > 1.64676 K_{n} - > 1.64676 Kn​−>1.64676，当 n − > ∞ n - > \infty n−>∞

1 / K n − > 0.607253 1/K_{n} - > 0.607253 1/Kn​−>0.607253，当 n − > ∞ n - > \infty n−>∞

7 旋转模式

CORDIC 方法有两种模式，旋转模式和向量模式。工作模式决定了控制算子 d i d_{i} di​ 的条件。在旋转模式中，旋转剩余角度初始变量 z 0 = θ z_{0} = \theta z0​=θ，经过 n 次旋转后，使 z 0 = 0 z_{0} = 0 z0​=0。整个旋转过程可以表示为一系列旋转角度不断偏摆，从而逼近所需旋转角度的迭代过程。n 次迭代后可以得到：

x n = K n ( x 0 cos ⁡ z 0 − y 0 sin ⁡ z 0 ) y n = K n ( y 0 cos ⁡ z 0 + x 0 sin ⁡ z 0 ) z n = 0 \begin{align*} x_{n} &= K_{n}(x_{0} \cos z_{0} - y_{0} \sin z_{0}) \\ y_{n} &= K_{n}(y_{0} \cos z_{0} + x_{0} \sin z_{0}) \\ z_{n} &= 0 \end{align*} xn​yn​zn​​=Kn​(x0​cosz0​−y0​sinz0​)=Kn​(y0​cosz0​+x0​sinz0​)=0​

通过设置 x 0 = 1 / K n = 0.6073 x_{0} = 1/K_{n} = 0.6073 x0​=1/Kn​=0.6073 和 y 0 = 0 y_{0} = 0 y0​=0 可以计算 cos ⁡ z 0 \cos z_{0} cosz0​和 sin ⁡ z 0 \sin z_{0} sinz0​。分析可知 CORDIC 算法在圆周系统旋转模式下可以用来计算一个输入角的正弦值、余弦值。

x n = K n ( x 0 cos ⁡ z 0 − y 0 sin ⁡ z 0 ) = cos ⁡ θ y n = K n ( y 0 cos ⁡ z 0 + x 0 sin ⁡ z 0 ) = sin ⁡ θ \begin{align*} x_{n} &= K_{n}(x_{0} \cos z_{0} - y_{0} \sin z_{0}) = \cos\theta \\ y_{n} &= K_{n}(y_{0} \cos z_{0} + x_{0} \sin z_{0}) = \sin\theta \end{align*} xn​yn​​=Kn​(x0​cosz0​−y0​sinz0​)=cosθ=Kn​(y0​cosz0​+x0​sinz0​)=sinθ​

img

图 3 旋转模式角度逼近

8 象限判断

对于任意角度 θ \theta θ，都可以通过旋转一系列小角度来 i i i 完成。第一次迭代旋转 4 5 ∘ 45^{\circ} 45∘，第二次迭代旋转 26. 6 ∘ 26.6^{\circ} 26.6∘，第三次迭代旋转 1 4 ∘ 14^{\circ} 14∘。很显然，每次旋转的方向都影响到最终的累计角度。在 − 99. 7 ∘ < θ < 99. 7 ∘ - 99.7^{\circ}< \theta < 99.7^{\circ} −99.7∘<θ<99.7∘ 的范围内任意角度都可以旋转。对于该角度范围之外的角度，可以通过三角函数等式转化成该范围内的角度。因此设计时通常把角度 θ \theta θ 转化到第一象限角，根据 θ \theta θ 所在的圆坐标系象限获得 x n x_{n} xn​ 和 y n y_{n} yn​。

9 溢出处理

通过表 2 可以看出，把初始 θ \theta θ 角度设置为 90°，由于 x 0 x_{0} x0​ 赋的初值存在误差（一般是偏大），最终获得的 x n x_{n} xn​ 和 y n y_{n} yn​ 是有可能大于 1 的。所以设计时需要做溢出保护。或者把 x 0 x_{0} x0​ 初值减小，这可能会牺牲精度。

表 2 θ \theta θ 接近 90° 数据溢出

idthetazyxSINCOS
014589.997000.6073019900
1126.644.99700.60730.60731990019900
211418.39700.91100.3037298509950
317.14.39700.98690.0759323382488
4-13.6-2.70300.9964-0.047432648-1555
511.80.89700.99930.014832746486
6-10.9-0.90300.9998-0.016432761-537
7-10.4-0.00301.0000-0.000832769-26
810.20.39701.00000.007032769230
910.10.19701.00010.003132770102
1010.0560.09701.00010.00123277038
1110.0270.04101.00010.0002327716

10 Verilog 代码实现

CORDIC 算法代码

module MyCORDIC (
	input clk,
	input rst_n,
	input  [15:0] theta,
	output reg [15:0] sin_theta,
	output reg [15:0] cos_theta
);

parameter Kn  = 'd19898;    // 0.607253*2^15
parameter iKn = 'd53961;    // 1.64676*2^15 

parameter arctan_0 = 8192    ;  //arctan (1/2)
parameter arctan_1 = 4836    ;  //arctan (1/2^1)
parameter arctan_2 = 2555    ;  //arctan (1/2^2)
parameter arctan_3 = 1297    ;  //arctan (1/2^3)
parameter arctan_4 = 651     ;  //arctan (1/2^4)
parameter arctan_5 = 326     ;  //arctan (1/2^5)
parameter arctan_6 = 163     ;  //arctan (1/2^6)
parameter arctan_7 = 81      ;  //arctan (1/2^7)
parameter arctan_8 = 41      ;  //arctan (1/2^8)
parameter arctan_9 = 20      ;  //arctan (1/2^9)
parameter arctan_10 = 10     ;  //arctan (1/2^10)
parameter arctan_11 = 5      ;  //arctan (1/2^11)

reg signed [15:0] x [11:0];
reg signed [15:0] y [11:0];
reg signed [15:0] z [11:0];

wire [15:0] x_tmp;
wire [15:0] y_tmp;

reg signed [15:0] theta_1;
wire [2:0] Quadrant;//theta 角所在的象限

// 象限判断
assign Quadrant = theta [15:14] + 1;

always@(*)
begin
	begin
		theta_1 <= {2'b00,theta [13:0]};
		if (Quadrant==1)
		begin
			theta_1 <= theta;
		end
		else if (Quadrant==2)
		begin
			theta_1 <= 32768 - theta;
		end
		else if (Quadrant==3)
		begin
			theta_1 <= theta - 32768;
		end
		else if (Quadrant==4)
		begin
			theta_1 <= 65536 - theta;
		end
	end
end

always@(posedge clk or negedge rst_n)
begin
	if (!rst_n)
	begin
		x [0] <= 16'd0;
		y [0] <= 16'd0;
		z [0] <= 16'd0;
	end
	else
	begin
		x [0] <= Kn;
		y [0] <= 16'd0;
		z [0] <= theta_1;
	end
end

always@(posedge clk or negedge rst_n) //i=0
begin
	if (!rst_n)
	begin
		x [1] <= 16'd0;
		y [1] <= 16'd0;
		z [1] <= 16'd0;
	end
	else
	begin
		if (z [0][15]) // 剩余角度为负数，顺时针旋转，d=-1
		begin
			x [1] <= x [0] + y [0];
			y [1] <= y [0] - x [0];
			z [1] <= z [0] + arctan_0;
		end          // 剩余角度为正数，顺时针旋转，d=+1
		else
		begin
			x [1] <= x [0] - y [0];
			y [1] <= y [0] + x [0];
			z [1] <= z [0] - arctan_0;
		end
	end
end

// >>> 符号表示算术右移，不改变符号位
always@(posedge clk or negedge rst_n) //i=1
begin
	if (!rst_n)
	begin
		x [2] <= 16'd0;
		y [2] <= 16'd0;
		z [2] <= 16'd0;
	end
	else
	begin
		if (z [1][15]) // 剩余角度为负数，顺时针旋转，d=-1
		begin
			x [2] <= x [1] + (y [1] >>> 1);
			y [2] <= y [1] - (x [1] >>> 1);
			z [2] <= z [1] + arctan_1;
		end          // 剩余角度为正数，逆时针旋转，d=+1
		else
		begin
			x [2] <= x [1] - (y [1] >>> 1);
			y [2] <= y [1] + (x [1] >>> 1);
			z [2] <= z [1] - arctan_1;
		end
	end
end

always@(posedge clk or negedge rst_n) //i=2
begin
	if (!rst_n)
	begin
		x [3] <= 16'd0;
		y [3] <= 16'd0;
		z [3] <= 16'd0;
	end
	else
	begin
		if (z [2][15]) // 剩余角度为负数，顺时针旋转，d=-1
		begin
			x [3] <= x [2] + (y [2] >>> 2);
			y [3] <= y [2] - (x [2] >>> 2);
			z [3] <= z [2] + arctan_2;
		end          // 剩余角度为正数，逆时针旋转，d=+1
		else
		begin
			x [3] <= x [2] - (y [2] >>> 2);
			y [3] <= y [2] + (x [2] >>> 2);
			z [3] <= z [2] - arctan_2;
		end
	end
end

always@(posedge clk or negedge rst_n) //i=3
begin
	if (!rst_n)
	begin
		x [4] <= 16'd0;
		y [4] <= 16'd0;
		z [4] <= 16'd0;
	end
	else
	begin
		if (z [3][15]) // 剩余角度为负数，顺时针旋转，d=-1
		begin
			x [4] <= x [3] + (y [3] >>> 3);
			y [4] <= y [3] - (x [3] >>> 3);
			z [4] <= z [3] + arctan_3;
		end          // 剩余角度为正数，逆时针旋转，d=+1
		else
		begin
			x [4] <= x [3] - (y [3] >>> 3);
			y [4] <= y [3] + (x [3] >>> 3);
			z [4] <= z [3] - arctan_3;
		end
	end
end


always@(posedge clk or negedge rst_n) //i=4
begin
	if (!rst_n)
	begin
		x [5] <= 16'd0;
		y [5] <= 16'd0;
		z [5] <= 16'd0;
	end
	else
	begin
		if (z [4][15]) // 剩余角度为负数，顺时针旋转，d=-1
		begin
			x [5] <= x [4] + (y [4] >>> 4);
			y [5] <= y [4] - (x [4] >>> 4);
			z [5] <= z [4] + arctan_4;
		end          // 剩余角度为正数，逆时针旋转，d=+1
		else
		begin
			x [5] <= x [4] - (y [4] >>> 4);
			y [5] <= y [4] + (x [4] >>> 4);
			z [5] <= z [4] - arctan_4;
		end
	end
end

always@(posedge clk or negedge rst_n) //i=5
begin
	if (!rst_n)
	begin
		x [6] <= 16'd0;
		y [6] <= 16'd0;
		z [6] <= 16'd0;
	end
	else
	begin
		if (z [5][15]) // 剩余角度为负数，顺时针旋转，d=-1
		begin
			x [6] <= x [5] + (y [5] >>> 5);
			y [6] <= y [5] - (x [5] >>> 5);
			z [6] <= z [5] + arctan_5;
		end          // 剩余角度为正数，逆时针旋转，d=+1
		else
		begin
			x [6] <= x [5] - (y [5] >>> 5);
			y [6] <= y [5] + (x [5] >>> 5);
			z [6] <= z [5] - arctan_5;
		end
	end
end

always@(posedge clk or negedge rst_n) //i=6
begin
	if (!rst_n)
	begin
		x [7] <= 16'd0;
		y [7] <= 16'd0;
		z [7] <= 16'd0;
	end
	else
	begin
		if (z [6][15]) // 剩余角度为负数，顺时针旋转，d=-1
		begin
			x [7] <= x [6] + (y [6] >>> 6);
			y [7] <= y [6] - (x [6] >>> 6);
			z [7] <= z [6] + arctan_6;
		end          // 剩余角度为正数，逆时针旋转，d=+1
		else
		begin
			x [7] <= x [6] - (y [6] >>> 6);
			y [7] <= y [6] + (x [6] >>> 6);
			z [7] <= z [6] - arctan_6;
		end
	end
end


always@(posedge clk or negedge rst_n) //i=7
begin
	if (!rst_n)
	begin
		x [8] <= 16'd0;
		y [8] <= 16'd0;
		z [8] <= 16'd0;
	end
	else
	begin
		if (z [7][15]) // 剩余角度为负数，顺时针旋转，d=-1
		begin
			x [8] <= x [7] + (y [7] >>> 7);
			y [8] <= y [7] - (x [7] >>> 7);
			z [8] <= z [7] + arctan_7;
		end          // 剩余角度为正数，逆时针旋转，d=+1
		else
		begin
			x [8] <= x [7] - (y [7] >>> 7);
			y [8] <= y [7] + (x [7] >>> 7);
			z [8] <= z [7] - arctan_7;
		end
	end
end

always@(posedge clk or negedge rst_n) //i=8
begin
	if (!rst_n)
	begin
		x [9] <= 16'd0;
		y [9] <= 16'd0;
		z [9] <= 16'd0;
	end
	else
	begin
		if (z [8][15]) // 剩余角度为负数，顺时针旋转，d=-1
		begin
			x [9] <= x [8] + (y [8] >>> 8);
			y [9] <= y [8] - (x [8] >>> 8);
			z [9] <= z [8] + arctan_8;
		end          // 剩余角度为正数，逆时针旋转，d=+1
		else
		begin
			x [9] <= x [8] - (y [8] >>> 8);
			y [9] <= y [8] + (x [8] >>> 8);
			z [9] <= z [8] - arctan_8;
		end
	end
end

always@(posedge clk or negedge rst_n) //i=9
begin
	if (!rst_n)
	begin
		x [10] <= 16'd0;
		y [10] <= 16'd0;
		z [10] <= 16'd0;
	end
	else
	begin
		if (z [9][15]) // 剩余角度为负数，顺时针旋转，d=-1
		begin
			x [10] <= x [9] + (y [9] >>> 9);
			y [10] <= y [9] - (x [9] >>> 9);
			z [10] <= z [9] + arctan_9;
		end          // 剩余角度为正数，逆时针旋转，d=+1
		else
		begin
			x [10] <= x [9] - (y [9] >>> 9);
			y [10] <= y [9] + (x [9] >>> 9);
			z [10] <= z [9] - arctan_9;
		end
	end
end

always@(posedge clk or negedge rst_n) //i=10
begin
	if (!rst_n)
	begin
		x [11] <= 16'd0;
		y [11] <= 16'd0;
		z [11] <= 16'd0;
	end
	else
	begin
		if (z [10][15]) // 剩余角度为负数，顺时针旋转，d=-1
		begin
			x [11] <= x [10] + (y [10] >>> 10);
			y [11] <= y [10] - (x [10] >>> 10);
			z [11] <= z [10] + arctan_10;
		end          // 剩余角度为正数，逆时针旋转，d=+1
		else
		begin
			x [11] <= x [10] - (y [10] >>> 10);
			y [11] <= y [10] + (x [10] >>> 10);
			z [11] <= z [10] - arctan_10;
		end
	end
end

// 溢出判断
assign x_tmp = x [11][15]==1 ? 16'h7FFF : x [11];
assign y_tmp = y [11][15]==1 ? 16'h7FFF : y [11];
//assign x_tmp = x [11];
//assign y_tmp = y [11];

always@(posedge clk or negedge rst_n) //i=11
begin
	if (!rst_n)
	begin
		sin_theta <= 16'd0;
		cos_theta <= 16'd0;
	end
	else
	begin
		if (Quadrant == 3'd1)
		begin
			sin_theta <= y_tmp;
			cos_theta <= x_tmp;
		end
		else if (Quadrant == 3'd2)
		begin
			sin_theta <= y_tmp;
			cos_theta <= -x_tmp;
		end
		else if (Quadrant == 3'd3)
		begin
			sin_theta <= -y_tmp;
			cos_theta <= -x_tmp;
		end
		else if (Quadrant == 3'd4)
		begin
			sin_theta <= -y_tmp;
			cos_theta <= x_tmp;
		end
		else
		begin
			sin_theta <= sin_theta;
			cos_theta <= cos_theta;
		end
	end
end
endmodule

testbench 文件

`timescale 1 ns/ 1 ns

module MyCORDIC_tb (
);

integer i;

reg clk,rst_n;
reg [15:0] theta,theta_n;
wire [15:0] sin_theta,cos_theta;

reg [15:0] cnt,cnt_n;

MyCORDIC u0 (
	.clk       (clk      ),
	.rst_n     (rst_n    ),
	.theta     (theta    ),
	.sin_theta (sin_theta),
	.cos_theta (cos_theta)
);

initial
begin
    #0 clk = 1'b0;
    #10 rst_n = 1'b0;
    #10 rst_n = 1'b1;
    #10000000$stop;
end 

initial
begin
    #0 theta = 16'd20;
   for (i=0;i<10000;i=i+1)
	begin
		#400;
		theta <= theta + 16'd200;
	end
end 

always #10
begin
    clk = ~clk;
end

endmodule

ModelSim 仿真

img

参考文献

  1. XILINX FPGAs for DSP CORDIC 算法 [OL].
  2. 张剑锋。基于改进 CORDIC 算法的 DDFS 和 FFT 研究与实现 [D].
    国防科学技术大学，2011.
  3. Cordic 算法的原理介绍 [OL].
    https://www.cnblogs.com/touchblue/p/3535968.html
  4. cordic 算法详解 [OL].
    https://blog.csdn.net/u010712012/article/details/77755567

via:

  • Cordic 算法的原理介绍 - 乐富道 - 博客园 乐富道 2014-01-28 23:05
    https://www.cnblogs.com/touchblue/p/3535968.html

  • 基于 FPGA 的 CORDIC 算法实现 ——Verilog 版_cordic 算法详解 一休哥 - CSDN 博客 善良的一休君于 2017-08-21 20:07:34 发布.
    https://blog.csdn.net/qq_39210023/article/details/77456031

  • CORDIC 算法原理详解及其 Verilog 实现_verilog cordic 算法 - CSDN 博客 爱吃蛋挞的Dolly
    于 2021-01-11 16:54:10 发布
    https://blog.csdn.net/gemengxia/article/details/112475073

  • CORDIC 算法详解及 FPGA 实现 - CSDN 博客 ngany 于 2021-05-30 18:52:04 发布
    https://blog.csdn.net/ngany/article/details/117401494

本文来自互联网用户投稿，该文观点仅代表作者本人，不代表本站立场。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如若转载，请注明出处：http://www.coloradmin.cn/o/2262898.html

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈，一经查实，立即删除！

相关文章

鸿蒙学习笔记：用户登录界面

鸿蒙学习笔记：用户登录界面

文章目录 1. 提出任务2. 完成任务2.1 创建鸿蒙项目2.2 准备图片资源2.3 编写首页代码2.4 启动应用 3. 实战小结 1. 提出任务 本次任务聚焦于运用 ArkUI 打造用户登录界面。需呈现特定元素&#xff1a;一张图片增添视觉感&#xff0c;两个分别用于账号与密码的文本输入框&#…
阅读更多...
着色器 （三）

着色器 （三）

今天&#xff0c;是我们介绍opengl着色器最后一章&#xff0c;着色器(Shader)是运行在GPU上的小程序。这些小程序为图形渲染管线的某个特定部分而运行。从基本意义上来说&#xff0c;着色器只是一种把输入转化为输出的程序。着色器也是一种非常独立的程序&#xff0c;因为它们之…
阅读更多...
leetcode:3285. 找到稳定山的下标（python3解法）

leetcode:3285. 找到稳定山的下标（python3解法）

难度&#xff1a;简单 有 n 座山排成一列&#xff0c;每座山都有一个高度。给你一个整数数组 height &#xff0c;其中 height[i] 表示第 i 座山的高度&#xff0c;再给你一个整数 threshold 。 对于下标不为 0 的一座山&#xff0c;如果它左侧相邻的山的高度 严格大于 thresho…
阅读更多...
深度学习之超分辨率算法——SRGAN

深度学习之超分辨率算法——SRGAN

更新版本 实现了生成对抗网络在超分辨率上的使用 更新了损失函数&#xff0c;增加先验函数 SRresnet实现 import torch import torchvision from torch import nnclass ConvBlock(nn.Module):def __init__(self, kernel_size3, stride1, n_inchannels64):super(ConvBlock…
阅读更多...
集成方案 | Docusign + 金蝶云，实现合同签署流程自动化！

集成方案 | Docusign + 金蝶云，实现合同签署流程自动化！

本文将详细介绍 Docusign 与金蝶云的集成步骤及其效果&#xff0c;并通过实际应用场景来展示 Docusign 的强大集成能力&#xff0c;以证明 Docusign 集成功能的高效性和实用性。 在当今商业环境中&#xff0c;流程的无缝整合与数据的实时性对于企业的成功至关重要。金蝶云&…
阅读更多...
数据结构----链表头插中插尾插

数据结构----链表头插中插尾插

一、链表的基本概念 链表是一种线性数据结构&#xff0c;它由一系列节点组成。每个节点包含两个主要部分&#xff1a; 数据域&#xff1a;用于存储数据元素&#xff0c;可以是任何类型的数据&#xff0c;如整数、字符、结构体等。指针域&#xff1a;用于存储下一个节点&#…
阅读更多...
Service Discovery in Microservices 客户端/服务端服务发现

Service Discovery in Microservices 客户端/服务端服务发现

原文链接 Client Side Service Discovery in Microservices - GeeksforGeeks 原文链接 Server Side Service Discovery in Microservices - GeeksforGeeks 目录 服务发现介绍 Server-Side 服务发现 实例&#xff1a; Client-Side 服务发现 实例&#xff1a; 服务发现介绍…
阅读更多...
Git连接远程仓库（超详细）

Git连接远程仓库（超详细）

目录 一、Gitee 远程仓库连接 1. HTTPS 方式 2. SSH公钥方式 &#xff08;1&#xff09;账户公钥 &#xff08;2&#xff09;仓库公钥 仓库的 SSH Key 和账户 SSH Key 的区别&#xff1f;​ 二、GitHub远程仓库连接 1. HTTPS方式 2.SSH公钥方式 本文将介绍如何通过 H…
阅读更多...
系列4:基于Centos-8.6 Kubernetes多网卡节点Calico选择网卡配置

系列4:基于Centos-8.6 Kubernetes多网卡节点Calico选择网卡配置

每日禅语 不动心”是一个人修养和定力的体现&#xff0c;若一个人心无定力&#xff0c;就会被外界环境左右&#xff0c;随外界的境遇而动摇。佛家认为&#xff0c;心是一切的基础&#xff0c;一个人如果想要真正入定&#xff0c;必须先从修心开始。修心即是净心&#xff0c;心灵…
阅读更多...
Docker：Dockerfile(补充四)

Docker：Dockerfile(补充四)

这里写目录标题 1. Dockerfile常见指令1.1 DockerFile例子 2. 一些其他命令 1. Dockerfile常见指令 简单的dockerFile文件 FROM openjdk:17LABEL authorleifengyangCOPY app.jar /app.jarEXPOSE 8080ENTRYPOINT ["java","-jar","/app.jar"]# 使…
阅读更多...
98. 验证二叉搜索树(java)

98. 验证二叉搜索树(java)

题目描述&#xff1a; 给你一个二叉树的根节点 root &#xff0c;判断其是否是一个有效的二叉搜索树。 有效 二叉搜索树定义如下&#xff1a; 节点的左 子树 只包含 小于 当前节点的数。节点的右子树只包含 大于 当前节点的数。所有左子树和右子树自身必须也是二叉搜索树。 …
阅读更多...
微软 Phi-4：小型模型的推理能力大突破

微软 Phi-4：小型模型的推理能力大突破

在人工智能领域&#xff0c;语言模型的发展日新月异。微软作为行业的重要参与者&#xff0c;一直致力于推动语言模型技术的进步。近日&#xff0c;微软推出了最新的小型语言模型 Phi-4&#xff0c;这款模型以其卓越的复杂推理能力和在数学领域的出色表现&#xff0c;引起了广泛…
阅读更多...
libaom 源码分析：熵编码模块介绍

libaom 源码分析：熵编码模块介绍

AV1 熵编码原理介绍 关于AV1 熵编码原理介绍可以参考：AV1 编码标准熵编码技术概述libaom 熵编码相关源码介绍 函数流程图 核心函数介绍 av1_pack_bitstream 函数：该函数负责将编码后的数据打包成符合 AV1 标准的比特流格式；包括写入序列头 OBU 的函数 av1_write_obu_header…
阅读更多...
JAVA基于百度AI人脸识别签到考勤系统（开题报告+作品+论文）

JAVA基于百度AI人脸识别签到考勤系统（开题报告+作品+论文）

博主介绍&#xff1a;黄菊华老师《Vue.js入门与商城开发实战》《微信小程序商城开发》图书作者&#xff0c;CSDN博客专家&#xff0c;在线教育专家&#xff0c;CSDN钻石讲师&#xff1b;专注大学生毕业设计教育、辅导。 所有项目都配有从入门到精通的基础知识视频课程&#xff…
阅读更多...
go 中使用redis 基础用法

go 中使用redis 基础用法

1、安装redis 参考链接&#xff1a;https://www.codeleading.com/article/98554130215/ 1.1 查看是否有redis yum 源 yum install redis没有可用的软件包&#xff0c;执行1.2 1.2下载fedora的epel仓库 yum install epel-release --下载fedora的epel仓库1.3启动redis s…
阅读更多...
postman添加cookie

postman添加cookie

点击cookies 输入域名&#xff0c;添加该域名下的cookies 发送改域名下的请求&#xff0c;cookie会自动追加上
阅读更多...
简易记事本开发-（SSM+Vue）

简易记事本开发-（SSM+Vue）

目录 前言 一、项目需求分析 二、项目环境搭建 1.创建MavenWeb项目&#xff1a; 2.配置 Spring、SpringMVC 和 MyBatis SpringMVC 配置文件 (spring-mvc.xml)&#xff1a; 配置视图解析器、处理器映射器&#xff0c;配置了CORS&#xff08;跨源资源共享&#xff09;&#x…
阅读更多...
vsCode 报错[vue/no-v-model-argument]e‘v-model‘ directives require no argument

vsCode 报错[vue/no-v-model-argument]e‘v-model‘ directives require no argument

在vue3中使用ui库中的组件语法v-model:value时会提示[vue/no-multiple-template-root]The template root requires exactly one element. 引入组件使用单标签时会提示[vue/no-multiple-template-root]“The template root requires exactly one element. 原因&#xff1a; 1.可…
阅读更多...
初学stm32 -- SysTick定时器

初学stm32 -- SysTick定时器

以delay延时函数来介绍SysTick定时器的配置与使用 首先是delay_init()延时初始化函数&#xff0c;这个函数主要是去初始化SysTick定时器&#xff1b; void delay_init() {SysTick_CLKSourceConfig(SysTick_CLKSource_HCLK_Div8); //选择外部时钟 HCLK/8fac_usSystemCoreCloc…
阅读更多...
Gitlab 数据备份全攻略：命令、方法与注意事项

Gitlab 数据备份全攻略：命令、方法与注意事项

文章目录 1、备份命令2、备份目录名称说明3、手工备份配置文件3.1 备份配置文件3.2 备份ssh文件 4、备份注意事项4.1 停止puma和sicdekiq组件4.2 copy策略需要更多磁盘空间 5、数据备份方法5.1 docker命令备份5.2 kubectl命令备份5.3 参数说明5.4、选择性备份5.5、非tar备份5.6…
阅读更多...
推荐文章
最新文章

Copyright @ 2022~2023