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

Я просматривал эту фантастическую статью: http://blogs.zynaptiq.com/bernsee/pitch-shifting-using-the-ft/

Хотя это фантастика, это очень тяжело и тяжело. Этот материал меня действительно напрягает.

Я извлек математические данные из модуля кода Стефана, который вычисляет точную частоту для заданного интервала. Но я не Я понимаю последний расчет. Может кто-нибудь объяснить мне математическую конструкцию в конце?

Прежде чем углубиться в код, позвольте мне установить сцену:

• Допустим, мы установили fftFrameSize = 1024, поэтому мы имеем дело с 512 + 1 ячейками

• Например, идеальная частота Bin [1] соответствует одной волне в кадре. При частоте дискретизации 40 кГц tOneFrame = 1024 / 40K секунд = 1 / 40s, поэтому Bin [1] в идеале будет собирать сигнал 40 Гц.

• Установив osamp (overSample) = 4, мы продвигаемся по входному сигналу пошагово of 256. Итак, первый анализ исследует байты с нуля по 1023, затем с 256 по 1279 и т. д. Обратите внимание, что каждое число с плавающей запятой обрабатывается 4 раза.

...

void calcBins( 
              long fftFrameSize, 
              long osamp, 
              float sampleRate, 
              float * floats, 
              BIN * bins
              )
{
    /* initialize our static arrays */
    static float gFFTworksp[2*MAX_FRAME_LENGTH];
    static float gLastPhase[MAX_FRAME_LENGTH/2+1];

    static long gInit = 0;
    if (! gInit) 
    {
        memset(gFFTworksp, 0, 2*MAX_FRAME_LENGTH*sizeof(float));
        memset(gLastPhase, 0, (MAX_FRAME_LENGTH/2+1)*sizeof(float));
        gInit = 1;
    }

    /* do windowing and re,im interleave */
    for (long k = 0; k < fftFrameSize; k++) 
    {
        double window = -.5*cos(2.*M_PI*(double)k/(double)fftFrameSize)+.5;
        gFFTworksp[2*k] = floats[k] * window;
        printf("sinValue: %f", gFFTworksp[2*k]);
        gFFTworksp[2*k+1] = 0.;
    }

    /* do transform */
    smbFft(gFFTworksp, fftFrameSize, -1);

    printf("\n");

    /* this is the analysis step */
    for (long k = 0; k <= fftFrameSize/2; k++) 
    {
        /* de-interlace FFT buffer */
        double real = gFFTworksp[2*k];
        double imag = gFFTworksp[2*k+1];

        /* compute magnitude and phase */
        double magn = 2.*sqrt(real*real + imag*imag);
        double phase = atan2(imag,real);

        /* compute phase difference */
        double phaseDiff = phase - gLastPhase[k];
        gLastPhase[k] = phase;

        /* subtract expected phase difference */
        double binPhaseOffset = M_TWOPI * (double)k / (double)osamp;
        double deltaPhase = phaseDiff - binPhaseOffset;

        /* map delta phase into [-Pi, Pi) interval */
        // better, but obfuscatory...
        //    deltaPhase -= M_TWOPI * floor(deltaPhase / M_TWOPI + .5);

        while (deltaPhase >= M_PI)
            deltaPhase -= M_TWOPI;
        while (deltaPhase < -M_PI)
            deltaPhase += M_TWOPI;

(EDIT :) Теперь бит, который я не получаю:

        // Get deviation from bin frequency from the +/- Pi interval 
        // Compute the k-th partials' true frequency    

        // Start with bin's ideal frequency
        double bin0Freq = (double)sampleRate / (double)fftFrameSize;
        bins[k].idealFreq = (double)k * bin0Freq;

        // Add deltaFreq
        double sampleTime = 1. / (double)sampleRate;
        double samplesInStep = (double)fftFrameSize / (double)osamp;
        double stepTime = sampleTime * samplesInStep;
        double deltaTime = stepTime;        

        // Definition of frequency is rate of change of phase, i.e. f = dϕ/dt
        // double deltaPhaseUnit = deltaPhase / M_TWOPI; // range [-.5, .5)
        double freqAdjust = (1. / M_TWOPI) * deltaPhase / deltaTime; 

        // Actual freq <-- WHY ???
        bins[k].freq = bins[k].idealFreq + freqAdjust;
    }
}

Я просто не могу этого ясно видеть, хотя кажется, что он смотрит прямо в лицо. Может кто-нибудь объяснить этот процесс с нуля, шаг за шагом?