首页 | 新闻 | 新品 | 文库 | 方案 | 视频 | 下载 | 商城 | 开发板 | 数据中心 | 座谈新版 | 培训 | 工具 | 博客 | 论坛 | 百科 | GEC | 活动 | 主题月 | 电子展
返回列表 回复 发帖

三角函数计算,Cordic算法入门

三角函数计算,Cordic算法入门

三角函数的计算是个复杂的主题,有计算机之前,人们通常通过查找三角函数表来计算任意角度的三角函数的值。这种表格在人们刚刚产生三角函数的概念的时候就已经有了,它们通常是通过从已知值(比如sin(π/2)=1)开始并重复应用半角和和差公式而生成。
现在有了计算机,三角函数表便推出了历史的舞台。但是像我这样的喜欢刨根问底的人,不禁要问计算机又是如何计算三角函数值的呢。最容易想到的办法就是利用级数展开,比如泰勒级数来逼近三角函数,只要项数取得足够多就能以任意的精度来逼近函数值。除了泰勒级数逼近之外,还有其他许多的逼近方法,比如切比雪夫逼近、最佳一致逼近和Padé逼近等。
所有这些逼近方法本质上都是用多项式函数来近似我们要计算的三角函数,计算过程中必然要涉及到大量的浮点运算。在缺乏硬件乘法器的简单设备上(比如没有浮点运算单元的单片机),用这些方法来计算三角函数会非常的费时。为了解决这个问题,J. Volder于1959年提出了一种快速算法,称之为CORDIC(COordinate Rotation DIgital Computer) 算法,这个算法只利用移位和加减运算,就能计算常用三角函数值,如Sin,Cos,Sinh,Cosh等函数。 J. Walther在1974年在这种算法的基础上进一步改进,使其可以计算出多种超越函数,更大的扩展了Cordic 算法的应用。因为Cordic 算法只用了移位和加法,很容易用纯硬件来实现,因此我们常能在FPGA运算平台上见到它的身影。不过,大多数的软件程序员们都没有听说过这种算法,也更不会主动的去用这种算法。其实,在嵌入式软件开发,尤其是在没有浮点运算指令的嵌入式平台(比如定点型DSP)上做开发时,还是会遇上可以用到Cordic 算法的情况的,所以掌握基本的Cordic算法还是有用的。
从二分查找法说起
先从一个例子说起,知道平面上一点在直角坐标系下的坐标(X,Y)=(100,200),如何求的在极坐标系下的坐标(ρ,θ)。用计算器计算一下可知答案是(223.61,63.435)。

图 1 直角坐标系到极坐标系的转换

为了突出重点,这里我们只讨论X和Y都为正数的情况。这时θ=atan(y/x)。求θ的过程也就是求atan 函数的过程。Cordic算法采用的想法很直接,将(X,Y)旋转一定的度数,如果旋转完纵坐标变为了0,那么旋转的度数就是θ。坐标旋转的公式可能大家都忘了,这里把公式列出了。设(x,y)是原始坐标点,将其以原点为中心,顺时针旋转θ之后的坐标记为(x’,y’),则有如下公式:

也可以写为矩阵形式:

如何旋转呢,可以借鉴二分查找法的思想。我们知道θ的范围是0到90度。那么就先旋转45度试试。

旋转之后纵坐标为70.71,还是大于0,说明旋转的度数不够,接着再旋转22.5度(45度的一半)。

这时总共旋转了45+22.5=67.5度。结果纵坐标变为了负数,说明θ<67.5度,这时就要往回转,还是二分查找法的思想,这次转11.25度。

这时总共旋转了45+22.5-11.25=56.25度。又转过头了,接着旋转,这次顺时针转5.625度。

这时总共旋转了45+22.5-11.25+5.625=61.875度。这时纵坐标已经很接近0了。我们只是说明算法的思想,因此就不接着往下计算了。计算到这里我们给的答案是 61.875±5.625。二分查找法本质上查找的是一个区间,因此我们给出的是θ值的一个范围。同时,坐标到原点的距离ρ也求出来了,ρ=223.52。与标准答案比较一下计算的结果还是可以的。旋转的过程图示如下。

可能有读者会问,计算中用到了 sin 函数和 cos 函数,这些值又是怎么计算呢。很简单,我们只用到很少的几个特殊点的sin 函数和 cos 函数的值,提前计算好存起来,用时查表。
下面给出上面方法的C 语言实现。
#include
#include

double my_atan2(double x, double y);
int main(void)
{
double z = my_atan2(100.0, 200.0);
printf("\n z = %f \n", z);

return 0;
}

