当前位置:   article > 正文

基于FPGA的CORDIC算法实现——Verilog版_cordic算法详解 一休哥

cordic算法详解 一休哥

目前,学习与开发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的几何原理介绍

假设在xy坐标系中有一个点P1(x1,y1),将P1点绕原点旋转θ角后得到点P2(x2,y2)。
这里写图片描述
于是可以得到P1和P2的关系。
x2 = x1cosθ – y1sinθ = cosθ(x1 – y1tanθ)
y2 = y1cosθ + x1sinθ = cosθ(y1 +x1tanθ)
以上就是CORDIC的几何原理部分,而我们该如何深入理解这个几何原理的真正含义呢?
从原理中,我们可以知道,当已知一个点P1的坐标,并已知该点P1旋转的角度θ,则可以根据上述公式求得目标点P2的坐标。然后,麻烦来了,我们需要用FPGA去执行上述运算操作,而FPGA的Verilog语言根本不支持三角函数运算。因此,我们需要对上述式子进行简化操作,将复杂的运算操作转换为一种单一的“迭代位移”算法。那么,接下来我们介绍优化算法部分。

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θ1 tanθ0)
y2 = cosθ1(y1 + d1x1tanθ1) = cosθ1cosθ0(y0 + d0x0tanθ0 + d1x0tanθ1 – d1d0y0tanθ1 tanθ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 – d iy i2-i,提取了cosθi,2-i等效替换了tanθi之后
yi+1 = yi + d ix i2-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元乘积的倒数
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)
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24
  • 25
  • 26

这里写图片描述
从上表也可以看出,当迭代次数为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的值我们也同样预先存入系统中。如此,我们就实现了正弦和余弦操作了。

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

1、CORDIC的旋转模式

由于算法较复杂,一休哥再总结一些具体的操作流程。
1、 设置迭代次数为16,则x0 = 0.607253,y0 = 0,并输入待计算的角度θ,θ在[-99.7°,99.7°]范围内。
2、 根据三个迭代公式进行迭代,i从0至15:
xi+1 = xi – d iy i2-i
yi+1 = yi + d ix i2-i
zi+1 = zi - diθi
注:z0 = θ,di与zi同符号。
3、 经过16次迭代计算后,得到的x16 和y16分别为cosθ和sinθ。
至此,关于CORDIC的三角函数cosθ和sinθ的计算原理讲解结束。
关于CORDIC算法计算三角函数cosθ和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)
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24

2、CORDIC的向量模式

讲完了旋转模式后,我们接着讲讲向量模式下的圆坐标系。
在这里,我们需从头来过了,假设在xy坐标系中有一个点P0(x0,y0)&

声明:本文内容由网友自发贡献,不代表【wpsshop博客】立场,版权归原作者所有,本站不承担相应法律责任。如您发现有侵权的内容,请联系我们。转载请注明出处:https://www.wpsshop.cn/w/盐析白兔/article/detail/625447
推荐阅读
相关标签
  

闽ICP备14008679号