CORIDIC算法学习记录

目录

• 问题
• 问题分析
• CORDIC算法原理
  • 逼近方法及步骤
  • 逼近过程中的符号确定
  • 根据角度计算正切值
  • 举个例子逼近\(\theta=50^{\degree}\)并求其正切值

CORDIC算法叫坐标旋转数字计算法，由J.Volder在1959年提出，可以快速且简单的计算角度的数值。

问题

已知\(y,x\)，如何快速计算角度\(\theta\)？

问题分析

许多人，直接求解\(\arctan(\frac{y}{x})\)，但是反正切求解没有那么容易，如果要查表，由不得又要一堆操作下来，费事费力。虽然可以借助计算机求解，计算机只能加减乘除，计算反正切，需要级数展开，根据泰勒公式\(\arctan(x)=x-\frac{x^3}{3}+\frac{x^5}{5}-\frac{x^7}{7}+\ldots\)，可见包含了大量的乘法和除法，也不会快到哪里。

CORDIC算法原理

逼近方法及步骤

• 选择一系列标准角度，A=[45 26.5651 14.0362 7.12502 3.57633 1.78991 0.895174 0.447614 0.223811 0.111906 0.0559529 0.0279765 0.0139882 0.0069941]
• 将任何要计算的角度\(\theta\)与比较对象\(\beta\)比较，
• 第一次\(\theta\)与\(\beta\)取\(A_{(1)}=45^{\degree}\)比较，
• 若\(\theta>\beta=45^{\degree}\)时，增大\(\beta\)，将\(\beta\)加上\(26.5651^{\degree}\)
• 再次比较\(\theta\)和\(\beta\)大小，当\(\theta > \beta\)就继续增加\(A\)的下一项，当\(\theta < \beta\)就继续减去\(A\)的下一项
• 直到\(\theta\)和\(\beta\)足够接近即可停止

例如：待测角度为\(53^{\degree}\)，逼近过程用表格表示如下：

序号	初始值\(\beta\)	比较结果	变化的值	符号	备注
1	45	1	26.5651	+
2	71.5651	0	14.0362	-
3	57.5289	0	7.12502	-
4	50.4039	1	3.57633	+
5	53.9802	0	1.78991	-
6	52.1903	1	0.895174	+
7	53.0855	0	0.447614	-
8	52.611186	1	0.223811	+
9	52.834997	1	0.111906	+
10	52.946903	1	0.0559529	+
11	53.0028559	0	0.0279765	-
12	52.9748794	1	0.0139882	+
13	52.9888676	1	0.0069941	+

最终为52.9958617，可以看到距离53已经很接近了，如果不够还可以继续。

可以看到，上边的比较过程，用眼观察并判断下一项的正负号比较容易，有人觉得直接比较大小，随着小数位数的增加，比较大小并不可靠。

> if(0.999999999999999999999999<1)
> a=10;
> end

以上代码中，我这里octave9.4中，小数点后16以前比较大小正确，超过就可能出错了，编程中小数比较大小不一定可靠。

逼近过程中的符号确定

