C语言中的无穷正弦生成

Question

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

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

目前我们只是在做：

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

出现两个问题/疑问：

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

Answer 1

引人入胜的线程。 Walter Mitty's answer 中提到的 Minsky 算法让我想起了在 Electronics & Wireless World 上发表并保留的画圆的方法。 （信用：https://www.electronicsworld.co.uk/magazines/）。我把它贴在这里是为了感兴趣。

但是，对于我自己的类似项目（用于音频合成），我使用查找表，其中包含足够多的点，线性插值足够准确（做数学！）

Answer 2

我想直接解决您代码中的嵌入式编程问题 - @0___________'s answer 是在微控制器上执行此操作的正确方法，我不会重蹈覆辙。

• 表示时间的变量从不是浮点数。如果您的增量不是 2 的幂，则错误总是会累积。即使是这样，最终您的增量将小于最小增量并且计时器将停止。始终使用整数表示时间。您可以选择一个足够大的整数来忽略翻转 - 表示毫秒的无符号 32 位整数需要 50 天才能翻转，而无符号 64 位整数需要超过 5 亿年。
• 在您不关心信号相位的情况下生成任何周期性信号不需要时间变量。相反，您可以保留一个在周期结束时重置为 0 的内部计数器。 （当您将 DMA 与查找表一起使用时，这正是您所做的 - 计数器是 DMA 控制器的下一次读取指针。）
• 每当您在微控制器中使用超越函数（例如正弦）时，您的第一个想法应该是“我可以为此使用查找表吗？”您无法在 4 GHz+ 多核处理器上使用现代操作系统来优化您的负载。您经常处理一个单线程，该线程将等待您的 200 MHz 微控制器使 FPU 脱离待机状态并执行近似算法。超越函数的成本很高。 LUTs 也有一定的成本，但如果您经常使用该功能，那么您很有可能会更喜欢 LUT 的权衡取舍。

Answer 3

有一种替代方法可以计算一系列正弦（和余弦）值，这些值增加了一些非常小的角度。它本质上是计算圆的 X 和 Y 坐标，然后将 Y 值除以某个常数得到正弦，再将 X 值除以相同的常数得到余弦。

如果您满足于生成“非常圆的椭圆”，则可以使用以下技巧，该技巧归因于 1960 年代的Marvin Minsky。它比计算正弦和余弦要快得多，尽管它在序列中引入了一个非常小的错误。这是 Hakmem 文档第 149 项的摘录。概述了明斯基圆算法。

---

第 149 项（明斯基）：循环算法 这是在点绘图显示器上绘制几乎圆形的优雅方法：

NEW X = OLD X - epsilon * OLD Y
NEW Y = OLD Y + epsilon * NEW(!) X

这会生成一个以原点为中心的非常圆的椭圆，其大小由初始点决定。 epsilon 决定了循环点的角速度，对偏心率有轻微的影响。如果 epsilon 是 2 的幂，那么我们甚至不需要乘法，更不用说平方根、正弦和余弦了！ “圆”将非常稳定，因为点很快就会变成周期性的。

当我试图在显示黑客中保存一个寄存器时，我错误地发明了圆形算法！ Ben Gurley 只用了大约六到七条指令就完成了一个惊人的显示技巧，这真是一个奇迹。但它基本上是面向线的。我突然想到有曲线会很令人兴奋，我试图用最少的指令来获得曲线显示技巧。

---

这里是 hakmem 的链接：http://inwap.com/pdp10/hbaker/hakmem/hacks.html

Answer 4

我认为可以使用模数，因为sin() 是周期性的。

那么你就不用担心这些问题了。

double time = 0;
long unsigned int timesteps = 0;
double sine;
void controll_loop()
{
  timesteps++;
  time += 30e-6;
  if( time > 1 )
  {
    time -= 1;
  }
  sine = sin( 2 * M_PI * 280 * time );
  ...
}

Answer 5

不要将正弦计算为时间的函数，而是维护一个正弦/余弦对并通过复数乘法将其推进。这不需要任何三角函数或查找表；只有四次乘法和偶尔的重新归一化：

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 没有二进制的精确表示。但是，通过按照建议重置生成器，我们将得到预期的精确周期性正弦。）

---

说明

有些人要求在 cmets 中解释为什么这样做。最简单的解释是使用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，因此我们只需要稍微修正一下就可以了。围绕1（一阶泰勒近似）扩展d 的表达式，我们得到代码中的公式：

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

TODO：要完全理解这个近似值的有效性，需要证明它可以比累积误差更快地补偿舍入误差，从而确定需要应用它的频率。

Answer 6

我对现有的答案感到相当震惊。你发现的第一个问题很容易解决，当你解决第一个问题时，下一个问题就神奇地消失了。

您需要对数学有基本的了解才能了解它是如何工作的。回想一下，sin(x+2pi) 在数学上只是sin(x)。当您的 sin(float) 实现切换到另一种算法时，您看到的时间会大幅增加，而您确实希望避免这种情况。

请记住，float 只有 6 个有效数字。 100000.0f*M_PI+x 将这 6 位数字用于 100000.0f*M_PI，因此 x 没有任何剩余。

