公式渲染引擎由Mathjax换为Katex,文章可能存在渲染问题,如有发现问题请留言。

前世今生

以下源自参考资料[6]:CORDIC技术并不是什么新鲜的东西。事实上它可以追溯到1957年由J.Volder发表的一篇文章。在上个世纪五十年代,在大型实际的计算机中的实行移位相加受到了当时技术上的限制,所以使用CORDIC变得非常必要。到了七十年代,Hewlett Packard和其他公司出产了手持计算器,许多计算器使用一个内部CORDIC单元来计算所有的三角函数(了解这件事的人们一定还记得,那时求一个角度的正切值需要延迟大约1秒中)。二十世纪八十年代,随着高速度乘法器与带有大存储量的通用处理器的出现,CORDIC算法变得无关紧要了。然而在二十一世纪的今天,对于FPGA来说,CORDIC一定是在DSP应用中(诸如多输入多输出(MIMO),波束形成以及其他自适应系统)计算三角函数的备选技术。

本文将详细介绍CORDIC算法,所有参考资料见文章尾部。

Cordic详解

圆周坐标系

从坐标旋转开始

如下图,一个直角坐标点(x1,y1)(x_1,y_1)逆时针旋转θθ角度到点(x2,y2)(x_2,y_2)如何计算(x2,y2)(x_2,y_2)呢?

cordic1

公式如下:

x2=x1cosθy1cosθy2=x1sinθ+y1cosθx_2 = x_1cosθ-y_1cosθ \\ y_2 = x_1sinθ+y_1cosθ

其实只需要一点高中的知识就能推导这个公式了。

cordic2

假设底角为αα,旋转半径为11,则

x2=cos(α+θ)=cos(α)cos(θ)sin(α)sin(θ)=x1cosθy1cosθy2=sin(α+θ)=sin(α)cos(θ)+cos(α)sin(θ)=x1sinθ+y1cosθx_2=cos(α+θ) = cos(α)cos(θ) - sin(α)sin(θ) =x_1cosθ-y_1cosθ \\ y_2=sin(α+θ) = sin(α)cos(θ) + cos(α)sin(θ) =x_1sinθ+y_1cosθ

上面的公式也可以写成矩阵的形式,即

[x2y2]=[cosθsinθsinθcosθ][x1y1]\left[ \begin{matrix} x_2\\ y_2 \end{matrix} \right] = \left[ \begin{matrix} cosθ&-sinθ\\ sinθ&cosθ \end{matrix} \right] \left[ \begin{matrix} x_1\\ y_1 \end{matrix} \right]

以上公式对应逆时针旋转。若是顺时针旋转,则公式应为

[x1y1]=[cosθsinθsinθcosθ][x2y2]\left[ \begin{matrix} x_1\\ y_1 \end{matrix} \right] = \left[ \begin{matrix} cosθ&sinθ\\ -sinθ&cosθ \end{matrix} \right] \left[ \begin{matrix} x_2\\ y_2 \end{matrix} \right]

那么已知旋转角θθ,如何通过FPGA来计算x2x_2y2y_2呢?答案就是反复迭代

二分角度旋转

假设(x1,y1)=(100,200)(x_1,y_1)=(100,200),要求其极坐标系下的坐标(ρ,θ)(ρ,θ)。当然,求θθ的过程也就是求arctan(y/x)arctan(y/x)的过程。首先通过计算器得到结果(ρ,θ)=(223.61,63.435)(ρ,θ) = (223.61,63.435)。下面计算先只关注θθ的旋转变化,不关注ρρ的长度变化。这时最直观简单的想法就是先旋转一个角度进行一次尝试。

已知角度在(0°,90°)(0°,90°)这个范围内,借鉴二分法的思想,先顺时针旋转45°45°进行尝试。根据之前推导的公式可知

[xy]=[cos(45°)sin(45°)sin(45°)cos(45°)][100200]=[212.1370.711]\left[ \begin{matrix} x\\ y \end{matrix} \right] = \left[ \begin{matrix} cos(45°)&sin(45°)\\ -sin(45°)&cos(45°) \end{matrix} \right] \left[ \begin{matrix} 100\\ 200 \end{matrix} \right] = \left[ \begin{matrix} 212.13\\ 70.711 \end{matrix} \right]

发现纵坐标y>0y>0,则旋转的角度不够,则在此基础上继续顺时针旋转45/2=22.5°45/2=22.5°,此时角度为45+22.5=67.5°45+22.5=67.5°。根据推导的公式可知

