c - 如何改进 sin/cos lut 函数？

Question

我正在尝试为我的 MCU 优化 sin/cos 以计算地理距离。这部分公式特别用到了三角函数:

double e = (MyTan( lat2 / 2 + quarter_pi ) / MyTan( lat1 / 2 + quarter_pi ));

所以我尝试为 -PI 到 PI 构建自己的 sin/cos 查找表，如下所示:

#define PARTPERDEGREE 10
double mysinlut[PARTPERDEGREE * 90 + 1];
double mycoslut[PARTPERDEGREE * 90 + 1];
void MySinCosCreate()
{
    int i;
    double angle, angleinc;

    // Each degree also divided into 10 parts
    angleinc = (M_PI / 180) / PARTPERDEGREE;
    for (i = 0, angle = 0.0; i <= (PARTPERDEGREE * 90 + 1); ++i, angle += angleinc)
    {
        mysinlut[i] = sin(angle);
    }

    angleinc = (M_PI / 180) / PARTPERDEGREE;
    for (i = 0, angle = 0.0; i <= (PARTPERDEGREE * 90 + 1); ++i, angle += angleinc)
    {
        mycoslut[i] = cos(angle);
    }
}



double MySin(double rad)
{
   int ix;
   int sign = 1;


   if(rad > (M_PI / 2))
        rad = M_PI / 2 - (rad - M_PI / 2);

   if(rad < -(M_PI / 2))
       rad = -M_PI / 2 - (rad + M_PI / 2);

   if(rad < 0)
   {
        sign = -1;
        rad *= -1;
   }

   ix = (rad * 180) / M_PI * PARTPERDEGREE;

   return sign * mysinlut[ix];
}

double MyCos(double rad)
{
   int ix;
   int sign = 1;


    if(rad > M_PI / 2)
    {
        rad = M_PI / 2 - (rad - M_PI / 2);
        sign = -1;
    }
    else if(rad < -(M_PI / 2))
    {
        rad = M_PI / 2 + (rad + M_PI / 2);
        sign = -1;
    }
    else if(rad > -M_PI / 2 && rad < M_PI / 2)
    {   
        rad = abs(rad);
        sign = 1;
    }

    ix = (rad * 180) / M_PI * PARTPERDEGREE;

    return sign * mycoslut[ix];
}

double MyTan(double rad)
{
    return MySin(rad) / MyCos(rad);
}

你可以看到表格的分辨率是每度10个部分。我可以稍微增加一点，但帮助不大，看起来我需要一些插值。谁能建议对我的功能进行一些实际改进以获得更好的结果？下面是 e 的 234 个不同结果的图表。蓝色系列具有理想的 sin/cos，红色系列来自 LUT。

Answer 1

您的查找表似乎太粗糙了。如果您不能使您的表更精细，使用导数来近似值应该可以得到更好的结果。我们有

sin (x+h) ≈ sin x + h*cos x
cos (x+h) ≈ cos x - h*sin x

对于小h。 (您可以使用更高的导数或使用您的(计算的)角度所在的表中的两个值来获得更好的近似值，但这需要更长的时间，而且我收集速度是 LUT 的首要原因。)

所以在你归一化角度之后，用

ix = (rad * 180) / M_PI * PARTPERDEGREE;

使用

double h = rad - ix*angleinc;
return sign*(mysinlut[ix] + h*mycoslut[ix]);

回复

return sign*(mycoslut[ix] - h*mysinlut[ix]);

这不应该太慢，并且应该在 LUT 点之间提供明显更好的近似值。