从FFT Bin中提取精确的频率，利用帧之间的相位变化。

Question

从FFT Bin中提取精确的频率，利用帧之间的相位变化。

cmathsignal-processingfftphase

27

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

虽然这篇文章很棒，但是它非常难以理解和消化。这些内容真的让我感到困扰。

我从Stefan的代码模块中提取了计算给定bin的准确频率的数学公式。但是我不理解最后的计算。有人能向我解释一下最后的数学构造吗？

在深入研究代码之前，让我先介绍一下背景：

假设我们设置fftFrameSize = 1024，因此我们正在处理512 + 1个bin
举个例子，Bin [1] 的理想频率在帧中完整地容纳了一个波形。以40KHz的采样率为例，tOneFrame = 1024 / 40K 秒= 1/40秒，因此Bin [1] 理论上应该采集40Hz的信号。
设置osamp（overSample）= 4，我们以256的步长沿着输入信号前进。因此，第一次分析检查0到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;
    }
}

虽然似乎它就在面前盯着我，但我仍然无法清晰地看到它。请问有人可以从头开始，一步一步地解释这个过程吗？

- P i

如何获得“BIN * bins”，它代表什么？ - Mord Fustang

6个回答

11

我已经在 Performous 中实现了这个算法。当你在一个时间偏移处再次进行FFT时，你期望相位会根据偏移改变，即两个相距256个样本的FFT应该对所有频率的信号都有256个样本的相位差异（假设信号本身是稳定的，在短周期内如256个样本是一个很好的假设）。

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

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

这个修正方法的原理是假设某个频率区间的信号具有该区间中心频率，然后计算该信号预期的相位偏移量。将预期偏移量从实际偏移量中减去，得到误差值。然后对误差值取模（范围为-π到π），得到余数，最后通过该区间中心频率加上修正值来计算最终频率。

// 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 越小，修正距离可能越远（但这也会增加所需的处理能力和由于不准确性而产生的不精确度）。

希望这有所帮助 :)

- Tronic

我现在有几个“论文风格”的理由要看。但是我不够聪明，无法从这些解释中自己构思出数学公式。我需要一些能够逐行生成数学公式的解释。一份数学证明。 - P i

也许这个可以帮助你？http://www.sengpielaudio.com/calculator-timedelayphase.htm（那里的时间延迟以毫秒为单位，但我认为你可以将256个样本转换成正确的时间量） - Tronic

7

最终我弄清楚了，真的是从零开始推导出来的。我知道一定有某种简单的方法可以推导出来，但我的错误通常是试图跟随他人的逻辑而不是运用自己的常识。

这个谜题需要两个关键来解锁它。

第一个关键是理解过采样如何在频率相位上引入旋转。
第二个关键来自于这里的图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;
}

- P i

请参考我的回答：http://math.stackexchange.com/questions/9416/extracting-exact-frequencies-from-fft-output/783273#783273，以了解二进制旋转。 - P i

7

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

如果你观察一个（固定频率和振幅）正弦波在时间上的单个点，相位会随着时间按比例增加。或者你可以反过来做：如果你测量正弦波的相位在任意单位时间内发生了多少变化，你就可以计算出该正弦波的频率。

相位变换器使用两个FFT来估计相位，并参考两个FFT窗口，两个FFT之间的偏移量是2个相位测量之间的时间距离。从此，你就可以得到该FFT bin的频率估计值（FFT bin大致上是用于隔离适合于该bin的正弦波分量或其他足够窄带信号的滤波器）。

为了使这种方法有效，使用的FFT bin附近的频谱必须相当稳定，例如不改变频率等。这是相位变换器所需的假设。

- hotpaw2

2

也许这可以帮助你。将FFT频率范围看作指定小时钟或转子的位置，每个转子以频率为速度旋转。对于稳定的信号，可以使用未知部分中的数学来预测转子的下一个（理论）位置。根据这个“应该是”的（理想）位置，您可以计算出几个有用的东西：（1）与相邻帧的一个频率范围内的相位差异，由 相位变换器 用于更好地估算频率范围，或者（2）更一般的相位偏差，这是音频中音符起始或其他事件的正指标。

- Bill

1

落在频率格点上的信号频率会使得格点相位提前2π的整数倍。由于FFT的周期性，对应于频率格点的相位也是2π的整数倍，因此在这种情况下不会发生相位变化。你提到的文章也解释了这一点。

- tahome

如果FFT步长与FFT大小相同，那就是正确的。然而，在这里，步骤会变得更小（osamp因子），然后即使对于中心频率，相位也不再保持不变。例如，考虑仅为一个样本的FFT步骤。对于较低的频率，基本上不会有相移，而对于非常高的频率，可能会有高达PI相位差。 - Tronic

我已经回答了自己的问题。但是如果我把赏金给我的答案，它就会丢失。我本来想把它给Tronic，因为他有一个很棒的开源项目（Performous），但他已经有很多积分了！所以...尽情享受吧 ;) - P i

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Paul R · Accepted Answer

基本原理非常简单。如果给定的组件恰好匹配一个 bin 频率，则其相位将不会从一个 FT 改变到下一个。然而，如果频率与 bin 频率不完全对应，则在连续的 FT 之间将会有相位变化。频率差值就是：

基本原理非常简单。如果给定的组件恰好匹配一个 bin 频率，则其相位将不会从一个傅里叶变换（FT）到下一个发生改变。但是，如果频率与 bin 频率不完全对应，则在连续的 FT 之间将会发生相位变化。频率差值就是：

delta_freq = delta_phase / delta_time

然后组成部分的精细估计频率将为：

freq_est = bin_freq + delta_freq