double my_atan2(double x, double y)
{
const double sine[] = {0.7071067811865,0.3826834323651,0.1950903220161,0.09801714032956,
0.04906767432742,0.02454122852291,0.01227153828572,0.006135884649154,0.003067956762966
,0.001533980186285,7.669903187427045e-4,3.834951875713956e-4,1.917475973107033e-4,
9.587379909597735e-5,4.793689960306688e-5,2.396844980841822e-5
};

const double cosine[] = {0.7071067811865,0.9238795325113,0.9807852804032,0.9951847266722,
0.9987954562052,0.9996988186962,0.9999247018391,0.9999811752826,0.9999952938096,
0.9999988234517,0.9999997058629,0.9999999264657,0.9999999816164,0.9999999954041,
0.999999998851,0.9999999997128
};

int i = 0;
double x_new, y_new;
double angleSum = 0.0;
double angle = 45.0;

for(i = 0; i < 15; i++)
{
if(y > 0)
{
x_new = x * cosine + y * sine;
y_new = y * cosine - x * sine;
x = x_new;
y = y_new;
angleSum += angle;
}
else
{
x_new = x * cosine - y * sine;
y_new = y * cosine + x * sine;
x = x_new;
y = y_new;
angleSum -= angle;
}
printf("Debug: i = %d angleSum = %f, angle = %f\n", i, angleSum, angle);
angle /= 2;
}
return angleSum;
}

程序运行的输出结果如下:
Debug: i = 0 angleSum = 45.000000, angle = 45.000000
Debug: i = 1 angleSum = 67.500000, angle = 22.500000
Debug: i = 2 angleSum = 56.250000, angle = 11.250000
Debug: i = 3 angleSum = 61.875000, angle = 5.625000
Debug: i = 4 angleSum = 64.687500, angle = 2.812500
Debug: i = 5 angleSum = 63.281250, angle = 1.406250
Debug: i = 6 angleSum = 63.984375, angle = 0.703125
Debug: i = 7 angleSum = 63.632813, angle = 0.351563
Debug: i = 8 angleSum = 63.457031, angle = 0.175781
Debug: i = 9 angleSum = 63.369141, angle = 0.087891
Debug: i = 10 angleSum = 63.413086, angle = 0.043945
Debug: i = 11 angleSum = 63.435059, angle = 0.021973
Debug: i = 12 angleSum = 63.424072, angle = 0.010986
Debug: i = 13 angleSum = 63.429565, angle = 0.005493
Debug: i = 14 angleSum = 63.432312, angle = 0.002747

z = 63.432312
减少乘法运算
现在已经有点 Cordic 算法的样子了,但是我们看到没次循环都要计算 4 次浮点数的乘法运算,运算量还是太大了。还需要进一步的改进。改进的切入点当然还是坐标变换的过程。我们将坐标变换公式变一下形。

可以看出 cos(θ)可以从矩阵运算中提出来。我们可以做的再彻底些,直接把 cos(θ) 给省略掉。

省略cos(θ)后发生了什么呢,每次旋转后的新坐标点到原点的距离都变长了,放缩的系数是1/cos(θ)。不过没有关系,我们求的是θ,不关心ρ的改变。这样的变形非常的简单,但是每次循环的运算量一下就从4次乘法降到了2次乘法了。
还是给出 C 语言的实现:
double my_atan3(double x, double y)
{
const double tangent[] = {1.0,0.4142135623731,0.1989123673797,0.09849140335716,0.04912684976947,
0.02454862210893,0.01227246237957,0.006136000157623,0.003067971201423,
0.001533981991089,7.669905443430926e-4,3.83495215771441e-4,1.917476008357089e-4,
9.587379953660303e-5,4.79368996581451e-5,2.3968449815303e-5
};

int i = 0;
double x_new, y_new;
double angleSum = 0.0;
double angle = 45.0;

for(i = 0; i < 15; i++)
{
if(y > 0)
{
x_new = x + y * tangent;
y_new = y - x * tangent;
x = x_new;
y = y_new;
angleSum += angle;
}
else
{
x_new = x - y * tangent;
y_new = y + x * tangent;
x = x_new;
y = y_new;
angleSum -= angle;
}
printf("Debug: i = %d angleSum = %f, angle = %f, ρ = %f\n", i, angleSum, angle, hypot(x,y));
angle /= 2;
}
return angleSum;
}

计算的结果是:
Debug: i = 0 angleSum = 45.000000, angle = 45.000000, ρ = 316.227766
Debug: i = 1 angleSum = 67.500000, angle = 22.500000, ρ = 342.282467
Debug: i = 2 angleSum = 56.250000, angle = 11.250000, ρ = 348.988177
Debug: i = 3 angleSum = 61.875000, angle = 5.625000, ρ = 350.676782
Debug: i = 4 angleSum = 64.687500, angle = 2.812500, ρ = 351.099697
Debug: i = 5 angleSum = 63.281250, angle = 1.406250, ρ = 351.205473
Debug: i = 6 angleSum = 63.984375, angle = 0.703125, ρ = 351.231921
Debug: i = 7 angleSum = 63.632813, angle = 0.351563, ρ = 351.238533
Debug: i = 8 angleSum = 63.457031, angle = 0.175781, ρ = 351.240186
Debug: i = 9 angleSum = 63.369141, angle = 0.087891, ρ = 351.240599
Debug: i = 10 angleSum = 63.413086, angle = 0.043945, ρ = 351.240702
Debug: i = 11 angleSum = 63.435059, angle = 0.021973, ρ = 351.240728
Debug: i = 12 angleSum = 63.424072, angle = 0.010986, ρ = 351.240734
Debug: i = 13 angleSum = 63.429565, angle = 0.005493, ρ = 351.240736
Debug: i = 14 angleSum = 63.432312, angle = 0.002747, ρ = 351.240736

z = 63.432312
返回列表