如何随着时间的推移准确计算正弦波

Question

用例是为数字合成生成正弦波，因此，我们需要计算 sin(d t) 的所有值，其中：

t是一个整数，代表样本数。这是可变的。一小时 CD 音质的范围是 0 到 158,760,000。

d 是双精度的，表示角度的增量。这是恒定的。范围是：大于0，小于pi。

目标是使用传统的 int 和 double 数据类型实现高精度。性能并不重要。

天真的实现是：

double next()
{
    t++;
    return sin( ((double) t) * (d) );
}

但是，问题在于，当 t 增加时，准确度会降低，因为向“sin”函数提供了大量数字。

改进的版本如下：

double next()
{
    d_sum += d;
    if (d_sum >= (M_PI*2)) d_sum -= (M_PI*2);

    return sin(d_sum);
}

在这里，我确保将 0 到 2*pi 范围内的数字提供给“sin”函数。

但是，现在的问题是，当 d 很小时，会有很多小的添加，每次都会降低准确度。

这里的问题是如何提高准确性。

---

附录 1

“准确性降低，因为向“sin”函数提供了大量数字”：

#include <stdio.h>
#include <math.h>

#define TEST      (300000006.7846112)
#define TEST_MOD  (0.0463259891528704262050786960234519968548937998410258872449766)
#define SIN_TEST  (0.0463094209176730795999323058165987662490610492247070175523420)

int main()
{
    double a = sin(TEST);
    double b = sin(TEST_MOD);

    printf("a=%0.20f \n" , a);
    printf("diff=%0.20f \n" , a - SIN_TEST);
    printf("b=%0.20f \n" , b);
    printf("diff=%0.20f \n" , b - SIN_TEST);
    return 0;
}

输出：

a=0.04630944601888796475
diff=0.00000002510121488442
b=0.04630942091767308033
diff=0.00000000000000000000

Answer 1

对于超准确性，OP 有两个问题：

1. 将d 乘以n 并保持比double 更高的精度。这在下面的第一部分中得到了回答。

2. 执行期间的mod。简单的解决方案是使用度数，然后使用mod 360，这很容易做到。做大角度的2*π 是很棘手的，因为它需要2*π 的值，比(double) 2.0 * M_PI 多27 位的精度

---

用2个doubles代表d。

让我们假设 32 位 int 和 binary64 double。所以double 的准确度是 53 位。

0 <= n <= 158,760,000 大约是 2^27.2。由于double 可以连续准确地处理 53 位无符号整数，即 53-28 --> 25，因此任何只有 25 个有效位的double 都可以乘以n，并且仍然是精确的。

将 d 分段为 2 个 doubles dmsb,dlsb，最高 25 位和最低 28 位。

int exp;
double dmsb = frexp(d, &exp);  // exact result
dmsb = floor(dmsb * POW2_25);  // exact result
dmsb /= POW2_25;               // exact result
dmsb *= pow(2, exp);           // exact result
double dlsb = d - dmsb;        // exact result

那么dmsb*n 的每次乘法（或连续加法）都是精确的。 （这是重要的部分。）dlsb*n 只会在其最少的几个位上出错。

double next()
{
    d_sum_msb += dmsb;  // exact
    d_sum_lsb += dlsb;
    double angle = fmod(d_sum_msb, M_PI*2);  // exact
    angle += fmod(d_sum_lsb, M_PI*2);
    return sin(angle);
}

注意：fmod(x,y) 结果应与x,y 完全相同。

---

#include <stdio.h>
#include <math.h>

#define AS_n 158760000
double AS_d = 300000006.7846112 / AS_n;
double AS_d_sum_msb = 0.0;
double AS_d_sum_lsb = 0.0;
double AS_dmsb = 0.0;
double AS_dlsb = 0.0;

double next() {
  AS_d_sum_msb += AS_dmsb;  // exact
  AS_d_sum_lsb += AS_dlsb;
  double angle = fmod(AS_d_sum_msb, M_PI * 2);  // exact
  angle += fmod(AS_d_sum_lsb, M_PI * 2);
  return sin(angle);
}

#define POW2_25 (1U << 25)

