引言

目前在计算机中，三角函数就大多是采用级数展开的方法来计算的，例如使用泰勒展开等逼近。只要展开的项数足够多，就能达到所需要的计算精度。标准C语言的math库的三角函数，就大多是使用泰勒展开的方法实现的。但是使用级数展开的方法，基本需要依赖于浮点运算。尽管现代的很多计算机设备都有浮点计算功能，但还有很大一些嵌入式设备是不支持浮点运算的，使用级数展开的方法会造成较大的计算负担。

在嵌入式系统和数字信号处理（DSP）应用中，三角函数计算是非常常见的需求。然而，传统的浮点运算在资源受限的嵌入式平台上往往效率低下。CORDIC（Coordinate Rotation Digital Computer）算法提供了一种高效的解决方案，它仅使用移位和加减法就能实现三角函数、双曲函数等多种数学运算。

本文的主要内容是介绍CORDIC算法，实现简单的定点的三角函数计算，并提供简单的arctan()、cos()、sin()函数的实现代码。使用的编程语言是C语言。

CORDIC 算法原理

CORDIC（Coordinate Rotation Digital Computer）算法即坐标旋转数字计算方法，是J.D.Volder1于1959年首次提出的。

CORDIC算法本质上是一个二分查找法，利用已知的一些数据进行多次的比较，不断缩小所求的目标的区间，计算出达到所要求精度的目标值。

CORDIC算法如其名，是通过旋转的方法来求出目标函数的值。为了方便解释，以下说明的要计算的角度都是处于第一象限的。实际上其他象限角度的三角函数可以通过三角变换转换到第一象限进行计算。

咱们从计算反正切函数arctan() 入手。

设 （x,y）是原始坐标点，将其以原点为中心顺时针旋转 θ角，新的坐标记为 （x',y'）,则有如下公式：

$$
\left[ \begin{matrix} x’ \\ y’ \\ \end{matrix} \right] =\left[ \begin{matrix} cos(\theta) & sin(\theta) \\ -sin(\theta) & cos(\theta) \\ \end{matrix} \right] \left[ \begin{matrix} x \\ y \\ \end{matrix} \right]
$$

我们化简一下，把 cos()提取出来。实际上，在这里我们比较关心的是 x、y、tan()的值，cos()在矩阵外面就是个倍数，影响的主要是旋转角的半径，不是我们关注的量。因此可以干脆点，直接省略掉 cos()，这样矩阵就只剩下 tan()函数，变成如下：

$$
\left[ \begin{matrix} x’ \\ y’ \\ \end{matrix} \right] = \left[ \begin{matrix} 1 & tan(\theta) \\ -tan(\theta) & 1 \\ \end{matrix} \right] \left[ \begin{matrix} x \\ y \\ \end{matrix} \right]
$$

我们所求的 arctan()，实际上就是已知 tan(θ)的值，要求 θ。换句话讲，就是上面的矩阵中，x、y、tan()中就都是我们已知的值，我们要求出 θ。CORDIC算法就是通过多次顺时针旋转，每次旋转都经过上面的矩阵计算出新的坐标值 x'、y'，累积旋转过的角度。我们的目标就是要让这个角旋转到0度角，也就是纵坐标为0的时候，那么此时多次旋转累积的角度即为所求的目标角度。 当旋转之后的纵坐标小于零，说明我们旋转过头了，下一次旋转就要进行逆时针旋转；当旋转之后的纵坐标大于零，说明还没转到零度角，下一次还是继续进行逆时针旋转。

那么新的问题出现了，我们希望能尽量减少计算量，每次旋转的角度确定为多少比较合适呢？这里就结合了二分查找的思想和定点型计算的特点。

在计算机的各种运算中，移位运算相比于乘除法运算来说速度快，因此我们这里使用左移代替乘法（乘以2^N），右移代替除法（除以2^N）。所以在每次旋转时，我们就要保证变化的量（此处是指 x、y）都是乘以或者除以2的N次方。在上面的矩阵中，x、y每次旋转都是乘以 tan(θ)，所以我们要保证每次旋转的 tan(θ)的值是1、1/2、1/4、1/8……那么我们每次旋转的角度就是当 tan(θ)=1、1/2、1/4、1/8……的 θ角度，也就是45°、26.5°、14.03°、7.125°……这里就是类似于二分查找的过程，而纵坐标为0就是要查找到的目标值。

CORDIC算法的实现思路

CORDIC算法的大致过程现在已经了解了，接下来就得来完善一下代码实现的一些细节问题。

1.确定放大倍数

由于没有浮点运算和浮点数，我们需要使用放大的方式来表示上面谈到的各种角度和数值。放大的倍数可以根据自己所需要的计算精度来确定。另外，两个放大后的数进行乘除法的时候，会被放大两次或者放大倍数被抵消，因此需要特殊处理一下。

在这里，设定数值放大的倍数为256，也就是Q8。角度同样也是放大256的，也就是浮点数1=定点256，角度值1°=定点256。

2.确定每次旋转的角度及其tan值

上面已经说到我们每次旋转时的角度都是45°、26.5°、14.03°、7.125°……对应角度以及tan值如下表。这些数据我们就当作一些常量数据放起来，计算三角函数时就利用这些已知的值进行计算。

