c - 正弦函数频率突然跳变

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 中使用 double 有帮助
• 将参数从 -2pi 减少到 2pi 有帮助
• 苹果M1和树莓派跳转相同
• phase 参数似乎没有溢出
• 有趣的阅读 https://randomascii.wordpress.com/2014/10/09/intel-underestimates-error-bounds-by-1-3-quintillion/但不太可能导致频率突然跳跃

Answer 1

一旦 phase 变得足够大，它的增量就会通过浮点舍入进行量化，也许足以使四舍五入到最接近尾数的四舍五入不能平衡，取而代之的是意味着您始终少推进阶段。参见 wikipedia's article on single-precision IEEE float

值限于 4 位有效数字的以 10 为底的类比将是 101.2 + 0.111 只会上升到 101.3，而不是 101.311。 (计算机使用二进制 float ，所以实际上像 101.125 这样的东西可以精确表示，而 101.1 则不能。)

我怀疑这就是正在发生的事情。您不会溢出到 +Inf，但考虑到 float 的相对精度有限，将一个小数字添加到一个巨大的数字最终根本不会移动它:最接近 phase + phase_increment 的可表示 float 将是原始 phase。

检验这个假设的一种快速方法是使用double phase(不做任何其他更改)，看看它是否移动了您获得频率的点脚步。 (要大得多，可能会超过您生成的内容)。

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

您仍然可以保留单精度 sinf 和 cosf，尽管更改它们是另一回事:正如@Weather Vane 指出的那样，一些 sin 实现在缩小范围方面失去了显着的精度。但是您会认为这是一个噪声较大的信号，而不是频率的一致变化(这需要相位的比例因子。)

您还可以使用 float 让它运行更长时间，看看您是否在频率上又降低了一步，或者一直改变到 DC(恒定相位，频率 = 0)。

或者使用 float，手动展开，这样您就有两个计数器偏移 phase_increment，并且每个计数器递增 2*phase_increment。因此，phase 与添加内容的相对大小不会变得如此极端。不过，它们仍然会达到相同的绝对幅度，相对于原始 phase_increment 而言较大，因此仍然会有一些舍入。

---

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

为了这个方法的性能，如果你的编译器没有优化你的调用，你可能希望 sincosf 同时计算两个值。 (并希望通过调用 SIMD sincosf 进行自动矢量化。)