Я пишу источник сигнала для цифровой обработки сигналов и обнаружил странное поведение. Вот код:
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);
Это создаст файл со сложным сигналом, смещенным на -3355 Гц от центра. Поэтому, когда я запускаю БПФ над этим файлом, я ожидаю сигнал на частоте -3355 Гц. С некоторым незначительным частотным распределением вокруг этого числа из-за квантования. Но на самом деле я получаю следующее:
Есть довольно значительный скачок частоты (~2000 Гц) в районе 50 секунд.
Кто-нибудь знает, почему?
Мои эксперименты показали:
phase
помогаетphase
нет переполнения@Sedenion, это воображаемый юнит, импортированный из complex.h
Тема генерации бесконечного синуса обсуждалась на stackoverflow.com/q/69729326/4989451. Это не ответ на ваш вопрос, но вы можете найти его полезным
MS справочная страница для sin()
говорит: «Если Икс больше или равно 263 или меньше или равно -263, происходит потеря значимости в результате».
Как только phase
становится достаточно большим, приращения к нему квантуются путем округления с плавающей запятой, может быть, достаточно, чтобы округление до ближайшего с тай-брейком округления до четной мантиссы не уравновешивалось, а вместо этого означало, что вы последовательно продвигаете фазу меньше. См. статья в Википедии о плавающей запятой одинарной точности IEEE
Аналогия с основанием 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
и посмотреть, получите ли вы еще один шаг вниз по частоте или изменение до постоянного тока (постоянная фаза, частота = 0).
Или с помощью float
разверните вручную, чтобы у вас было два счетчика, смещенных на phase_increment
, и вы увеличили каждый на 2*phase_increment
. Таким образом, относительная величина phase
по сравнению с тем, что добавляется, не становится такой экстремальной. Однако они по-прежнему будут иметь ту же абсолютную величину, большую по сравнению с исходной phase_increment
, поэтому некоторое округление все же будет.
Я не знаю, есть ли лучшие стратегии для создания сложных синусоидальных волн, я просто пытаюсь объяснить поведение, которое вы видите.
Для производительности этого метода вы, вероятно, захотите sincosf
вычислять оба значения одновременно, если ваш компилятор не оптимизирует ваши вызовы для этого. (И, надеюсь, автоматизировать векторизацию, вызвав SIMD sincosf
.)
Понятно. Благодарю вас! Кстати: -O2 создаст sincosf для этого кода.
То, что вы наблюдаете, - это проблемы с точностью с плавающей запятой. Точность поплавков ухудшается с величиной.
К счастью, для фаз или подобных угловых значений есть простой обходной путь: оберните значение вокруг 2π, чтобы величина никогда не была слишком большой.
Добавьте следующий код после phase += phase_increment;
и посмотрите, поможет ли это.
if ( phase < (float)( 2.0 * M_PI ) )
continue;
phase -= (float)( 2.0 * M_PI );
Да, редукция работает хорошо: ``` if (phase < -(float)( 2.0 * M_PI )) { Phase += (float)( 2.0 * M_PI ); } if ( фаза > (с плавающей запятой) ( 2.0 * M_PI )) { фаза -= (с плавающей запятой) ( 2.0 * M_PI ); }```
в
output[j] = cosf(phase) + sinf(phase) * I;
является ли эта большаяI
другой переменной, чем маленькаяi
, и если да, то что это за переменнаяI
?