正弦函数频率突然跳变

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与raspberry pi跳转相同
phase 参数似乎没有溢出
有趣的阅读https://randomascii.wordpress.com/2014/10/09/intel-underestimates-error-bounds-by-1-3-quintillion/但不太可能导致频率突然跳跃

Answer 1

一旦 phase 变得足够大，它的增量将通过 floating-point 舍入量化，也许足以舍入到最接近的 tie-break 甚至尾数的舍入不平衡，而是意味着您始终少推进阶段。参见 wikipedia's article on single-precision IEEE float

一个以 10 为基数的类比，其值限制为 4 个有效数字 101.2 + 0.111 只会上升到 101.3，而不是 101.311。（计算机使用二进制浮点数，所以实际上像 101.125 这样的东西可以精确表示，而 101.1 则不能。）

我怀疑这就是正在发生的事情。你不会溢出到 +Inf，但鉴于 float 的相对精度有限，将一个小数字添加到一个巨大的数字最终根本不会移动它：最接近 [=16 的可表示浮点数=] 将是原来的 phase.

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

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

您仍然可以保留 single-precision sinf 和 cosf，尽管更改它们是另一回事：正如@Weather Vane 指出的那样，一些 sin 实现在 range-reduction 中失去了显着的精度。但是您会认为这是一个噪声较大的信号，而不是频率的一致变化（这需要相位的比例因子。）

您还可以通过 float 让它运行更长，看看您是否可以进一步降低频率，或者一直更改为 DC（恒定相位，频率 = 0）。

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

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

对于此方法的性能，您可能希望 sincosf 同时计算两个值，如果您的编译器没有优化您的调用。（希望 auto-vectorize 通过调用 SIMD sincosf。）

Answer 2

您观察到的是浮点精度问题。浮点数的精度随幅度降低。

幸运的是，对于相位或类似的 angular 值，有一个简单的解决方法，将值环绕在 2π 左右，这样幅度永远不会太大。

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

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

正弦函数频率突然跳变

Sudden jump in sin function frequency

c

floating-point

trigonometry

signal-processing

simd