Извлечение точных частот из ячеек БПФ с использованием изменения фазы между кадрами

Основной принцип очень прост. Если данный компонент точно соответствует частоте бина, то его фаза не будет меняться от одного FT к другому. Однако, если частота не соответствует точно частоте бина, тогда будет фазовый сдвиг между последовательными FT. Дельта частоты равна:

delta_freq = delta_phase / delta_time

и тогда уточненная оценка частоты компонента будет:

freq_est = bin_freq + delta_freq

Я сам реализовал этот алгоритм для Performous. Когда вы берете другое БПФ со смещением по времени, вы ожидаете, что фаза изменится в соответствии со смещением, т. е. два БПФ, взятые с интервалом 256 отсчетов, должны иметь разность фаз 256 отсчетов для всех частот, присутствующих в сигнале (это предполагает, что сами сигналы стабильны, что является хорошим предположением для коротких периодов, таких как 256 выборок).

Теперь фактические значения фазы, которые вы получаете из БПФ, представлены не в выборках, а в фазовом угле, поэтому они будут разными в зависимости от частоты. В следующем коде значение PhaseStep является коэффициентом преобразования, необходимым для каждого бина, т. е. для частоты, соответствующей бину x, фазовый сдвиг будет равен x * PhaseStep. Для центральных частот бина x будет целым числом (номер бина), но для фактически обнаруженных частот это может быть любое действительное число.

const double freqPerBin = SAMPLE_RATE / FFT_N;
const double phaseStep = 2.0 * M_PI * FFT_STEP / FFT_N;

Коррекция работает, предполагая, что сигнал в бине имеет центральную частоту бина, а затем вычисляя для этого ожидаемый фазовый сдвиг. Это ожидаемое смещение вычитается из фактического сдвига, оставляя ошибку. Берется остаток (по модулю 2 пи) (диапазон от -пи до пи), и окончательная частота рассчитывается с центром бина + поправкой.

// process phase difference
double delta = phase - m_fftLastPhase[k];
m_fftLastPhase[k] = phase;
delta -= k * phaseStep;  // subtract expected phase difference
delta = remainder(delta, 2.0 * M_PI);  // map delta phase into +/- M_PI interval
delta /= phaseStep;  // calculate diff from bin center frequency
double freq = (k + delta) * freqPerBin;  // calculate the true frequency

Обратите внимание, что многие соседние бины часто корректируются до одной и той же частоты, потому что дельта-коррекция может составлять до 0,5 * бинов FFT_N / FFT_STEP в любом случае, поэтому чем меньший FFT_STEP вы используете, тем дальше будут возможны коррекции (но это увеличивает вычислительную мощность). необходимости, а также неточности из-за неточностей).

Надеюсь, это поможет :)

Это метод оценки частоты, используемый в методах фазового вокодера.

Если вы посмотрите на одну точку синусоидальной волны (фиксированной частоты и фиксированной амплитуды) во времени, фаза будет увеличиваться со временем на величину, пропорциональную частоте. Или вы можете сделать обратное: если вы измерите, насколько фаза синусоиды изменяется в любую единицу времени, вы можете вычислить частоту этой синусоиды.

Фазовый вокодер использует два БПФ для оценки фазы относительно двух окон БПФ, а смещение двух БПФ представляет собой расстояние между двумя измерениями фазы во времени. Отсюда у вас есть оценка частоты для этого бина БПФ (бинар БПФ, грубо говоря, является фильтром для выделения синусоидальной составляющей или другого достаточно узкополосного сигнала, который помещается в этот бин).

Чтобы этот метод работал, спектр вблизи используемого бина БПФ должен быть достаточно стационарным, т.е. не меняется частота и т. д. Это предположение требуется для фазового вокодера.

Наконец я понял это; действительно мне пришлось вывести его с нуля. Я знал, что будет какой-то простой способ вывести его, моя (обычная) ошибка заключалась в попытке следовать логике других людей, а не использовать собственный здравый смысл.

Эта головоломка требует два ключа, чтобы разблокировать ее.

• Первый ключ – понять, как избыточная выборка приводит к чередованию фаз бинов.

• Второй ключ взят из диаграмм 3.3 и 3.4 здесь: http://www.dspdimension.com/admin/pitch-shifting-using-the-ft/

...

for (int k = 0; k <= fftFrameSize/2; k++) 
{
    // compute magnitude and phase 
    bins[k].mag = 2.*sqrt(fftBins[k].real*fftBins[k].real + fftBins[k].imag*fftBins[k].imag);
    bins[k].phase = atan2(fftBins[k].imag, fftBins[k].real);

    // Compute phase difference Δϕ fo bin[k]
    double deltaPhase;
    {
        double measuredPhaseDiff = bins[k].phase - gLastPhase[k];
        gLastPhase[k] = bins[k].phase;

        // Subtract expected phase difference <-- FIRST KEY
        // Think of a single wave in a 1024 float frame, with osamp = 4
        //   if the first sample catches it at phase = 0, the next will 
        //   catch it at pi/2 ie 1/4 * 2pi
        double binPhaseExpectedDiscrepancy = M_TWOPI * (double)k / (double)osamp;
        deltaPhase = measuredPhaseDiff - binPhaseExpectedDiscrepancy;

        // Wrap delta phase into [-Pi, Pi) interval 
        deltaPhase -= M_TWOPI * floor(deltaPhase / M_TWOPI + .5);
    }

    // say sampleRate = 40K samps/sec, fftFrameSize = 1024 samps in FFT giving bin[0] thru bin[512]
    // then bin[1] holds one whole wave in the frame, ie 44 waves in 1s ie 44Hz ie sampleRate / fftFrameSize
    double bin0Freq = (double)sampleRate / (double)fftFrameSize;
    bins[k].idealFreq = (double)k * bin0Freq;

    // Consider Δϕ for bin[k] between hops.
    // write as 2π / m.
    // so after m hops, Δϕ = 2π, ie 1 extra cycle has occurred   <-- SECOND KEY
    double m = M_TWOPI / deltaPhase;

    // so, m hops should have bin[k].idealFreq * t_mHops cycles.  plus this extra 1.
    // 
    // bin[k].idealFreq * t_mHops + 1 cycles in t_mHops seconds 
    //   => bins[k].actualFreq = bin[k].idealFreq + 1 / t_mHops
    double tFrame = fftFrameSize / sampleRate;
    double tHop = tFrame / osamp;
    double t_mHops = m * tHop;

    bins[k].freq = bins[k].idealFreq + 1. / t_mHops;
}

Может быть, это поможет. Думайте о ячейках БПФ как о маленьких часах или роторах, каждый из которых вращается с частотой ячейки. Для стабильного сигнала (теоретическое) следующее положение ротора можно предсказать, используя математику в бите, который вы не получаете. По отношению к этому «должно быть» (идеальному) положению вы можете вычислить несколько полезных вещей: (1) разницу с фазой в ячейке соседнего кадра, которая используется фазовым вокодером для лучшего оценка частоты бина или (2) в более общем случае отклонение фазы, которое является положительным индикатором начала ноты или какого-либо другого события в аудио.

Частоты сигналов, попадающие точно на частоту элемента разрешения, опережают фазу элемента разрешения на целое число, кратное 2π. Поскольку фазы бинов, соответствующие частотам бинов, кратны 2π из-за периодического характера БПФ, в этом случае фазового изменения нет. Упомянутая вами статья также объясняет это.