使用帧之间的相变从FFT箱中提取精确的频率

https://stackoverflow.com/questions/4633203

08-10-2019
|

题

我一直在研究这篇精彩的文章： http://blogs.zynaptiq.com/bernsee/pitch-shifting-using-the-ft/

虽然很棒，但它非常困难和沉重。这种材料真的在伸展我。

我从Stefan的代码模块中提取了数学，该模块计算给定垃圾箱的确切频率。但是我不明白最后一个计算。有人可以向我解释最后的数学结构吗？

在挖掘代码之前，让我设置场景：

假设我们设置了fftFrameize = 1024，所以我们正在处理512+1个垃圾箱
例如，bin [1]的理想频率适合框架中的单个波。以40kHz的样本速率，ToneFrame = 1024/40K秒= 1/40s，因此bin [1]理想地将收集40Hz信号。
设置OSAMP（超级样本）= 4，我们以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;

（编辑:)现在我没有得到的一点：

        // 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;
    }
}

我只是看不到它，即使它似乎盯着脸。有人可以逐步从头开始解释这个过程吗？

解决方案 4

最后，我已经弄清楚了。真的，我不得不从头开始得出它。我知道会有一些简单的方法来得出它，我（通常）的错误是尝试遵循他人的逻辑，而不是使用我自己的常识。

这个难题 两个钥匙 解锁。

第一个键 是要了解过度采样如何在箱相中引入旋转。
第二个键 来自图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;
}

其他提示

基本原理非常简单。如果给定的组件完全匹配垃圾箱频率，则其相位不会从一个ft变为另一个ft。但是，如果频率与垃圾箱频率不完全相对应，则连续的FT之间将发生相变。频率delta只是：

delta_freq = delta_phase / delta_time

然后，组件频率的精制估计将为：

freq_est = bin_freq + delta_freq

我已经实施了此算法性能我。当您一次偏移另一个FFT时，您希望阶段根据偏移发生变化，即两个FFT取256个样本分开的相位差应为信号中所有频率的256个样本（这假设信号本身就是信号本身是稳定的，这是256个样本（例如256个样本）的好假设。

现在，您从FFT获得的实际相值不在样品中，而是相位角度，因此根据频率的不同。在以下代码中，phaseStep值是每个bin所需的转换因子，即与bin x相对应的频率，相移将为x * phaseStep。对于bin中心频率，x将是一个整数（箱号），但对于实际检测到的频率，可能是任何实际数字。

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

校正是通过假设垃圾箱中的信号具有垃圾箱中心频率然后计算预期相移的方法来起作用的。这种预期的偏移是从实际偏移中提取的，留下了错误。取出剩余的（Modulo 2 Pi）（-PI至PI范围），并使用bin Center +校正计算最终频率。

// 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 bin，以便您使用的较小的fft_step，您使用的较小的FFT_STEP将可以校正（但这会增加处理能力由于不准确而需要不精确）。

我希望这有帮助：）

这是相位声码器方法使用的频率估计技术。

如果您在时间上查看（固定频率和固定幅度）正弦波的单个点，则相位将随时间推进，与频率成正比。或者，您可以进行相反：如果测量正弦曲线的相位在任何时间单位上变化多少，则可以计算正弦的频率。

相位声码编码器使用两个FFT来估计两个FFT窗口，而两个FFT的偏移是及时两相测量之间的距离。从那里开始，您可以对FFT箱进行频率估计（FFT箱大致是隔离正弦分量或其他足够窄带信号的滤波器，该信号适合该箱内）。

为了使这种方法起作用，使用的FFT箱附近的光谱必须相当固定，例如不变频率等。这是相位声码器所需的假设。

也许这会有所帮助。将FFT箱视为指定的小钟或转子，每个时钟都以垃圾箱的频率旋转。对于稳定的信号，可以使用数学在您无法获得的位置预测转子的下一个位置。在此“应该是”（理想的）位置上，您可以计算几个有用的事情：（1）相邻框架箱中的相位的差异 相位声码器 更好地估计箱频率，或（2）更普遍 相位偏差, ，这是音符发作或音频中其他事件的积极指标。

通过2π的整数倍数，准确地落在bin频率提前箱相的信号频率。由于与bin频率相对应的箱相是由于FFT的周期性，因此在这种情况下没有相变。您提到的文章也解释了这一点。

许可以下： CC-BY-SA 和归因

不隶属于 StackOverflow