前言
本文探讨学习 CORDIC 算法的向量模式与向量模式下的 $\arctan$ 实现。
关于 CORDIC 算法的原理和旋转模式可以参考前两篇文章,下文不再重复推导:
推导过程
对于平面向量 $\vec{V}=(x,y)$ 和 $x$ 轴的夹角 $\theta$ 有 $\tan(\theta) = y/x$
反过来就是 $\arctan(y/x) = \theta$
我们可以把 $\arctan(y)$ 看作 $\arctan(y/1)$,构造出向量 $V=(1,y)$
然后不断以微小角度旋转向量,直到 $V$ 接近对齐 $x$ 轴正方向,也就是 $y \approx 0$,此时累积的旋转角度 $z$ 就是向量 $V$ 和 $x$ 轴的夹角 $\theta$ 的近似了。
旋转过程
我们可以复用前文推导出的 CORDIC 旋转模式的迭代公式(无缩放):
$$ \begin{cases} x_{i+1} = x_i - y_i \cdot d_i \cdot 2^{-i} \\ y_{i+1} = y_i + x_i \cdot d_i \cdot 2^{-i} \end{cases} $$
其中,$d_i$ 的选择取决于 $y_i$。当 $y_i>0$ 时,向量在第一象限,旋转方向顺时针。
反之,$y_i<0$ 则说明向量在第四象限,旋转方向为逆时针。
$$ d_i = \begin{cases} -1 &\text{if } y_i \gt 0 \\ 1 &\text{if } y_i \lt 0 \end{cases} $$
角度累加器 $z$ 和旋转模式一样,采用查表法和 $d_i$ 确定符号。
值得注意的是,缩放因子 K 并未出现在上述公式里,因为我们不关心 x,y 的取值,只关心累加器 z 的结果。在几何意义上来说,K 其实是影响了向量的长度(因为乘上K实际是应用了一个向量的标量乘法),而不影响向量的角度。
边界情况
前文旋转模式实现中提到,CORDIC 算法旋转的极限大约是 99°。
由于我们构造的向量中,$x$ 分量固定为 1,所以向量一定会落在第一象限或第四象限内。
$\arctan$ 的值域也是 $[-\frac{pi}{2},\frac{pi}{2}]$,在 CORDIC 算法可以处理的范围内。
$\arctan$ 的定义域是实数 $\arctan(x), x \in \mathbb{R}$,实际实现中则受制于精度,可以近似地看作所有 Q16 类型取值都是合法的输入。
实现
实现中一个值得关注的点在于 $z$ 的符号。我们的旋转累积是向 $x$ 轴方向的,和旋转模式正好相反。所以 $z$ 累积的值应该做一个取反操作。
举例来说,旋转模式从0度转到45度,那么向量模式就是从45度转到零度。他们记录的旋转累积一个是45一个是-45,正好相反。
func ATan(v Q16) (Q16, error) {
var (
x Q16 = MustFromInt(1)
y Q16 = v
d Q16
z Q16
err error
)
for i := range atanLookupTable {
var (
nx Q16
ny Q16
)
if y > 0 {
d = -1
} else if y < 0 {
d = 1
} else {
break
}
intermediate := y >> i
nx, err = Sub(x, d*intermediate)
if err != nil {
return 0, err
}
intermediate = x >> i
ny, err = Add(y, d*intermediate)
if err != nil {
return 0, err
}
x = nx
y = ny
z += d * atanLookupTable[i]
}
return -z, nil
}