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

三角函数计算,Cordic算法入门(2)

三角函数计算,Cordic算法入门(2)

消除乘法运算
我们已经成功的将乘法的次数减少了一半,还有没有可能进一步降低运算量呢?还要从计算式入手。

第一次循环时,tan(45)=1,所以第一次循环实际上是不需要乘法运算的。第二次运算呢?
Tan(22.5)=0.4142135623731,很不幸,第二次循环乘数是个很不整的小数。是否能对其改造一下呢?答案是肯定的。第二次选择22.5度是因为二分查找法的查找效率最高。如果选用个在22.5到45度之间的值,查找的效率会降低一些。如果稍微降低一点查找的效率能让我们有效的减少乘法的次数,使最终的计算速度提高了,那么这种改进就是值得的。

我们发现tan(26.565051177078)=0.5,如果我们第二次旋转采用26.565051177078度,那么乘数变为0.5,如果我们采用定点数运算的话(没有浮点协处理器时为了加速计算我们会大量的采用定点数算法)乘以0.5就相当于将乘数右移一位。右移运算是很快的,这样第二次循环中的乘法运算也被消除了。
类似的方法,第三次循环中不用11.25度,而采用 14.0362434679265 度。
Tan(14.0362434679265)= 1/4

乘数右移两位就可以了。剩下的都以此类推。
tan(45)= 1
tan(26.565051177078)= 1/2
tan(14.0362434679265)= 1/4
tan(7.1250163489018)= 1/8
tan(3.57633437499735)= 1/16
tan(1.78991060824607)= 1/32
tan(0.8951737102111)= 1/64
tan(0.4476141708606)= 1/128
tan(0.2238105003685)= 1/256

还是给出C语言的实现代码,我们采用循序渐进的方法,先给出浮点数的实现(因为用到了浮点数,所以并没有减少乘法运算量,查找的效率也比二分查找法要低,理论上说这个算法实现很低效。不过这个代码的目的在于给出算法实现的示意性说明,还是有意义的)。
double my_atan4(double x, double y)
{
const double tangent[] = {1.0, 1 / 2.0, 1 / 4.0, 1 / 8.0, 1 / 16.0,
1 / 32.0, 1 / 64.0, 1 / 128.0, 1 / 256.0, 1 / 512.0,
1 / 1024.0, 1 / 2048.0, 1 / 4096.0, 1 / 8192.0, 1 / 16384.0
};
const double angle[] = {45.0, 26.565051177078, 14.0362434679265, 7.1250163489018, 3.57633437499735,
1.78991060824607, 0.8951737102111, 0.4476141708606, 0.2238105003685, 0.1119056770662,
0.0559528918938, 0.027976452617, 0.01398822714227, 0.006994113675353, 0.003497056850704
};

int i = 0;
double x_new, y_new;
double angleSum = 0.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));
}
return angleSum;
}

程序运行的输出结果如下:
Debug: i = 0 angleSum = 45.000000, angle = 45.000000, ρ = 316.227766
Debug: i = 1 angleSum = 71.565051, angle = 26.565051, ρ = 353.553391
Debug: i = 2 angleSum = 57.528808, angle = 14.036243, ρ = 364.434493
Debug: i = 3 angleSum = 64.653824, angle = 7.125016, ρ = 367.270602
Debug: i = 4 angleSum = 61.077490, angle = 3.576334, ρ = 367.987229
Debug: i = 5 angleSum = 62.867400, angle = 1.789911, ρ = 368.166866
Debug: i = 6 angleSum = 63.762574, angle = 0.895174, ρ = 368.211805
Debug: i = 7 angleSum = 63.314960, angle = 0.447614, ρ = 368.223042
Debug: i = 8 angleSum = 63.538770, angle = 0.223811, ρ = 368.225852
Debug: i = 9 angleSum = 63.426865, angle = 0.111906, ρ = 368.226554
Debug: i = 10 angleSum = 63.482818, angle = 0.055953, ρ = 368.226729
Debug: i = 11 angleSum = 63.454841, angle = 0.027976, ρ = 368.226773
Debug: i = 12 angleSum = 63.440853, angle = 0.013988, ρ = 368.226784
Debug: i = 13 angleSum = 63.433859, angle = 0.006994, ρ = 368.226787
Debug: i = 14 angleSum = 63.437356, angle = 0.003497, ρ = 368.226788

z = 63.437356
有了上面的准备,我们可以来讨论定点数算法了。所谓定点数运算,其实就是整数运算。我们用256 表示1度。这样的话我们就可以精确到1/256=0.00390625 度了,这对于大多数的情况都是足够精确的了。256 表示1度,那么45度就是 45*256 = 115200。其他的度数以此类推。
序号
度数
×256
取整
1
45.0
11520
11520
2
26.565051177078
6800.65310133196
6801
3
14.0362434679265
3593.27832778918
3593
4
7.1250163489018
1824.00418531886
1824
5
3.57633437499735
915.541599999322
916
6
1.78991060824607
458.217115710994
458
7
0.8951737102111
229.164469814035
229
8
0.4476141708606
114.589227740302
115
9
0.2238105003685
57.2954880943458
57
10
0.1119056770662
28.647853328949
29
11
0.0559528918938
14.3239403248137
14
12
0.027976452617
7.16197186995294
7
13
0.01398822714227
3.58098614841984
4
14
0.006994113675353
1.79049310089035
2
15
0.003497056850704
0.8952465537802
1

C 代码如下:
int my_atan5(int x, int y)
{
const int angle[] = {11520, 6801, 3593, 1824, 916, 458, 229, 115, 57, 29, 14, 7, 4, 2, 1};

int i = 0;
int x_new, y_new;
int angleSum = 0;

x *= 1024;// 将 X Y 放大一些,结果会更准确
y *= 1024;

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

计算结果如下:
Debug: i = 0 angleSum = 11520, angle = 11520
Debug: i = 1 angleSum = 18321, angle = 6801
Debug: i = 2 angleSum = 14728, angle = 3593
Debug: i = 3 angleSum = 16552, angle = 1824
Debug: i = 4 angleSum = 15636, angle = 916
Debug: i = 5 angleSum = 16094, angle = 458
Debug: i = 6 angleSum = 16323, angle = 229
Debug: i = 7 angleSum = 16208, angle = 115
Debug: i = 8 angleSum = 16265, angle = 57
Debug: i = 9 angleSum = 16236, angle = 29
Debug: i = 10 angleSum = 16250, angle = 14
Debug: i = 11 angleSum = 16243, angle = 7
Debug: i = 12 angleSum = 16239, angle = 4
Debug: i = 13 angleSum = 16237, angle = 2
Debug: i = 14 angleSum = 16238, angle = 1

z = 16238
16238/256=63.4296875度,精确的结果是63.4349499度,两个结果的差为0.00526,还是很精确的。
到这里 CORDIC 算法的最核心的思想就介绍完了。当然,这里介绍的只是CORDIC算法最基本的内容,实际上,利用CORDIC 算法不光可以计算 atan 函数,其他的像 Sin,Cos,Sinh,Cosh 等一系列的函数都可以计算,不过那些都不在本文的讨论范围内了。另外,每次旋转时到原点的距离都会发生变化,而这个变化是确定的,因此可以在循环运算结束后以此补偿回来,这样的话我们就同时将(ρ,θ)都计算出来了。

想进一步深入学习的可以阅读 John Stephen Walther 于2000年发表在 Journal of VLSI signal processing systems for signal, image and video technology上的综述性文章“The Story of Unified Cordic”。
返回列表