定点数三角函数:CORDIC 算法(3)

作于:2025-08-25 14:30 ,预计阅读时间 4 分钟

前言

本文探讨学习 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
}

/定点数/ /算法/ /数学/ /数值计算/