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

Question

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

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

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

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

在这里，我们也必须存储和更新余弦值并存储常数（对于给定的增量）因子 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

PS 一些 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           

Answer 2

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);

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

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

For hyper accuracy, OP has 2 problems:

1. multiplying d by n and maintaining more precision than double. That is answered in the first part below.

2. Performing a mod of the period. The simple solution is to use degrees and then mod 360, easy enough to do exactly. To do 2*π of large angles is tricky as it needs a value of 2*π with about 27 more bits of accuracy than (double) 2.0 * M_PI

---

Use 2 doubles to represent d.

Let us assume 32-bit int and binary64 double. So double has 53-bits of accuracy.

0 <= n <= 158,760,000 which is about 2^27.2. Since double can handle 53-bit unsigned integers continuously and exactly, 53-28 --> 25, any double with only 25 significant bits can be multiplied by n and still be exact.

Segment d into 2 doubles dmsb,dlsb, the 25-most significant digits and the 28- least.

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

Then each multiplication (or successive addition) of dmsb*n will be exact. (this is the important part.) dlsb*n will only error in its least few bits.

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);
}

Note: fmod(x,y) results are expected to be exact give exact 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);
}

Output

0.04630942695385031893

---

Use degrees

Recommend using degrees as 360 degrees is the exact period and M_PI*2 radians is an approximation. C cannot represent π exactly.

If OP still wants to use radians, for further insight on performing the mod of π, see Good to the Last Bit

Answer 5

以下代码将 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));
}