[xy]=[cos(22.5°)sin(22.5°)sin(22.5°)cos(22.5°)][212.1370.711]=[223.0415.85]\left[ \begin{matrix} x\\ y \end{matrix} \right] = \left[ \begin{matrix} cos(22.5°)&sin(22.5°)\\ -sin(22.5°)&cos(22.5°) \end{matrix} \right] \left[ \begin{matrix} 212.13\\ 70.711 \end{matrix} \right] = \left[ \begin{matrix} 223.04\\ -15.85 \end{matrix} \right]

发现纵坐标y<0y<0,则旋转角度超过了范围,则在此基础上逆时针回转22.5/2=11.25°22.5/2=11.25°,此时角度为45+22.511.25=56.25°45+22.5-11.25=56.25°。根据推导的公式可知

[xy]=[cos(11.25°)sin(11.25°)sin(11.25°)cos(11.25°)][223.0415.85]=[221.8527.967]\left[ \begin{matrix} x\\ y \end{matrix} \right] = \left[ \begin{matrix} cos(11.25°)&-sin(11.25°)\\ sin(11.25°)&cos(11.25°) \end{matrix} \right] \left[ \begin{matrix} 223.04\\ -15.85 \end{matrix} \right] = \left[ \begin{matrix} 221.85\\ 27.967 \end{matrix} \right]

发现纵坐标y>0y>0,则旋转的角度又不够,则此基础上继续顺时针旋转11.25/2=5.625°11.25/2=5.625°,此时角度为45+22.511.25+5.625=61.875°45+22.5-11.25+5.625=61.875°。根据推导的公式可知

[xy]=[cos(5.625°)sin(5.625°)sin(5.625°)cos(5.625°)][221.8527.967]=[223.526.0874]\left[ \begin{matrix} x\\ y \end{matrix} \right] = \left[ \begin{matrix} cos(5.625°)&sin(5.625°)\\ -sin(5.625°)&cos(5.625°) \end{matrix} \right] \left[ \begin{matrix} 221.85\\ 27.967 \end{matrix} \right] = \left[ \begin{matrix} 223.52\\ 6.0874 \end{matrix} \right]

这时纵坐标已经比较接近00了。按照这个模式继续算下去,可以不断逼近真实的结果。比如现在的结果范围为(ρ,θ)=(223.52,61.875±5.625)(ρ,θ)=(223.52,61.875±5.625),已经接近实际结果(ρ,θ)=(223.61,63.435)(ρ,θ)=(223.61,63.435)了。旋转过程演示图如下。

rotate

因此通过多次迭代,可以逐渐逼近结果。但是如何计算上面的cos(45),sin(45),cos(22.5),sin(22.5),cos(11.25),sin(11.25)cos(45),sin(45),cos(22.5),sin(22.5),cos(11.25),sin(11.25)等值?其实如果仔细观察,这些点都是固定的,因此他们的cos()cos()sin()sin()也是固定的,所以只需要把要用的点的sin()sin()cos()cos()提前计算好存起来,用时查表即可。Python实现如下。

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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
import math
import numpy as np

ITERATION_TIMES = 64
i = 0
d = 45
cos = [0] * ITERATION_TIMES
sin = [0] * ITERATION_TIMES
degree = [0] * ITERATION_TIMES

x = 100
y = 200
x_n = 0
y_n = 0
z = 0

#查找表
for i in range(ITERATION_TIMES):
cos[i] = np.cos(d*np.pi/180)
sin[i] = np.sin(d*np.pi/180)
degree[i] = d;
d = d/2

#迭代计算
for i in range(ITERATION_TIMES):
if(y>0):
x_n = x * cos[i] + y * sin[i]
y_n = y * cos[i] - x * sin[i]
z = z + degree[i]
x = x_n
y = y_n
#print("rotate_angle=",degree[i],"y=",y_n,"z=",z)
else:
x_n = x * cos[i] - y * sin[i]
y_n = y * cos[i] + x * sin[i]
z = z - degree[i]
x = x_n
y = y_n
#print("rotate_angle=",degree[i],"y=",y_n,"z=",z)
print(z)

经过不同迭代次数,得到不同的结果。迭代次数越高,精度自然就越高。

迭代次数 计算角度(实际结果63.43494882292201)
16 63.433685302734375
32 63.43494881177321
64 63.43494882292201

长度缩放旋转

二分角度旋转的计算存在一个问题:计算量巨大。可以发现每一次旋转变换都需要进行4次浮点乘法运算,对于FPGA来说这个运算量过大。改进的方法就是通过变换坐标旋转的公式。

每次的旋转公式可以总结为

