突然跳跃罪名频率
我正在编写数字信号处理的信号源，并发现了奇怪的行为。这是代码：
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）。
有人知道为什么吗？
的实验表明：
- 使用double for
epase帮助 - 将参数从-2pi减少到2PI，帮助
- Apple M1和Raspberry Pi也有同样的跳跃
- 我
阶段no Overflow参数 - 有趣阅读 https://randomascii.wordpress.com/2014/10/10/10/10/10/intel-underestimates-erestimates-error-bounds-bounds-bounds-by-1-3-quintillion/ ，但不太可能引起突然引起的跳跃
I'm writing signal source for digital signal processing and found strange behaviour. Here is the code:
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);
It will create a file with complex signal shifted by -3355hz from the centre. So when I run FFT over this file I expect a signal at -3355hz. With some minor frequency distribution around that number due to quantisation. But in reality I'm getting the following:
There is a quite significant frequency jump (~2000hz) around 50 seconds.
Anyone knows why?
My experiments showed:
- using double for the
phasehelps - doing argument reduction from -2pi to 2pi helps
- Apple M1 and raspberry pi have the same jump
- there seems to be no overflow in the
phaseargument - Fun reading https://randomascii.wordpress.com/2014/10/09/intel-underestimates-error-bounds-by-1-3-quintillion/ but unlikely causing sudden jump of a frequency
如果你对这篇内容有疑问，欢迎到本站社区发帖提问 参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。
绑定邮箱获取回复消息
由于您还没有绑定你的真实邮箱，如果其他用户或者作者回复了您的评论，将不能在第一时间通知您！
发布评论
评论（2）
一旦
阶段变得足够大，可以通过浮点舍入来量化它，也许足以使圆形到最接近的圆形圆形到曼蒂萨（Mantissa）都无法平衡，而是意味着您始终如一地较少阶段。参见href = “ https://en.wikipedia.org/wiki/single-precision_floation-point-point_format _format#precision_limitations_on_decimal_values_(Betefiend_1_1_1_and_and_and_and_and_and_and_and_and_and_and_and_and_and_and_and_16777216）“ -10与值限制为4个重要数字的类比是
101.2 + 0.111仅到101.3，而不是101.311。 （计算机使用二进制浮点部，因此实际上在101.1不在101.1的情况下，实际上是可以代表的。）我怀疑这就是发生的事情。您不会溢出到
+Inf，但是给定float的相对精度有限，将少量数字添加到一个巨大的数字最终将完全不会移动：最接近代表的浮点与阶段 + aphe_increment将是原始apeas。测试此假设的一种快速方法是使用
double Phase（不进行任何其他更改），看看是否会移动您获得频率步骤的点。 （要变得更大，可能已经超过了您生成的内容）。（哦，只是注意到您已经做到了，并且确实有帮助。如果您的意思是完全消除了问题，那可能就是这样
。 COSF ，尽管更改它们是另一个要尝试的事情：正如@Weather Vane指出的那样，一些
sin实现在范围减少方面失去了显着精确度。但是您希望这是一个嘈杂的信号，不是频率的一致变化（这将是该阶段的比例因素。）您还可以让它使用
float，看看是否是否频率下一步降低，或者一直更改为DC（恒定阶段，频率= 0）。或使用
float，手动展开，因此您的两个计数器由sthepe_increment偏移，然后通过2*stape_increment来递增每个计数器。因此，阶段与所添加的内容的相对大小并没有那么极端。但是，相对于原始的apeas_increment，它们仍然可以达到相同的绝对幅度，因此仍然会有一些舍入。我不知道是否有更好的策略来产生复杂的正弦波，我只是想解释您看到的行为。
为了执行此方法，您可能希望
sincosf同时计算两个值，如果您的编译器未对此进行优化。 （并希望通过调用SIMDsincosf来自动矢量化。）Once
phasegets large enough, increments to it are quantized by floating-point rounding, maybe enough that rounding to nearest with tie-break of rounding to even mantissa doesn't balance out, and instead means you consistently advance the phase less. See wikipedia's article on single-precision IEEE floatA base-10 analogy with values limited to 4 significant figures would be
101.2 + 0.111only going up to101.3, not101.311. (Computers use binary floats, so actually stuff like 101.125 is exactly representable where 101.1 isn't.)I suspect that's what's happening. You won't overflow to
+Inf, but given the limited relative precision of afloat, adding a small number to a huge number will eventually not move it at all: the closest representable float tophase + phase_incrementwill be the originalphase.A quick way to test this hypothesis would be to use
double phase(without making any other changes), and see if that moves the point at which you get frequency steps. (To be vastly greater, probably past the end of what you're generating).(Oh, just noticed you already did that, and it did help. If you mean it completely removed the problem, then that's probably it.)
You can still keep the single-precision
sinfandcosf, although changing them would be the other thing to try: as @Weather Vane points out, somesinimplementations lose significant precision in range-reduction. But you'd expect that to be a noisier signal, not a consistent change in frequency (which would take a scale factor for the phase.)You could also let it run longer with
floatand see if you get another step down in frequency, or a change all the way to DC (constant phase, frequency = 0).Or with
float, manually unroll so you have two counters offset byphase_incrementand you increment each one by2*phase_increment. So the relative magnitude ofphasevs. what's being added doesn't get so extreme. They'll still get to the same absolute magnitude, though, large relative to the originalphase_increment, so there will still be some rounding.I don't know if there are better strategies for generating complex sine waves, I'm just trying to explain the behaviour you're seeing.
For performance of this method, you probably want
sincosfto compute both values at the same time, if your compiler doesn't optimize your calls into that. (And hopefully auto-vectorize by calling a SIMDsincosf.)您观察到的是浮点精度问题。浮子的精度大小降解。
幸运的是，对于阶段或类似的角度值，有一个简单的解决方法，将值包裹在2π上，因此幅度永远不会太大。
在您的
apeas += sthepe_increment;之后添加以下代码，看看是否有帮助。What you observing is floating point precision issues. Precision of floats degrades with magnitude.
Fortunately, for phases or similar angular values there’s an easy workaround, wrap the value around 2π so the magnitude is never too large.
Add the following code after your
phase += phase_increment;and see if that helps.