因此，最简单的解决方案是自己跟踪x。在t=0，您将x 初始化为0.0f。每 30 个我们，您增加 x+= M_PI * 280 * 30e-06;。时间没有出现在这个公式中！最后，如果x>2*M_PI，则递减x-=2*M_PI;（因为sin(x)==sin(x-2*pi)

您现在有一个x，它很好地保持在0 到6.2834 的范围内，其中sin 速度很快，6 位精度都很有用。

Answer 7

使用查找表。您在与 Eugene Sh. 的讨论中的评论：

与正弦频率的小偏差（如 280.1Hz）是可以的。

在这种情况下，控制间隔为 30 µs，如果您有一个包含 119 个样本的表格，您会一遍又一遍地重复，您将得到一个 280.112 Hz 的正弦波。由于您有一个 12 位 DAC，如果您将其直接输出到 DAC，则只需要 119 * 2 = 238 个字节来存储它。如果您将其用作 cmets 中提到的进一步计算的输入，您可以根据需要将其存储为 float 或 double。在具有嵌入式静态 RAM 的 MCU 上，从内存加载最多只需要几个周期。

Answer 8

其他人的基于模的概念的变体怎么样：

int t = 0;
int divisor = 1000000;
void control_loop() {
    t += 30 * 280;
    if (t > divisor) t -= divisor;
    double sine = sin(2 * M_PI * t / (double)divisor));
    ...
}

它以整数计算模数，然后不会导致舍入误差。

Answer 9

如何生成可爱的正弦。

DAC 是 12 位的，因此您只有 4096 个级别。每个周期发送超过 4096 个样本是没有意义的。在现实生活中，生成高质量波形所需的样本要少得多。

1. 使用查找表创建 C 文件（使用您的 PC）。将输出重定向到文件 (https://helpdeskgeek.com/how-to/redirect-output-from-command-line-to-text-file/)。

#define STEP   ((2*M_PI) / 4096.0)

int main(void)
{
    double alpha = 0;

    printf("#include <stdint.h>\nconst uint16_t sine[4096] = {\n");
    for(int x = 0; x < 4096 / 16; x++)
    {
        for(int y = 0; y < 16; y++)
        {
            printf("%d, ", (int)(4095 * (sin(alpha) + 1.0) / 2.0));
            alpha += STEP;
        }
        printf("\n");
    }
    printf("};\n");
}

https://godbolt.org/z/e899d98oW

1. 配置定时器每秒触发溢出4096*280=1146880次。设置定时器以生成 DAC 触发事件。对于 180MHz 定时器时钟，它并不精确，频率为 279.906449045Hz。如果您需要更高的精度，请更改样本数以匹配您的定时器频率或/和更改定时器时钟频率（H7 定时器可以运行到 480MHz）

2. 将 DAC 配置为使用 DMA，并在触发事件时将值从步骤 1 中创建的查找表传输到 DAC。

3. 使用示波器享受美丽的正弦波。请注意，您的微控制器内核根本不会被加载。您将拥有它用于其他任务。如果要更改周期，只需重新配置计时器。您可以根据需要每秒执行多次。要重新配置定时器，请使用定时器 DMA 突发模式 - 这将在更新事件时自动重新加载 PSC 和 ARR 寄存器，而不会干扰生成的波形。

我知道这是高级 STM32 编程，需要寄存器级编程。我用它在我们的设备中生成复杂的波形。

这是正确的做法。没有控制循环，没有计算，没有核心负载。

Answer 10

如果您有几千字节的可用内存，则可以使用查找表完全消除此问题。

采样周期为 30 µs，2500 个样本的总持续时间为 75 ms。这正好等于 280 Hz 下 21 个周期的持续时间。

我还没有测试或编译以下代码，但它至少应该演示该方法：

double sin2500() {
  static double *table = NULL;
  static int n = 2499;
  if (!table) {
    table = malloc(2500 * sizeof(double));
    for (int i=0; i<2500; i++) table[i] = sin(2 * M_PI * 280 * i * 30e-06);
  }
  n = (n+1) % 2500;
  return table[n];
}

Answer 11

函数可以改写为

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

这与问题中的代码完全相同，但问题完全相同。但现在可以简化：

280 * 30e-6 = 280 * 30 / 1000000 = 21 / 2500 = 8.4e-3

这意味着当n 达到 2500 时，您已经输出了正好 21 个周期的正弦波。这意味着您可以将 n 设置回 0。 结果代码是：

int n;
void control_loop() {
    n += 1;
    if (n == 2500)
        n = 0;
    double sine = sin(2 * M_PI * 8.4e-3 * n);
    ...
}

只要您的代码可以运行 21 个周期而不会出现问题，它就会永远运行而不会出现问题。

Answer 12

正如在某些 cmets 中所指出的，time 值会随着时间不断增长。这带来了两个问题：

1. sin 函数可能必须在内部执行模数才能将内部值置于支持的范围内。
2. 由于增加了更高的数字，时间分辨率会随着值的增加而变得越来越差。

进行以下更改应该会提高性能：

double time;
void control_loop() {
    time += 30.0e-6;
    if((1.0/280.0) < time)
    {
        time -= 1.0/280.0;
    }
    
    double sine = sin(2 * M_PI * 280 * time);
    ...
}

请注意，一旦进行此更改，您将不再有时间变量。