C语言中无休止的正弦生成 [英] Endless sine generation in C

问题描述

我正在进行一个项目，该项目将计算正弦波作为控制回路的输入。

正弦波的频率为280赫兹，控制回路每30微秒运行一次，Arm Cortex-M7的所有内容都用C语言编写。

目前我们只是在做：

double time;
void control_loop() {
    time += 30e-6;
    double sine = sin(2 * M_PI * 280 * time);
    ...
}

出现两个问题/问题：

1. 长时间运行，time变大。突然之间，正弦函数的计算时间急剧增加(见下图)。这是为什么？这些功能通常是如何实现的？由于速度对我们来说是一个巨大的因素，有没有办法绕过这一点(而不会造成明显的精度损失)？我们使用的是Math.h(手臂GCC)中的Sin。
2. 一般情况下，我如何处理时间？在长时间运行时，变量time不可避免地会达到双精度极限。即使使用计数器time = counter++ * 30e-6;也只能改善这一点，但不能解决它。因为我肯定不是第一个想要在很长一段时间内产生正弦波的人，所以肯定有一些想法/论文/…关于如何快速而准确地实现这一点。

推荐答案

不计算作为时间函数的正弦，而是保持正弦/余弦对并通过复数乘法推进它。这不需要任何三角函数或查找表；只需要四次乘法和偶尔的重新标准化：

static const double a = 2 * M_PI * 280 * 30e-6;
static const double dx = cos(a);
static const double dy = sin(a);
double x = 1, y = 0; // complex x + iy
int counter = 0;

void control_loop() {
    double xx = dx*x - dy*y;
    double yy = dx*y + dy*x;
    x = xx, y = yy;

    // renormalize once in a while, based on
    // https://www.gamedev.net/forums/topic.asp?topic_id=278849
    if((counter++ & 0xff) == 0) {
        double d = 1 - (x*x + y*y - 1)/2;
        x *= d, y *= d;
    }

    double sine = y; // this is your sine
}

如果需要，可以通过重新计算dx、dy来调整频率。

此外，这里的所有操作都可以相当轻松地在定点上完成。

---

合理性

作为@user3386109 points out below (+1)，280 * 30e-6 = 21 / 2500是有理数，因此正弦应该正好在2500个样本之后循环。我们可以通过每2500次(或5000次、或10000次等)重置我们的生成器(x=1,y=0)来将此方法与他们的方法结合起来。这将消除重整化的需要，并消除任何长期的相位误差。

(从技术上讲，任何浮点数都是二进有理数。然而，280 * 30e-6没有精确的二进制表示。然而，通过按照建议重置发电机，我们将获得预期的精确周期正弦。)

---

说明

有些人要求在评论中解释为什么会这样做。最简单的解释是使用angle sum trigonometric identities：

xx = cos((n+1)*a) = cos(n*a)*cos(a) - sin(n*a)*sin(a) = x*dx - y*dy
yy = sin((n+1)*a) = sin(n*a)*cos(a) + cos(n*a)*sin(a) = y*dx + x*dy

正确之后是归纳法。

如果我们根据Euler's formula将这些正弦/余弦对视为复数，则这实质上就是De Moivre's formula。

更有洞察力的方法可能是几何地看待它。复数乘以exp(ia)等同于旋转a弧度。因此，通过重复乘以dx + idy = exp(ia)，我们沿着单位圆递增地旋转起点1 + 0i。根据欧拉公式，y坐标是当前相位的正弦。

正常化

虽然相位随着每次迭代而继续前进，但由于舍入误差，x + iy的量级(又名范数)偏离了1。然而，我们感兴趣的是产生幅度1的正弦，因此我们需要对x + iy进行归一化以补偿数值漂移。当然，直截了当的方法是用它自己的标准来划分：

double d = 1/sqrt(x*x + y*y);
x *= d, y *= d;

这需要计算平方根的倒数。即使我们每X次迭代才标准化一次，避免它仍然是很酷的。幸运的是，|x + iy|已经接近1，因此我们只需稍作修正就可以将其遏制住。展开d围绕1(一阶泰勒近似)的表达式，我们得到代码中的公式：

d = 1 - (x*x + y*y - 1)/2

TODO：要完全理解此近似的有效性，需要证明它补偿舍入误差的速度比它们累积的速度更快--从而得到需要应用它的频率的界限。

这篇关于C语言中无休止的正弦生成的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！