反函数atan2的速度优化

前言

在计算雷达方位角和俯仰角时有时候需要使用atan(),asin(),acos()这些反三角函数，但是使用这些反三角函数，C语言的math.h和arm_math.h的头文件中都有对应的库函数可以提供用户使用。尽管使用系统自带的库函数简单并且误差很小，但是自带的三角函数的计算时间在某些重视速度并不需要特别高精确性的系统中是不能接受的，甚至可能导致算法的时间超过约束的时间，所以在这个情况下必须使用一些方法来提高计算的速度，但是这种提高往往是以牺牲计算的精度或者空间来换取运行的时间。

atan和atan2之间的关系

在实际中计算雷达的角度时，我使用了一个atan2()函数，来计算FFT结果的相位角度。有些同学可能没有用过atan2()这个函数，其实这个函数和atan()是差不多，只不过atan2()函数可以根据输入参数来确定角度的象限。两者具体的关系如下：

a t a n 2 (y, x) = {\begin{aligned} a c t a n (\frac{y}{x}) & x > 0 \\ a c t a n (\frac{y}{x}) + π & y \geq 0, x < 0 \\ a c t a n (\frac{y}{x}) - π & y < 0, x < 0 \\ + \frac{π}{2} & y > 0, x = 0 \\ - \frac{π}{2} & y < 0, x = 0 \\ u n d e f i e d & y = 0, x = 0 \end{aligned}

其实两者是相同，arctan(x)返回的是 $(- \frac{π}{2}, \frac{π}{2}]$ ， arctan2(y,x)返回的范围是 $(- π, π]$ .

actan2的快速进行计算

使用actan(x)的近似的近似的计算的一个理论根据是利用三角反函数的 泰勒展开，~~众所皆知~~actan(x)的泰勒展开式如下所示：

a c t a n (x) = x - \frac{x^{3}}{3} + \frac{x^{5}}{5} - \frac{x^{7}}{7} + . . . + \frac{x^{n}}{n} + . . .

根据泰勒展开，根据你需要精度，都可以无限的逼近求解出actan(x)的值，但是在x接近1时，泰勒收敛的速度并不是特别的迅速，所以可以利用一些快速的近似算法来进行求解。

线性近似，最大的误差为 $0.07 r a d = 4^{\circ}$

a c t a n (x) \approx \frac{π}{4} x - 1 \leq x \leq 1

二阶近似，最大的误差为 $0.0053 r a d = {0.3}^{\circ}$

$a c t a n (x) \approx \frac{π}{4} x + 0.285 x (1 - | x |) - 1 \leq x \leq 1$
搜索更好的系数，通过系数匹配后，最大的近似误差为 $0.0038 r a d = {0.22}^{\circ}$

$α x^{3} + β x$ 近似，最大近似误差 $0.005 r a d = {0.29}^{\circ}$

a c t a n (x) \approx \frac{π}{4} x + x (0.186982 - 0.191942 x^{2}) - 1 \leq x \leq 1

$x (x - 1) (α x - β)$ 形式的三阶近似，最大的近似误差 $0.0015 r a d = {0.086}^{\circ}$

a c t a n (x) \approx \frac{π}{4} x - x (| x | - 1) (0.2447 + 0.0663 | x |) - 1 \leq x \leq 1

$x / (1 + β x^{2})$ 形式的近似，最大近似误差 $0.0047 r a d = {0.27}^{\circ}$

a c t a n (x) \approx \frac{x}{1 + 0.28086 x^{2}} - 1 \leq x \leq 1

一个小例子，近似化处理

这是在一个固定的区间[-1,1]之间的快速近似，在实际中区间肯定会存在其他位置的取值的情况，我在网上查询资料时，发现这样一个算法，大体的实现过程如下。

/***************************************************************************************
@Copyright    : 
@Author       : Krone
@Data         : Do not edit
@LastEditor   : 
@LastData     : 
@Describe     : actan2近似计算
*****************************************************************************************/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "math.h"
#include "stdbool.h"

#define max(a,b) (((a) > (b)) ? (a) : (b))
#define min(a,b) (((a) < (b)) ? (a) : (b))

#define ESP                     0.00001f			 /* 测试的误差 */
#define STEP                    0.01f			 /* 误差范围 */
#define PI                      3.14159274f 
/***************************************************************************************
@funticon name       : easy_actan2
@berif               : 简单计算actan2函数
@author              : Krone
@data                : 
@note                : Need note condition
@param                {float} x --- 输出的分子
@param                {float} y --- 输出的分母
@result              : 
@testResult          : Function self test result ok or fault
*****************************************************************************************/
static float easy_actan2(float dy, float dx)
{   
    if(dy > 0.0f && dx == 0.0f)
    {
        return PI/2;
    }
    if (dy < 0.0f && dx == 0.0f)
    {
        return -PI/2;
    }
    if (dy == 0.0f && dx == 0.0f)
    {
        printf("undefined\n");
    }
    
    float ax = fabs(dx), ay = fabs(dy);
    float temp1 = min(ax, ay)/max(ax, ay);
    float temp2 = temp1*temp1;
    float result = ((-0.0464964749 * temp2 + 0.15931422) * temp2 - 0.327622764) * temp2 * temp1 + temp1;

	if (ay > ax)
	{
		result = (PI / 2) - result;
	}
	if (dx < 0.0f)
	{
		result = PI - result;
	}
	if (dy < 0.0f)
	{
		result = (-result);
	}
    return result;
}

int main(void)
{
    /* 测试代码，测试代码和原本的库函数进行对比 */
    float angle = -PI;			/* 单位圆的角度 */ 
    float x_value = 0.0f;
    float y_value = 0.0f;
    bool is_success = true;
	float max_esp = 0.0f;
	float esp = 0.0f;

    for (angle; angle < PI- STEP; angle+=STEP)
    {
		y_value = sin(angle);
		x_value = cos(angle);
		esp = fabs(atan2(y_value, x_value) - easy_actan2(y_value, x_value));
		if ((esp - max_esp) > 1e-8)
		{
			max_esp = esp;
		}
        if(esp > ESP)
        {
            is_success = false;
        }
    }
    if (is_success == false)
    {
        printf("超过误差范围\n");
    }
    else
    {
        printf("测试没有超过范围");
    }
	printf("最大误差%f", esp);
	return 0;
}