int main(void) {
  int exp;
  AS_dmsb = frexp(AS_d, &exp);         // exact result
  AS_dmsb = floor(AS_dmsb * POW2_25);  // exact result
  AS_dmsb /= POW2_25;                  // exact result
  AS_dmsb *= pow(2, exp);              // exact result
  AS_dlsb = AS_d - AS_dmsb;            // exact result

  double y;
  for (long i = 0; i < AS_n; i++)
    y = next();
  printf("%.20f\n", y);
}

输出

0.04630942695385031893

---

使用度数

建议使用度数，因为360 度数是精确 周期，M_PI*2 弧度是近似值。 C 不能准确地表示 π。

如果 OP 仍想使用弧度，有关执行 π 的 mod 的进一步了解，请参阅Good to the Last Bit

Answer 2

以下代码将 sin() 函数的输入保持在一个小范围内，同时由于可能非常小的相位增量而在一定程度上减少了小的加法或减法的数量。

double next() {
    t0 += 1.0;
    d_sum = t0 * d;
    if ( d_sum > 2.0 * M_PI ) {
        t0 -= (( 2.0 * M_PI ) / d );
    }
    return (sin(d_sum));
}

Answer 3

将周期放大到 2^64，并使用整数运算进行乘法：

// constants:
double uint64Max = pow(2.0, 64.0);
double sinFactor = 2 * M_PI / (uint64Max);

// scale the period of the waveform up to 2^64
uint64_t multiplier = (uint64_t) floor(0.5 + uint64Max * d / (2.0 * M_PI));

// multiplication with index (implicitly modulo 2^64)
uint64_t x = i * multiplier;

// scale 2^64 down to 2π
double value = sin((double)x * sinFactor);

只要你的周期不是数十亿个样本，multiplier 的精度就足够了。

Answer 4

sin(x) = sin(x + 2N∙π)，所以问题可以归结为准确地找到一个小的数等于一个大的数数字 x 模 2π。

例如，–1.61059759 ≅ 256 mod 2π，您可以计算 sin(-1.61059759)，其精度高于 sin(256)

所以让我们选择一些整数来处理，256。首先找到等于 256 的幂的小数，模 2π：

// to be calculated once for a given frequency
// approximate hard-coded numbers for d = 1 below:
double modB = -1.61059759;  // = 256  mod (2π / d)
double modC =  2.37724612;  // = 256² mod (2π / d)
double modD = -0.89396887;  // = 256³ mod (2π / d)

然后将您的索引拆分为以 256 为基数的数字：

// split into a base 256 representation
int a = i         & 0xff;
int b = (i >> 8)  & 0xff;
int c = (i >> 16) & 0xff;
int d = (i >> 24) & 0xff;

您现在可以找到一个小得多的数字 x，它等于 i 模 2π/d

// use our smaller constants instead of the powers of 256
double x = a + modB * b + modC * c + modD * d;
double the_answer = sin(d * x);

对于 d 的不同值，您必须计算不同的值 modB、modC 和 modD，它们等于 256 的幂, 但取模 (2π / d)。您可以为这两个计算使用高精度库。

Answer 5

您可以尝试一种使用的方法，即快速傅立叶变换的一些实现。三角函数的值是根据以前的值和delta计算的。

Sin(A + d) = Sin(A) * Cos(d) + Cos(A) * Sin(d)

在这里，我们也必须存储和更新余弦值，并存储常数（对于给定的 delta）因子 Cos(d) 和 Sin(d)。

现在关于精度：小 d 的 cosine(d) 非常接近 1，因此存在精度损失的风险（数字中只有少数有效数字，例如 0.99999987）。为了克服这个问题，我们可以将常数因子存储为

dc = Cos(d) - 1 =  - 2 * Sin(d/2)^2
ds = Sin(d) 

使用其他公式更新当前值
（此处sa = Sin(A) 为当前值，ca = Cos(A) 为当前值）

ts = sa //remember last values
tc = ca
sa = sa * dc + ca * ds
ca = ca * dc - ts * ds
sa = sa + ts
ca = ca + tc

附：一些 FFT 实现会定期（每 K 步）通过 trig 更新 sa 和 ca 值。避免错误累积的函数。

示例结果。双打计算。

d=0.000125
800000000 iterations
finish angle 100000 radians

                             cos               sin
described method       -0.99936080743598  0.03574879796994 
Cos,Sin(100000)         -0.99936080743821  0.03574879797202
windows Calc           -0.9993608074382124518911354141448 
                            0.03574879797201650931647050069581