tan(θ) θ Q8放大之后
1 45.0° 11520
1/2 26.565051177078° 6801
1/4 14.0362434679265° 3593
1/8 7.1250163489018° 1824
1/16 3.57633437499735° 916
1/32 1.78991060824607° 458
1/64 0.8951737102111° 229
1/128 0.4476141708606° 115
1/256 0.2238105003685° 57
1/512 0.1119056770662° 29
1/1024 0.0559528918938° 14
1/2048 0.027976452617° 7
1/4096 0.01398822714227° 4
1/8192 0.006994113675353° 2
1/16384 0.003497056850704° 1

因为到了 tan(θ)=1/16384的时候，放大后的数据已经是1了，定点数无法再缩小了，因此放大256倍的精度只能到这了。

CORDIC算法的代码实现

cordic.h

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#ifndef _CORDIC_H
#define _CORDIC_H


/**********************************************************************************
1. The file provides arctan, cos and sin functions using CORDIR algorithm.
2. In the file, 1 degree = 256, float 1 = 256.
***********************************************************************************/

typedef unsigned char u8, uint8_t;
typedef signed char s8, int8_t;
typedef unsigned short u16, uint16_t;
typedef signed short s16, int16_t;
typedef unsigned long u32, uint32_t;
typedef signed long s32, int32_t;


s32 my_atan(s32 x, s32 y);
s32 my_cos(s32 angle);
s32 my_sin(s32 angle);

#endif // _CORDIC_H

cordic.c

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
#include "cordic.h"


/**********************************************************************************
1. The file provides arctan, cos and sin functions using CORDIR algorithm.
2. In the file, 1 degree = 256, float 1 = 256.
***********************************************************************************/


const u32 tan_angle[] = {11520, 6801, 3593, 1824, 916, 458, 229, 115, 57, 29, 14, 7, 4, 2, 1};


/*
 * Brief   : the arctan function.
 * Param   :
     - x : the abscissa of angle.
     - y : the ordinate of angle.
 * Return  : the magnified result.
 */
s32 my_atan(s32 x, s32 y)
{
    u8 quad = 0; //quadrant
    if (x < 0) {
        quad = 2;
        x = -x; //x change to positive number.
        if (y < 0)
            quad = 3;
    }

    x = (x << 10);
    y = (y << 10);

    s32 x_new, y_new;
    s32 angleSum = 0;

    for (u8 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 += tan_angle[i];
        } else if (y < 0) {
            x_new = x - (y >> i);
            y_new = y + (x >> i);
            x = x_new;
            y = y_new;
            angleSum -= tan_angle[i];
        } else {
            break;
        }
    }

    if (quad == 2)
        angleSum = 46080 - angleSum;
    else if (quad == 3)
        angleSum = -46080 - angleSum;

    return angleSum;
}

/*
 * Brief   : the cos function.
 * Param   :
     - angle : input angle (magnified).
 * Return  : the magnified result.
 */
s32 my_cos(s32 angle)
{
    angle = angle % 92160;
    if (angle > 46080)
        angle -= 92160;
    else if (angle < -46080)
        angle += 92160;

    u8 quad = 0; //quadrant
    if (angle > 23040) {
        quad = 2;
        angle = 46080 - angle;
    }
    else if (angle < -23040) {
        quad = 3;
        angle = -46080 - angle;
    }

    s32 x = 156; //156 = 0.6073 * 256
    s32 y = 0;
    s32 x_new, y_new;

    for (u8 i = 0; i < 15; i++) {
        if (angle > 0) {
            x_new = x - (y >> i);
            y_new = y + (x >> i);
            x = x_new;
            y = y_new;
            angle -= tan_angle[i];
        } else if (angle < 0){
            x_new = x + (y >> i);
            y_new = y - (x >> i);
            x = x_new;
            y = y_new;
            angle += tan_angle[i];
        } else {
            break;
        }
    }

    if (quad == 2 || quad == 3)
        x = -x;

    return x;
}

/*
 * Brief   : the sin function.
 * Param   :
     - angle : input angle (magnified).
 * Return  : the magnified result.
 */
s32 my_sin(s32 angle)
{
    angle = angle % 92160;
    if (angle > 46080)
        angle -= 92160;
    else if (angle < -46080)
        angle += 92160;

    if (angle > 23040)
        angle = 46080 - angle;
    else if (angle < -23040)
        angle = -46080 - angle;

    s32 x = 156; //156 = 0.6073 * 256
    s32 y = 0;
    s32 x_new, y_new;

    for (u8 i = 0; i < 15; i++) {
        if (angle > 0) {
            x_new = x - (y >> i);
            y_new = y + (x >> i);
            x = x_new;
            y = y_new;
            angle -= tan_angle[i];
        } else if (angle < 0) {
            x_new = x + (y >> i);
            y_new = y - (x >> i);
            x = x_new;
            y = y_new;
            angle += tan_angle[i];
        } else {
            break;
        }
    }

    return y;
}

参考文献：

1、 https://blog.csdn.net/liyuanbhu/article/details/8458769

适用场景

✅ 适合使用 CORDIC 的场景：

  • 资源受限的 MCU（无 FPU）
  • 实时性要求高的 DSP 应用
  • 功耗敏感的应用

❌ 不适合使用 CORDIC 的场景：

  • 对精度要求极高
  • 有高性能 FPU 可用

总结

CORDIC 算法是一种优雅且高效的三角函数计算方法，特别适合资源受限的嵌入式系统。通过本文的实现，我们可以在没有浮点单元的平台上获得足够精度的三角函数计算结果。

关键要点：

  1. CORDIC 仅使用移位和加减法，无需乘法器
  2. 可以通过修改Q定点格式在精度和范围之间进行平衡

在实际项目中，可以根据具体需求调整迭代次数和定点格式，在精度和性能之间找到最佳平衡点。