在C++中创建正弦查找表

22

如何将以下伪代码重写为C++?

real array sine_table[-1000..1000]
    for x from -1000 to 1000
        sine_table[x] := sine(pi * x / 1000)
我需要创建一个正弦表查找表。

2
你需要在哪个任务上寻求帮助?是调用 sin 函数吗?还是应对需要负数索引的情况?或者是声明一个数组?或者是在 C++ 中编写 for 循环?亦或是找到代表实数类型的名称?还是全部都需要帮助? - Steve Jessop
6
当你完成查找表时,你可能需要对其进行性能测试,并验证它是否真的比调用sin()(或sinf())更快。现代CPU的运行速度比现代RAM要快得多,因此你可能会发现在查找表中查找结果所需的时间比重新计算结果要长(因为计算结果可以完全在CPU上完成)。 - Jeremy Friesner
不仅如此,你还有更多的代码,需要从CPU缓存中清除一些行。你需要在实际例程中进行性能测试,而不仅仅是在测试平台上。 - Loren Pechtel
4
这可能是针对嵌入式系统或FPGA,其中查找速度比计算 x — (x^3)/(3!) + (x^5)/(5!) — (x^7)/(7!) + ... 更快。 - Mike DeSimone
所选答案提供了一个解决方案,用于问题“如何在不使用语言内置的sin方法的情况下计算正弦函数的值”。在这种情况下,泰勒是最著名的软件实现算法,但是计算机ALU、计算器和其他数字系统通常通过硬件模拟正弦函数(以及许多其他函数),实现CORDIC示例)。它比完整的查找表具有更轻的占用空间。 - mins
6个回答

38

通过仅存储第一象限(即x在[0,pi/2])的值,您可以将表格的大小减小为原来的25%。

要做到这一点,您的查找程序只需要使用简单的三角恒等式将所有x的值映射到第一象限:

  • sin(x) = - sin(-x),从第四象限映射到第一象限
  • sin(x) = sin(pi - x),从第二象限映射到第一象限

要从第三象限映射到第一象限,请应用两个恒等式,即sin(x) = - sin (pi + x)

这种策略是否有所帮助取决于内存使用对您的案例有多重要。但是,存储四倍于您需要的值只是为了避免在查找期间进行比较和减法似乎是浪费的。

我赞同Jeremy的建议,即测量构建表格是否比仅使用std :: sin()更好。即使使用原始大表格,您每次查找表格时也必须花费周期来将参数转换为最接近pi / 1000的增量,并且在此过程中会失去一些精度。

如果您真的想要以速度为代价交换精度,则可以尝试仅使用Taylor级数展开的前几个项来逼近sin()函数。

  • sin(x) = x - x^3/3! + x^5/5! ...,其中^表示乘方运算,!表示阶乘。

当然,为了效率,您应该预先计算阶乘并利用较低次幂的x来计算更高次幂,例如在计算x^5时使用x^3。

最后一个要点是,上面截断的Taylor级数对于接近零的值更精确,因此在计算近似正弦之前将其映射到第一或第四象限仍然是有价值的。

补充说明: 基于两个观察到的潜在改进:
1.如果您可以在第一个八分之一[0,pi/4]中计算出正弦和余弦,则可以计算任何三角函数
2.以零为中心的Taylor级数展开在接近零时更精确

如果你决定使用截断的泰勒级数,那么你可以通过映射到正弦或余弦来改善精度(或者使用更少的项以获得类似的精度),从而使角度在区间[0,pi/4]内,使用诸如sin(x) = cos(pi/2-x)和cos(x) = sin(pi/2-x)之类的恒等式,除了上面提到的公式(例如,如果x> pi/4,并且你已经映射到第一象限)。
或者,如果你决定同时使用正弦和余弦的表格查找,你可以通过额外的比较和减法将查找映射到更小的范围,只覆盖[0,pi/4]范围的两个较小的表格。然后,你可以使用更少的内存来存储表格,或者使用相同的内存但提供更细的粒度和精度。

