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

Posted 2023-02-14

【中文标题】使用帧之间的相位变化从 FFT Bins 中提取精确频率

【英文标题】：Extracting precise frequencies from FFT Bins using phase change between frames

【发布时间】：2011-06-05 17:09:03

【问题描述】：

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

虽然很棒，但它非常艰难和沉重。这种材料真的让我很紧张。

我从 Stefan 的代码模块中提取了数学，该模块计算给定 bin 的确切频率。但我不明白最后的计算。有人能解释一下最后的数学结构吗？

在深入研究代码之前，让我先设置一下场景：

假设我们设置 fftFrameSize = 1024，所以我们要处理 512+1 个 bin

例如，Bin[1] 的理想频率适合帧中的单个波。在 40KHz 的采样率下，tOneFrame = 1024/40K 秒 = 1/40s，因此 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;
    

我只是看不清楚，尽管它似乎在盯着脸看。有人可以从头开始一步一步地解释这个过程吗？

【问题讨论】：

如何获得BIN * bins 它代表什么？

【参考方案1】：

基本原理很简单。如果给定的组件与 bin 频率完全匹配，则其相位不会从一个 FT 变为下一个 FT。然而，如果频率不完全对应于 bin 频率，则在连续 FT 之间会有相位变化。频率增量只是：

delta_freq = delta_phase / delta_time

然后对组件频率的精确估计将是：

freq_est = bin_freq + delta_freq

【讨论】：

抱歉我很笨，但我还是不明白为什么这是真的。使用这个数学我仍然觉得很不靠谱。

如果 2 个 FFT 的偏移量与正弦波的一个周期不同，那么即使正弦波频率以 bin 为中心，也会发生相位变化。

知道频率的一个定义是相位变化率，即f = dϕ/dt。

我会冒险有人嫉妒你的 l33tDSPsk1llz :p 好吧，不是我。我非常感谢您和 HotPaw 提供了全新的视角。现在我真的可以理解这一点了——终于！！！

@Ohmu：很高兴听到你取得了进展——如果你打算做更多这类事情，我建议阅读一本好的介绍性 DSP 书——Richard Lyons 的书，了解数字信号处理，非常好，比大多数人实用得多。

【参考方案2】：

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

如果您及时查看（固定频率和固定幅度）正弦波上的单个点，相位将随时间提前与频率成正比的量。或者你也可以反过来：如果你测量一个正弦曲线的相位在任何单位时间内变化了多少，你就可以计算出那个正弦曲线的频率。

相位声码器使用两个 FFT 参考两个 FFT 窗口估计相位，两个 FFT 的偏移量是 2 个相位测量值之间的时间距离。从那时起，您就有了对该 FFT bin 的频率估计（FFT bin 大致是一个滤波器，用于隔离正弦分量或其他适合该 bin 的足够窄带信号）。

要使此方法起作用，使用的 FFT bin 附近的频谱必须相当稳定，例如频率不变等。这是相位声码器需要的假设。

【参考方案3】：

恰好落在 bin 频率上的信号频率将 bin 相位提前 2π 的整数倍。由于 FFT 的周期性，对应于 bin 频率的 bin 相位是 2π 的倍数，因此在这种情况下没有相位变化。你提到的文章也解释了这一点。

【讨论】：

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

我已经回答了我自己的问题。但是，如果我对我的答案给予赏金，它将丢失。由于他很棒的开源项目（Performous），我打算把它给 Tronic，但他有很多积分！所以......享受;）

【参考方案4】：

我自己为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;

通过假设 bin 中的信号具有 bin 中心频率，然后计算预期相移来进行校正。从实际班次中减去这个预期班次，留下误差。取余数（模 2 pi）（-pi 到 pi 范围）并使用 bin 中心 + 校正计算最终频率。

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

请注意，许多相邻的 bin 通常最终被校正为相同的频率，因为 delta 校正可以高达 0.5 * FFT_N / FFT_STEP bin所需的处理能力以及由于不准确而导致的不精确性）。

我希望这会有所帮助:)

【讨论】：

我现在有一些“论文风格”的基本原理可供参考。但我不够聪明，无法从这些解释中自己制定数学。我经过一些解释，逐行生成数学。数学证明。

也许这会有所帮助？ sengpielaudio.com/calculator-timedelayphase.htm（时间延迟以毫秒为单位，但我想您可以将 256 个样本转换为适当的时间）

【参考方案5】：

我终于想通了；真的我不得不从头开始。我知道会有一些简单的方法来推导它，我的（通常）错误是试图遵循别人的逻辑而不是使用我自己的常识。

这个谜题需要两把钥匙来解锁。

第一个关键是了解过采样如何在 bin 相位上引入旋转。

第二个键来自此处的图 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;

【讨论】：

编辑：查看我在math.stackexchange.com/questions/9416/… 的回答以了解垃圾箱轮换

【参考方案6】：

也许这会有所帮助。将 FFT 箱视为指定小时钟或转子，每个都以箱的频率旋转。对于稳定的信号，转子的（理论）下一个位置可以使用您没有得到的位中的数学来预测。 针对这个“应该”（理想）位置，您可以计算几个有用的东西：(1) 与相邻帧的 bin 中的相位差，phase vocoder 使用它来更好地bin 频率的估计，或者 (2) 更一般地说 相位偏差，这是音频中音符开始或其他事件的积极指标。