• 将原角度对应的向量\((x,y)\)做相应角度的旋转，根据旋转后的纵坐标是否大于零，就知道下一次操作的正负号了
• 若旋转后的坐标\((x_1,y_1)\)，根据欧拉公式\(x_1=x\cos\alpha-y\sin\alpha\),\(y_1=x\sin\alpha+y\cos\alpha\)，其中\(\alpha\)为旋转的角度，方向与正负号有关
• 为了避免计算正余弦值，改造公式为：\(x_1=\cos\alpha(1-y\tan\alpha\)，\(y_1=\cos\alpha(x\tan\alpha+y)\)，我们只要能判断旋转后纵坐标\(y_1\)是否大于零，就能确定下次的旋转方向，无需计算\(\cos\alpha\)，计算\(x\tan\alpha\)和\(y\tan\alpha\)
• 前面我们选择角度时取值的特征就是\(\tan(45^{\degree})=1\)，依次为\(\frac{1}{2}，\frac{1}{4}，\frac{1}{8}，\frac{1}{16}\ldots\)，也就是说按照表里的角度值逼近时，无需计算正切值，每次都是前一个的\(\frac{1}{2}\)。是不是需要除法，计算机计算不需要，直接移位一次就能乘以或除以2了，是不是就很方便了。

根据角度计算正切值

• 对角度\(\theta\),则\(\tan\theta=\frac{\sin\theta}{\cos\theta}\)，计算出\(\sin\theta,\cos\theta\)；
• 由于\(\sin\theta=\frac{y}{r}\)，\(\cos\theta=\frac{x}{r}\)，单位圆内，知道\(x,y\)即可；

举个例子逼近\(\theta=50^{\degree}\)并求其正切值

使用的欧拉公式：

\[x_1=\cos\alpha(x-y\tan\alpha)\\ y_1=\cos\alpha(x\tan\alpha+y) \]

0. 向量\((1,0)\)开始，向量的角度\(\alpha\)为\(0^{\degree}\)
1. 小于目标，需要增大\(\alpha\)，增大\(45^{\degree}\),$\begin{bmatrix}x_1\y_1\end{bmatrix}=\cos{(45^{\degree})}\begin{bmatrix}x-y\tan 45^{\degree}\x\tan 45^{{\degree}+y\end{bmatrix}=\cos45}\begin{bmatrix}x_1^{'}\y_1\end{bmatrix} \(，这里只需要移位和加减即可得到\)\begin{bmatrix}x_1^{'}\y_1\end{bmatrix}$,前面的余弦值计算稍后处理
2. 比较新的\(\alpha\)，依旧小于目标，需要继续增加，虽然离目标只有\(5^{\degree}\),无法直接计算正切还有坐标，继续选择标准角度26.5651，经过旋转后，新坐标：$ \begin{bmatrix}x_1\y_1\end{bmatrix}=\cos{(45^{{\degree})}\cos{(26.5651})}\begin{bmatrix}x_1^{'}-y_1\tan 26.5651^{{\degree}\x_1}\tan 26.5651^{{\degree}+y_1}\end{bmatrix}=\cos45^{{\degree}\cos{(26.5651})}\begin{bmatrix}x_2^{'}\y_2\end{bmatrix} $此时的角度：\(71.5651\)，
3. 比较新的\(\alpha\)，大于目标，需要减小，选择角度-14.0362，就这样一直进行下去，最终角度会越来越接近50度，例如经过8次旋转，角度为：\(45+26.5651-14.0362-7.12502-3.56733+1.78911+0.895174+0.447614=49.9601\)，此时的误差已经很小了，如果不够继续即可；
4. \(\begin{bmatrix}x_8\\ y_8 \end{bmatrix}=\cos(45)^{\degree}\cos(26.5651)^{\degree}\cos(-14.0362)^{\degree}\cos(-7.12502)^{\degree}\cos(-3.56733)^{\degree}\cos(1.78911)^{\degree}\cos(0.895174)^{\degree}\cos(0.447614)^{\degree}\begin{bmatrix}x_8^{'}\\ y_8^{'} \end{bmatrix}\)

而且以上连乘序列最终收敛到0.6072529，由于\(\cos\alpha=\frac{1}{\sqrt{1+(\tan\alpha)^2}}\),也就是说以上序列就成了$\prod^{n_{i=0}\frac{1}{\sqrt{1+{(2})}^2}} $

验证如下：

N=20;
A=zeros(1,N);
for i=0:1:N-1A(i+1)=1./sqrt(1+2^(-2*i));
end
B=cumprod(A);i=1:N;
%双y轴做图
yyaxis left
%scatter(i,A,'markeredgecolor',[0 .7 .8],'markerfacecolor',[0 .5 .6],'linewidth',2.5);
scatter(i,A,'markeredgecolor','r','markerfacecolor','g','linewidth',2.5);
xlim([-1 N+1]);
ylim([0.55 1.05]);
title('difference A and B');
xlabel('i');
ylabel('A changes')yyaxis right
%scatter(i,B);
scatter(i,B,'markeredgecolor','b','markerfacecolor','r','linewidth',1.5);
ylim([0.55 1.05]);
ylabel('B,changes')

结果图就不放了，也就是说该多项式大概从第六项开始相乘时，就趋近与固定值0.6073了，以后无论怎么加项，角度在逼近，精度在提高，系统不再发生变化，随着N增加，A从第四项开始进入1%误差范围，从第七项开始，N每增加1，精度提高一个数量级，直到第27项开始matlab精度达到极值时，值为1.000000000000000不再变化。与A对应的B，从第4项开始进入到0.6系数中，从第7项开始进入0.6275开始进入精度稳定，直到第26项值为0.607252935008881，保持稳定不再变化。也就是说，才用CORDIC算法经过26次迭代，就可以实现几乎所有的精度要求。

参考：https://zhuanlan.zhihu.com/p/717126119