[xi+1yi+1]=[cos(θi)disin(θi)disin(θi)cos(θi)][xiyi]=cos(θi)[1ditan(θi)ditan(θi)1][xiyi]zi+1={zi+θi,if yi>0ziθi,if yi≤0where di={+1,if yi>01,if yi≤0 ,θi=45°2i\left[ \begin{matrix} x_{i+1}\\ y_{i+1} \end{matrix} \right] = \left[ \begin{matrix} cos(θ_i)&d_isin(θ_i)\\ -d_isin(θ_i)&cos(θ_i) \end{matrix} \right] \left[ \begin{matrix} x_i\\ y_i \end{matrix} \right] = cos(θ_i) \left[ \begin{matrix} 1&d_itan(θ_i)\\ -d_itan(θ_i)&1 \end{matrix} \right] \left[ \begin{matrix} x_i\\ y_i \end{matrix} \right] \\ z_{i+1} = \begin{cases} z_i + θ_i, & \text{if $y_i$>0} \\ z_i - θ_i, & \text{if $y_i$≤0} \end{cases} \\ where\ d_i = \begin{cases} +1, & \text{if $y_i$>0} \\ -1, & \text{if $y_i$≤0} \end{cases}\ ,θ_i=\frac{45°}{2^i}

这样进行转换后可以cos(θ)cos(θ)从矩阵运算中提取出来并省略掉。为什么可以将cos(θ)cos(θ)省略掉呢? 因为现在的计算只求θθ的值,即旋转角度变化,不关心ρρ长度的变化,省略掉的cos(θ)cos(θ)只是在长度上进行了放缩,并不影响其旋转的角度。因此转换公式变为

[xi+1yi+1]=[1ditan(θi)ditan(θi)1][xiyi]where di={+1,if yi>01,if yi≤0\left[ \begin{matrix} x_{i+1}\\ y_{i+1} \end{matrix} \right] = \left[ \begin{matrix} 1&d_itan(θ_i)\\ -d_itan(θ_i)&1 \end{matrix} \right] \left[ \begin{matrix} x_i\\ y_i \end{matrix} \right] \\ where\ d_i = \begin{cases} +1, & \text{if $y_i$>0} \\ -1, & \text{if $y_i$≤0} \end{cases}

这种旋转方式被称为伪旋转(Pseudo-Rotation)。这样乘法的个数就从每一次旋转变换的4次变为了两次,存查找表的值也只需要tan(45),tan(22.5)tan(45),tan(22.5)…这样的值,减少了一半。过程演示图如下。

rotate

由上方过程演示图对比可以看出,旋转向量的长度在增长,但其旋转角度变化保持与之前一致。Python实现如下。

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
27
28
29
30
31
32
33
34
35
36
37
38
import math
import numpy as np

ITERATION_TIMES = 64
i = 0
d = 45
tan = [0] * ITERATION_TIMES
degree = [0] * ITERATION_TIMES

x = 100
y = 200
x_n = 0
y_n = 0
z = 0

#查找表
for i in range(ITERATION_TIMES):
tan[i] = np.tan(d*np.pi/180)
degree[i] = d;
d = d/2

#迭代计算
for i in range(ITERATION_TIMES):
if(y>0):
x_n = x + y * tan[i]
y_n = y - x * tan[i]
z = z + degree[i]
x = x_n
y = y_n
#print("rotate_angle=",degree[i],"y=",y_n,"z=",z)
else:
x_n = x - y * tan[i]
y_n = y + x * tan[i]
z = z - degree[i]
x = x_n
y = y_n
#print("rotate_angle=",degree[i],"y=",y_n,"z=",z)
print(z)

结果与之前的二分角度旋转一致。

cordic旋转

经过了长度缩放的变换,算法成功将乘法运算减少了一半。那还能不能继续减少运算量呢?其实依旧可以从之前的长度缩放变换的转换公式入手。

[xi+1yi+1]=[1ditan(θi)ditan(θi)1][xiyi]where di={+1,if yi>01,if yi≤0\left[ \begin{matrix} x_{i+1}\\ y_{i+1} \end{matrix} \right] = \left[ \begin{matrix} 1&d_itan(θ_i)\\ -d_itan(θ_i)&1 \end{matrix} \right] \left[ \begin{matrix} x_i\\ y_i \end{matrix} \right] where\ d_i= \begin{cases} +1, & \text{if $y_i$>0} \\ -1, & \text{if $y_i$≤0} \end{cases}

第一次旋转45°时,转换公式为

[xi+1yi+1]=[1tan(45)tan(45)1][xiyi]=[1111][xiyi]\left[ \begin{matrix} x_{i+1}\\ y_{i+1} \end{matrix} \right] = \left[ \begin{matrix} 1&tan(45)\\ -tan(45)&1 \end{matrix} \right] \left[ \begin{matrix} x_i\\ y_i \end{matrix} \right] = \left[ \begin{matrix} 1&1\\ -1&1 \end{matrix} \right] \left[ \begin{matrix} x_i\\ y_i \end{matrix} \right]

可见此次转换其实不需要进行乘法运算,可以直接进行加减运算就可以得到结果。然而第二次旋转22.5°22.5°的时候,由于tan(22.5°)=0.4142135623731tan(22.5°)=0.4142135623731,这是一个不整的小数。因此需要复杂的乘法运算才能得到结果。

那有没有办法规避这种复杂乘法?当然是有的。注意一开始的二分法是通过对45°45°角进行二分得到了22.5°22.5°,但其实迭代旋转角度可以不选择22.5°22.5°,只要比45°45°角小就可以。假设选择tan(θ)=12tan(θ)=\frac{1}{2}θθ角(θ=26.565051177078°θ=26.565051177078°),则第二次旋转的公式变换为

[xi+1yi+1]=[1tan(26.565051177078°)tan(26.565051177078°)1][xiyi]=[112121][xiyi]\left[ \begin{matrix} x_{i+1}\\ y_{i+1} \end{matrix} \right] = \left[ \begin{matrix} 1&tan(26.565051177078°)\\ -tan(26.565051177078°)&1 \end{matrix} \right] \left[ \begin{matrix} x_i\\ y_i \end{matrix} \right] = \left[ \begin{matrix} 1&\frac{1}{2}\\ -\frac{1}{2}&1 \end{matrix} \right] \left[ \begin{matrix} x_i\\ y_i \end{matrix} \right]

同样的,每次选择tan(θ)=1,12,14tan(θ)=1,\frac{1}{2},\frac{1}{4}…依次类推,这样计算就在定点运算层面简化了。因为乘12\frac{1}{2}的操作等于将数右移一位,乘14\frac{1}{4}等于将数右移两位,以此类推。这样,复杂的乘法运算转化为了简单的移位和加减运算。通项转换公式变换为

[xi+1yi+1]=[1ditan(θi)ditan(θi)1][xiyi]=[1di2idi2i1][xiyi]zi+1={zi+θi,if yi>0ziθi,if yi≤0where di={+1,if yi>01,if yi≤0 ,θi=arctan(12i)\left[ \begin{matrix} x_{i+1}\\ y_{i+1} \end{matrix} \right] = \left[ \begin{matrix} 1&d_itan(θ_i)\\ -d_itan(θ_i)&1 \end{matrix} \right] \left[ \begin{matrix} x_i\\ y_i \end{matrix} \right] = \left[ \begin{matrix} 1&\frac{d_i}{2^i}\\ -\frac{d_i}{2^i}&1 \end{matrix} \right] \left[ \begin{matrix} x_i\\ y_i \end{matrix} \right] \\ z_{i+1} = \begin{cases} z_i + θ_i, & \text{if $y_i$>0} \\ z_i - θ_i, & \text{if $y_i$≤0} \end{cases} \\ where\ d_i = \begin{cases} +1, & \text{if $y_i$>0} \\ -1, & \text{if $y_i$≤0} \end{cases}\ ,θ_i=arctan(\frac{1}{2^i})

同样的,可以将θ=arctan(12i)θ=arctan(\frac{1}{2^i})的值预先计算出来存起来,用时查表就可以。

tan(θ) θ
1 45.0°
0.5 26.56505117707799°
0.25 14.036243467926477°
0.125 7.125016348901798°
0.0625 3.5763343749973515°
‭0.03125‬ 1.7899106082460694°
‭0.015625‬ 0.8951737102110744°

过程演示图如下。

rotate

由上方过程演示图对比可以看出,cordic变换向量的长度与长度缩放旋转的变化类似,而其旋转角度发生了根本性变化,使得横纵坐标的变化的小数更整。Python实现如下。

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
27
28
29
30
31
32
33
34
35
36
37
import math
import numpy as np

ITERATION_TIMES = 64
i = 0
d = 1
degree = [0] * ITERATION_TIMES

x = 100
y = 200
x_n = 0
y_n = 0
z = 0

#查找表
for i in range(ITERATION_TIMES):
degree[i] = math.atan(d)
d = d/2

#迭代计算
d = 1
for i in range(ITERATION_TIMES):
if(y>0):
x_n = x + (y >> d)
y_n = y - (x >> d)
z = z + degree[i]
x = x_n
y = y_n
else:
x_n = x - (y >> d)
y_n = y + (x >> d)
z = z - degree[i]
x = x_n
y = y_n
d = d >> 1

print(z*180/np.pi)

至此基本讲清了cordic算法的本质原理(在圆周坐标下)。接下来补充一些其他cordic算法的相关内容。

参数、模式与坐标系

伸缩因子Kn

上一节中只关注了角度θθ的求解,并没有关注到长度ρρ的求解。实际上,如果考虑长度变化的话,那么cordic旋转公式应还原为(也就是将省略掉的cos(θ)cos(θ)补回来):

[xi+1yi+1]=[cos(θi)disin(θi)disin(θi)cos(θi)][xiyi]=cos(θi)[1ditan(θi)ditan(θi)1][xiyi]=cos(θi)[1di2idi2i1][xiyi]zi+1=zi+diθiwhere di={+1,if yi>01,if yi≤0, θi=arctan(12i)\left[ \begin{matrix} x_{i+1}\\ y_{i+1} \end{matrix} \right] = \left[ \begin{matrix} cos(θ_i)&d_isin(θ_i)\\ -d_isin(θ_i)&cos(θ_i) \end{matrix} \right] \left[ \begin{matrix} x_i\\ y_i \end{matrix} \right] = cos(θ_i) \left[ \begin{matrix} 1&d_itan(θ_i)\\ -d_itan(θ_i)&1 \end{matrix} \right] \left[ \begin{matrix} x_i\\ y_i \end{matrix} \right] = cos(θ_i) \left[ \begin{matrix} 1&\frac{d_i}{2^i}\\ -\frac{d_i}{2^i}&1 \end{matrix} \right] \left[ \begin{matrix} x_i\\ y_i \end{matrix} \right] \\ z_{i+1} = z_i + d_iθ_i\\ \\ where\ d_i = \begin{cases} +1, & \text{if $y_i$>0} \\ -1, & \text{if $y_i$≤0} \end{cases},\ θ_i = arctan(\frac{1}{2^i})

还原为方程组的形式:

xi+1=cos(θi)(xi+diyitan(θi)=cos(θi)(xi+diyi2i)yi+1=cos(θi)(yidixitan(θi)=cos(θi)(yidixi2i)zi+1=zi+diθiwhere di={+1,if yi>01,if yi≤0, θi=arctan(12i)x_{i+1} = cos(θ_i)(x_i+d_iy_itan(θ_i)=cos(θ_i)(x_i+d_i\frac{y_i}{2^i})\\ y_{i+1} = cos(θ_i)(y_i-d_ix_itan(θ_i)=cos(θ_i)(y_i-d_i\frac{x_i}{2^i})\\ z_{i+1} = z_i + d_iθ_i\\ where\ d_i = \begin{cases} +1, & \text{if $y_i$>0} \\ -1, & \text{if $y_i$≤0} \end{cases}, \ θ_i = arctan(\frac{1}{2^i})\\

假设进行64次迭代,则结果计算为

x64=n=063cos(θn)(xn+dyn2n)y64=n=063cos(θn)(yndxn2n)z64=z0+n=063dθnwhere{d=+1,if yn>0d=1,if yn≤0, θn=arctan(12n)x_{64} = \prod_{n=0}^{63}cos(θ_n)(x_n+d\frac{y_n}{2^n})\\ y_{64} = \prod_{n=0}^{63}cos(θ_n)(y_n-d\frac{x_n}{2^n})\\ z_{64} = z_0 + \sum_{n=0}^{63}dθ_n \\ where \begin{cases} d = +1, & \text{if $y_n$>0} \\ d = -1, & \text{if $y_n$≤0} \end{cases},\ θ_n = arctan(\frac{1}{2^n})\\

可以发现其实迭代过程的计算分为两部分,一部分是cordic的旋转变换,一部分是cordic的长度伸缩变换。

cordic长度伸缩变换:

n=063cos(θn)\prod_{n=0}^{63}cos(θ_n)

cordic旋转变换:

n=063(xn+dyn2n)n=063(yndxn2n)where{d=+1,if yn>0d=1,if yn≤0\prod_{n=0}^{63}(x_n+d\frac{y_n}{2^n})\\ \prod_{n=0}^{63}(y_n-d\frac{x_n}{2^n})\\ where \begin{cases} d = +1, & \text{if $y_n$>0} \\ d = -1, & \text{if $y_n$≤0} \end{cases}

由于在cordic旋转变换中,已知θn=arctan(12n)θ_n=arctan(\frac{1}{2^n}),因此长度伸缩变换可以变化为:

n=063(cos(θn))=n=063(11+tan2(θn))=n=063(11+(12n)2)\prod_{n=0}^{63}(cos(θ_n))=\prod_{n=0}^{63}(\sqrt{\frac{1}{1+tan^2(θ_n)}})=\prod_{n=0}^{63}(\sqrt{\frac{1}{1+(\frac{1}{2^n})^2}})

nn\to∞,上面的式子是收敛的,当结束变换时,应当乘以收敛值的倒数以恢复长度变换带来的影响。恢复长度变换影响而乘的这个数就叫做伸缩因子KnK_n,是一个固定常数。

Kn=limn(0n(1+(12n)2))1.646761Kn0.6073K_n=\lim_{n\to∞}(\prod_{0}^{n}(\sqrt{1+(\frac{1}{2^n})^2}))≈1.64676 \\ \frac{1}{K_n}≈ 0.6073

通过下表可见伸缩因子随迭代次数的变化

迭代次数 伸缩因子Kn
10 1.646759211139822
20 1.6467602581200669
50 1.6467602581210652
100 1.6467602581210652

可见伸缩因子随着迭代次数增加而变化越来越小,最终收敛。这样可以利用cordic旋转后通过乘以伸缩因子来保证长度变换也是正确的,从而求到长度ρρ。最终的变换公式如下。

xn+1=Knn=0n(xn+dyn2n)yn+1=Knn=0n(yndxn2n)zn+1=z0+n=0ndθnwhere{d=+1,if yn>0d=1,if yn≤0, θn=arctan(12n)x_{n+1} = K_n\prod_{n=0}^{n}(x_n+d\frac{y_n}{2^n})\\ y_{n+1} = K_n\prod_{n=0}^{n}(y_n-d\frac{x_n}{2^n})\\ z_{n+1} = z_0 + \sum_{n=0}^{n}dθ_n\\ where \begin{cases} d = +1, & \text{if $y_n$>0} \\ d = -1, & \text{if $y_n$≤0} \end{cases},\ θ_n = arctan(\frac{1}{2^n})\\

两种旋转模式

在圆周坐标系下的坐标旋转其实存在两个问题:

  1. 已知x1x_1y1y_1如何计算旋转角θθ使得y=0y=0

  2. 已知旋转角θθx1x_1,且y1=0y_1=0如何计算x2x_2y2y_2

这是两个镜像问题,只是从两个方向来考虑。

  1. 通过不断进行cordic旋转变换逼近y=0y=0的点(顺时针,逼近y=0y=0点)

    使结果纵坐标y0y\to0cordic旋转变换称为向量模式(vector mode)

  2. 通过不断进行cordic旋转变换逼近旋转角θθ(逆时针,逼近θ=0θ=0的点)

    使结果角度θ0θ\to0cordic旋转变换被称为旋转模式(rotation mode)

  • 圆周坐标的向量模式(vector mode)的公式为(上一节已经推导)
  1. 迭代式:

    {xi+1=xi+diyi2iyi+1=yidixi2izi+1=zi+dθiwhere{di=+1,if yi>0di=1,if yi≤0,θi=arctan(12i)\begin{cases} x_{i+1} = x_i+d_i\frac{y_i}{2^i}\\ y_{i+1} = y_i-d_i\frac{x_i}{2^i}\\ z_{i+1} = z_i + dθ_i\\ \end{cases} \\ where \begin{cases} d_i = +1, & \text{if $y_i$>0} \\ d_i = -1, & \text{if $y_i$≤0} \end{cases},θ_i = arctan(\frac{1}{2^i})\\

  2. 结果式:

    xn+1=Knn=0n(xn+dyn2n)yn+1=Knn=0n(yndxn2n)zn+1=z0+n=0ndθnwhere{d=+1,if yn>0d=1,if yn≤0, θn=arctan(12n)x_{n+1} = K_n\prod_{n=0}^{n}(x_n+d\frac{y_n}{2^n})\\ y_{n+1} = K_n\prod_{n=0}^{n}(y_n-d\frac{x_n}{2^n})\\ z_{n+1} = z_0 + \sum_{n=0}^{n}dθ_n \\ where \begin{cases} d = +1, & \text{if $y_n$>0} \\ d = -1, & \text{if $y_n$≤0} \end{cases},\ θ_n = arctan(\frac{1}{2^n})\\

  • 圆周坐标下旋转模式(rotation mode)的公式为
  1. 迭代式:

    {xi+1=xidiyi2iyi+1=yi+dixi2izi+1=zidiθiwhere{di=+1,if zi>0di=1,if zi≤0,θi=arctan(12i)\begin{cases} x_{i+1} = x_i-d_i\frac{y_i}{2^i}\\ y_{i+1} = y_i+d_i\frac{x_i}{2^i}\\ z_{i+1} = z_i - d_iθ_i\\ \end{cases} \\ where \begin{cases} d_i = +1, & \text{if $z_i$>0} \\ d_i = -1, & \text{if $z_i$≤0} \end{cases},θ_i = arctan(\frac{1}{2^i})\\

  2. 结果式:

    xn+1=Knn=0n(xndyn2n)yn+1=Knn=0n(yn+dxn2n)zn+1=z0n=0ndθnwhere{d=+1,if zn>0d=1,if zn≤0, θn=arctan(12n)x_{n+1} = K_n\prod_{n=0}^{n}(x_n-d\frac{y_n}{2^n})\\ y_{n+1} = K_n\prod_{n=0}^{n}(y_n+d\frac{x_n}{2^n})\\ z_{n+1} = z_0 - \sum_{n=0}^{n}dθ_n\\ where \begin{cases} d = +1, & \text{if $z_n$>0} \\ d = -1, & \text{if $z_n$≤0} \end{cases},\ θ_n = arctan(\frac{1}{2^n})\\

三种坐标系

上文只讨论了圆周坐标系,实际cordic算法可用于三种坐标系,进行不同的复杂运算。

  • 圆周坐标系(circular rotations)
  • 线性坐标系(linear rotations)
  • 双曲线坐标系(hyperbolic rotations)

cordic3

圆周坐标系

之前的所有讨论均基于圆周坐标系。这里不再赘述。

线性坐标系

假设(x1,y1)=(250,100)(x_1,y_1)=(250,100),要求其在线性坐标系下的旋转角θθ。首先根据观察可知tanθ=25tanθ=\frac{2}{5} 。这时利用cordic旋转算法进行尝试。

类似的,先逆时针"旋转"45°45°进行尝试。根据之前推导的公式可知,此时显然x2=x1=100x_2=x_1=100y2=y1x1tan45°=y1x1y_2=y_1-x_1tan45°=y_1-x_1。写成矩阵的形式:

[xy]=[10tan(45°)1][250100]=[1011][250100]=[250150]\left[ \begin{matrix} x\\ y \end{matrix} \right] = \left[ \begin{matrix} 1&0\\ -tan(45°)&1 \end{matrix} \right] \left[ \begin{matrix} 250\\ 100 \end{matrix} \right] = \left[ \begin{matrix} 1&0\\ -1&1 \end{matrix} \right] \left[ \begin{matrix} 250\\ 100 \end{matrix} \right] = \left[ \begin{matrix} 250\\ -150 \end{matrix} \right]

发现纵坐标y<0y<0,则旋转角度超过了范围,则在此基础上顺时针"回转"tan(θ)=12tan(θ)=\frac{1}{2}θθ角(θ=26.565051177078°θ=26.565051177078°)。写成矩阵的形式:

[xy]=[10tan(26.565051177078°)1][250150]=[10121][250150]=[25025]\left[ \begin{matrix} x\\ y \end{matrix} \right] = \left[ \begin{matrix} 1&0\\ tan(26.565051177078°)&1 \end{matrix} \right] \left[ \begin{matrix} 250\\ -150 \end{matrix} \right] = \left[ \begin{matrix} 1&0\\ \frac{1}{2}&1 \end{matrix} \right] \left[ \begin{matrix} 250\\ -150 \end{matrix} \right] = \left[ \begin{matrix} 250\\ -25 \end{matrix} \right]

发现纵坐标y<0y<0,则旋转角度依然不够,则在此基础上顺时针"回转"tan(θ)=14tan(θ)=\frac{1}{4}θθ角(θ=14.036243467926477°θ=14.036243467926477°)。写成矩阵的形式:

[xy]=[10tan(14.036243467926477°)1][25025]=[10141][25025]=[25037.5]\left[ \begin{matrix} x\\ y \end{matrix} \right] = \left[ \begin{matrix} 1&0\\ tan(14.036243467926477°)&1 \end{matrix} \right] \left[ \begin{matrix} 250\\ -25 \end{matrix} \right] = \left[ \begin{matrix} 1&0\\ \frac{1}{4}&1 \end{matrix} \right] \left[ \begin{matrix} 250\\ -25 \end{matrix} \right] = \left[ \begin{matrix} 250\\ 37.5 \end{matrix} \right]

发现纵坐标y>0y>0,则旋转角度超过了范围,则在此基础上继续逆时针"旋转"tan(θ)=18tan(θ)=\frac{1}{8}θθ角(θ=7.125016348901798°θ=7.125016348901798°)。写成矩阵的形式:

[xy]=[10tan(7.125016348901798°)1][25037.5]=[10181][25037.5]=[2506.25]\left[ \begin{matrix} x\\ y \end{matrix} \right] = \left[ \begin{matrix} 1&0\\ -tan(7.125016348901798°)&1 \end{matrix} \right] \left[ \begin{matrix} 250\\ 37.5 \end{matrix} \right] = \left[ \begin{matrix} 1&0\\ -\frac{1}{8}&1 \end{matrix} \right] \left[ \begin{matrix} 250\\ 37.5 \end{matrix} \right] = \left[ \begin{matrix} 250\\ 6.25 \end{matrix} \right]

这时纵坐标已经比较接近00了。类似的,按照这个模式继续算下去,可以不断逼近真实的结果。比如现在的结果范围为(x,y)=(250,6.25±15.625)(x,y)=(250,6.25±15.625),已经接近实际结果(x,y)=(250,0)(x,y)=(250,0)了。旋转过程演示图如下。

linear-rotate

线性坐标系的"旋转"其实与圆周坐标系下的旋转含义并不相同了,因为这里的矩阵变化(所谓旋转)其实表现的就是一种迭代映射,其变换逻辑图如下:

linear-map

注意,上述讨论分析的是使结果纵坐标y0y\to0cordic旋转变换,因此是向量模式(vector mode)

  • 线性坐标的向量模式(vector mode)的公式为(根据上述内容推导,zztan(y1/x1)tan(y1/x1),注意cordic算法规定迭代n从1开始,且由于没有迭代乘法,并不存在伸缩因子KnK_n
  1. 迭代式:

    {xi+1=xiyi+1=yidixi2izi+1=zi+dixi2iwhere{di=+1,if yi>0di=1,if yi≤0\begin{cases} x_{i+1} = x_i\\ y_{i+1} = y_i - d_i\frac{x_i}{2^i}\\ z_{i+1} = z_i + d_i\frac{x_i}{2^i}\\ \end{cases} \\ where \begin{cases} d_i = +1, & \text{if $y_i$>0} \\ d_i = -1, & \text{if $y_i$≤0} \end{cases}

  2. 结果式:

    xn+1=xn=x0yn+1=y0n=0ndxn2nzn+1=z0+n=1nd2nwhere{d=+1,if yn>0d=1,if yn≤0x_{n+1} = x_n = x_0\\ y_{n+1} = y_0 - \sum_{n=0}^{n}\frac{dx_n}{2^n}\\ z_{n+1} = z_0 + \sum_{n=1}^{n}\frac{d}{2^n} \\ where \begin{cases} d = +1, & \text{if $y_n$>0} \\ d = -1, & \text{if $y_n$≤0} \end{cases}

  • 圆周坐标下旋转模式(rotation mode)的公式为
  1. 迭代式:

    {xi+1=xiyi+1=yi+dxi2izi+1=zidxi2iwhere{d=+1,if zi>0d=1,if zi≤0\begin{cases} x_{i+1} = x_i\\ y_{i+1} = y_i + d\frac{x_i}{2^i}\\ z_{i+1} = z_i - d\frac{x_i}{2^i}\\ \end{cases} \\ where \begin{cases} d = +1, & \text{if $z_i$>0} \\ d = -1, & \text{if $z_i$≤0} \end{cases}

  2. 结果式:

    xn+1=xn=x0yn+1=y0+n=0ndxn2nzn+1=z0n=1nd2nwhere{d=+1,if zn>0d=1,if zn≤0x_{n+1} = x_n = x_0\\ y_{n+1} = y_0 + \sum_{n=0}^{n}\frac{dx_n}{2^n}\\ z_{n+1} = z_0 - \sum_{n=1}^{n}\frac{d}{2^n}\\ where \begin{cases} d = +1, & \text{if $z_n$>0} \\ d = -1, & \text{if $z_n$≤0} \end{cases}

双曲坐标系

双曲坐标系的cordic旋转比较类似于圆周坐标系。

应用问题

cordic功能表

图源为参考文献[6] ,侵删。

cordic_function_list

cordic计算范围及转换

cordic算法的计算是有范围限制的,超范围的值需要进行转换后才能利用cordic算法进行计算。

cordic计算误差

近似误差和舍入误差

见参考文献[7]

FPGA实现

总结


未完待续…长期施工中…

参考资料

[1] 邹熙. 基于CORDIC的指数函数的FPGA实现[J]. 大众科技, 2008, 000(010):36-37.

[2] 三角函数计算,Cordic 算法入门

[3] CORDIC_百度百科

[4] 周晓青, 李合生, 陶荣辉, et al. 基于CORDIC算法的双曲正余弦函数FPGA实现[J]. 太赫兹科学与电子信息学报, 2010(2):211-214.

[5] cordic算法详解 转载

[6] Xilinx CORDIC算法

[7] Y. H. Hu, “The Quantization Effects of the CORDIC Algorithm”, in IEEE Trans. On Signal Processing, Vol 40,
No 4, April 1992