正弦函数频率突然跳跃

Question

我正在为数字信号处理编写信号源，但发现了奇怪的行为。以下是代码：

    float complex *output = malloc(sizeof(float complex) * 48000);
    FILE *f = fopen("test_signal.cf32", "wb");
    int64_t freq = -3355;
    float phase_increment = 2 * (float) M_PI * (float) freq / 48000;
    float phase = 0.0F;
    for (int i = 0; i < 150 * 5; i++) {
        for (int j = 0; j < 9600; j++) {
            output[j] = cosf(phase) + sinf(phase) * I;
            phase += phase_increment;
        }
        fwrite(output, sizeof(float complex), 9600, f);
    }
    fclose(f);

它将创建一个复杂信号的文件，该信号从中心向左偏移-3355Hz。因此，当我在此文件上运行FFT时，我期望在-3355Hz处获得信号。由于量化，可能会在该数字周围有一些次要频率分布。但实际情况是我得到了以下结果：

大约在50秒左右，频率出现了相当显著的跳变（~2000hz）。

有人知道原因吗？

我的实验表明：

• 使用双精度浮点数的phase有所帮助
• 将参数从-2pi到2pi进行减少有所帮助
• 苹果M1和树莓派有相同的跳变
• phase参数似乎没有溢出
• 有趣的阅读https://randomascii.wordpress.com/2014/10/09/intel-underestimates-error-bounds-by-1-3-quintillion/，但不太可能导致频率突然跳变

- dernasherbrezon

---

在 output[j] = cosf(phase) + sinf(phase) * I; 中，这个大写的 I 是不是和小写的 i 不同的变量？如果是的话，那这个 I 变量是什么？ - user8849929

@Sedenion，它是从complex.h导入的虚数单位。 - tstanisl

在https://dev59.com/FVEG5IYBdhLWcg3wPXtF上讨论了无限正弦波的生成主题。这不是你问题的答案，但你可能会发现它有用。 - tstanisl

sin() 的 MS man page 表示：“如果 x 大于或等于 263，或小于或等于 -263，则结果会出现精度损失。” - Weather Vane

Answer 1

一旦“phase”变得足够大，对其的增量将通过浮点舍入进行量化，可能会导致舍入到最近的舍入到偶数尾数无法平衡，而是意味着您始终将相位推进得更少。请参见wikipedia's article on single-precision IEEE float 一个基于10进制的类比，值被限制为4个有效数字，例如“101.2 + 0.111”只能上升到“101.3”，而不是“101.311”。（计算机使用二进制浮点数，因此实际上可以精确表示诸如101.125之类的内容，其中101.1则不行。）

我怀疑这就是发生的事情。您不会溢出到“+Inf”，但是考虑到float的有限相对精度，将小数加到大数中最终将不会移动它：与“phase + phase_increment”最接近的可表示浮点数将是原始的“phase”。

一个快速测试这个假设的方法是使用double phase（不进行任何其他更改），看看是否移动了您获得频率步骤的点。（可能大大超过您生成的结束点）。

（哦，刚注意到您已经这样做了，并且它确实有所帮助。如果您的意思是它完全解决了问题，那么可能就是这样。）

您仍然可以保留单精度的sinf和cosf，尽管更改它们是另一种尝试的方法：正如@Weather Vane指出的那样，某些sin实现在范围缩减中失去了显着的精度。但您会期望这是一个更嘈杂的信号，而不是频率的一致变化（需要相位的比例因子）。

您还可以让它在float下运行更长时间，看看是否会再次降低频率，或者在DC（常数相位，频率= 0）处发生变化。

使用 float，手动展开循环，使得你有两个计数器，它们的偏移量为 phase_increment，并且每个计数器都增加 2*phase_increment。这样，phase 相对于正在添加的内容的相对幅度就不会那么极端了。它们仍然会达到相同的绝对幅度，尽管相对于原始的 phase_increment 来说很大，因此仍然会有一些舍入误差。

---

我不知道生成复杂正弦波的更好策略是否存在，我只是试图解释您所看到的行为。

对于此方法的性能，如果您的编译器未将您的调用优化为同时计算两个值，则可能需要使用sincosf。 （并希望通过调用SIMD sincosf 自动向量化。）

- Peter Cordes

---

明白了，谢谢！顺便说一下：对于这段代码，-O2会生成sincosf。 - dernasherbrezon

Answer 2

你所观察到的是浮点数精度问题。随着数量级的增加，浮点数的精度会降低。

幸运的是，对于相位或类似的角度值，有一个简单的解决方法，将值包裹在2π周围，这样数量级就永远不会太大。

在你的phase += phase_increment;后添加以下代码，看看是否有帮助。

if( phase < (float)( 2.0 * M_PI ) )
    continue;
phase -= (float)( 2.0 * M_PI );

- Soonts

---

降维运算很好用：if( phase < -(float)( 2.0 * M_PI ) ) { phase += (float)( 2.0 * M_PI ); } if( phase > (float)( 2.0 * M_PI ) ) { phase -= (float)( 2.0 * M_PI ); } - dernasherbrezon