c语言 - 如何精确计算随时间变化的正弦波 - c - How to compute sine wave with accuracy over the time 小贝子编程网

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

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

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

在这里，我确保为"sin"函数提供从 0 到 2*pi 范围内的数字。

但是，现在的问题是，当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

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

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

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

现在关于精度：小d的余弦(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

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

例如，–1.61059759 ≅ 256mod2π，您可以比sin(256)更精确地计算sin(-1.61059759)

因此，让我们选择一些整数来使用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)。您可以使用高精度库进行这几次计算。

将周期放大到 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的精度就足够好了。

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

对于超精度，OP 有 2 个问题：

将d乘以n并保持比double更高的精度。这在下面的第一部分中得到了回答。
执行该时期的mod。简单的解决方案是使用度数，然后mod 360，很容易做到。要2*π大角度是很棘手的，因为它需要2*π的值，并且精度比(double) 2.0 * M_PI多约27位

使用 2doubles 表示d。

让我们假设 32 位int和二进制 64double. 所以double有53位的精度。

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

将d分割成 2 个doublesdmsb,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("%.20fn", y);
}

输出

0.04630942695385031893

使用学位

建议使用度，因为360度是确切的周期，M_PI*2弧度是近似值。 C 不能完全表示π。

如果 OP 仍然想使用弧度，有关执行 π mod 的进一步见解，请参阅好到最后位

c语言 - 如何精确计算随时间变化的正弦波

相关内容

最新更新

热门标签：