5
此外,Quarter-Sine表可以通过类似的象限重新映射方法用于寻找余弦值。 - Mike DeSimone
3
第一个恒等式有一个拼写错误,应该是 sin(-x) = -sin(x) - Ben Voigt
谢谢你发现这个问题,Ben。我已经相应地更正了我的回答。 - Alex Blakemore
1
“你可以通过仅存储第一象限的值,将表格的大小减小到原始大小的25%。”无论如何,要获得[-pi,+pi]的实际值,级数的顺序必须至少为10,级数6为[-pi/2,+pi/2]提供了良好的近似值(图表)。 - mins

5

还有一点需要注意:调用三角函数的代价很高。如果你想要为正弦函数准备具有恒定步长的查找表,你可以节省计算时间,但会损失一些潜在的精度。

假设你的最小步长为“a”。也就是说,你需要 sin(a),sin(2a),sin(3a) 等等。

那么你可以使用以下技巧:首先计算 sin(a) 和 cos(a)。然后对于每个连续的步骤,使用以下三角函数等式:

  • sin([n+1] * a) = sin(n*a) * cos(a) + cos(n*a) * sin(a)
  • cos([n+1] * a) = cos(n*a) * cos(a) - sin(n*a) * sin(a)

这种方法的缺点是在此过程中会积累舍入误差。


5
long double sine_table[2001];
for (int index = 0; index < 2001; index++)
{
    sine_table[index] = std::sin(PI * (index - 1000) / 1000.0);
}

  1. 别忘了 #include <cmath>
  2. 实际上这样只会创建少一个项目;在伪代码中相当于 sine_table[-1000..999]
- Mike DeSimone

3


double table[1000] = {0};
for (int i = 1; i <= 1000; i++)
{
    sine_table[i-1] = std::sin(PI * i/ 1000.0);
}

double getSineValue(int multipleOfPi){ 如果(multipleOfPi == 0) 返回0.0; int 符号 = 1; 如果(multipleOfPi < 0){ 符号 = -1;
} 返回符号sine_table[符号multipleOfPi - 1]; }

您可以通过技巧sin(pi/2 +/- angle) = +/- cos(angle)将数组长度缩小到500。 所以存储0到pi/4的sin和cos。 我不记得从头脑中提取的内容,但它加速了我的程序。


1
一个注意点。getSineValue()的参数名称可能会误导。它实际上代表Pi/1000的倍数,因此可以称为thousandthsOfPi。不是要卖弄学识,但这很容易让某人混淆getSineValue(500)表示sin(pi/2)。 - Alex Blakemore
我对将for循环从0开始优化到1感到惊讶。 :) 请注意,该循环存在一个错误,它会写入table [1000],这不是数组中的有效索引(它超出了范围)。 - Jeremy Friesner
很抱歉出现了循环错误。实际上我的意思是将1到1000存储为0是不明显的情况。我已经编辑了代码以进行更正。 - Master Yoda
我认为还不太对...例如,现在sine_table [0]的值将是sin(PI * 1 / 1000.0)的值,而sin(0)的正确值应该是零。 - Jeremy Friesner
这就是为什么它在这里:返回 signsine_table[signmultipleOfPi - 1]; - Master Yoda

2
你需要使用<cmath>中的std::sin()函数。

2

以下是从某本书或其他来源获得的另一种近似方法:

streamin ramp;
streamout sine;

float x,rect,k,i,j;

x = ramp -0.5;

rect = x * (1 - x < 0 & 2);
k = (rect + 0.42493299) *(rect -0.5) * (rect - 0.92493302) ;
i = 0.436501 + (rect * (rect + 1.05802));
j = 1.21551 + (rect * (rect - 2.0580201));
sine = i*j*k*60.252201*x;

完整讨论请点击此处:http://synthmaker.co.uk/forum/viewtopic.php?f=4&t=6457&st=0&sk=t&sd=a 我认为您知道,使用除法比乘以小数要慢得多,/5总是比*0.2慢。
这只是一个近似值。
另外:
streamin ramp;
streamin x;  // 1.5 = Saw   3.142 = Sin   4.5 = SawSin
streamout sine;
float saw,saw2;
saw = (ramp * 2 - 1) * x;
saw2 = saw * saw;

sine = -0.166667 + saw2 * (0.00833333 + saw2 * (-0.000198409 + saw2 * (2.7526e-006+saw2 * -2.39e-008)));
sine = saw * (1+ saw2 * sine);

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接