具有纯正弦音的 FFT 更精确的频率

Posted 2023-02-25

【中文标题】具有纯正弦音的 FFT 更精确的频率

【英文标题】：More precise frequency from FFT with pure sine tones

【发布时间】：2014-01-26 15:25:24

【问题描述】：

我目前正在使用此处的 FFT 代码： https://github.com/syedhali/EZAudio/tree/master/EZAudioExamples/ios/EZAudioFFTExample

以下是 2 个相关方法的代码：

-(void)createFFTWithBufferSize:(float)bufferSize withAudioData:(float*)data 

  // Setup the length
  _log2n = log2f(bufferSize);

  // Calculate the weights array. This is a one-off operation.
  _FFTSetup = vDSP_create_fftsetup(_log2n, FFT_RADIX2);

  // For an FFT, numSamples must be a power of 2, i.e. is always even
  int nOver2 = bufferSize/2;

  // Populate *window with the values for a hamming window function
  float *window = (float *)malloc(sizeof(float)*bufferSize);
  vDSP_hamm_window(window, bufferSize, 0);
  // Window the samples
  vDSP_vmul(data, 1, window, 1, data, 1, bufferSize);
  free(window);

  // Define complex buffer
 _A.realp = (float *) malloc(nOver2*sizeof(float));
 _A.imagp = (float *) malloc(nOver2*sizeof(float));



-(void)updateFFTWithBufferSize:(float)bufferSize withAudioData:(float*)data 

  // For an FFT, numSamples must be a power of 2, i.e. is always even
  int nOver2 = bufferSize/2;

  // Pack samples:
  // C(re) -> A[n], C(im) -> A[n+1]
  vDSP_ctoz((COMPLEX*)data, 2, &_A, 1, nOver2);

  // Perform a forward FFT using fftSetup and A
  // Results are returned in A
  vDSP_fft_zrip(_FFTSetup, &_A, 1, _log2n, FFT_FORWARD);

  // Convert COMPLEX_SPLIT A result to magnitudes
  float amp[nOver2];
  float maxMag = 0;

  for(int i=0; i<nOver2; i++) 
    // Calculate the magnitude
    float mag = _A.realp[i]*_A.realp[i]+_A.imagp[i]*_A.imagp[i];
    maxMag = mag > maxMag ? mag : maxMag;
  
  for(int i=0; i<nOver2; i++) 
    // Calculate the magnitude
    float mag = _A.realp[i]*_A.realp[i]+_A.imagp[i]*_A.imagp[i];
   // Bind the value to be less than 1.0 to fit in the graph
   amp[i] = [EZAudio MAP:mag leftMin:0.0 leftMax:maxMag rightMin:0.0 rightMax:1.0];
 

我已经修改了上面的 updateFFTWithBufferSize 方法，这样我就可以像这样获得以赫兹为单位的频率：

for(int i=0; i<nOver2; i++) 
    // Calculate the magnitude
    float mag = _A.realp[i]*_A.realp[i]+_A.imagp[i]*_A.imagp[i];
    if(maxMag < mag) 
        _i_max = i;
    
    maxMag = mag > maxMag ? mag : maxMag;


float frequency = _i_max / bufferSize * 44100;
NSLog(@"FREQUENCY: %f", frequency); 

我已经用 Audacity 在不同的频率下生成了一些纯正弦音来进行测试。我看到的问题是代码为两个值相对接近的不同正弦音返回相同的频率。

例如： 以 19255Hz 生成的正弦音将从 FFT 显示为 19293.750000Hz。以 19330Hz 产生的正弦音也是如此。计算中一定有问题。

对于如何修改上述代码以获得更精确的纯正弦音 FFT 频率读数的任何帮助，我们将不胜感激。 谢谢！

【问题讨论】：

由于我的声音的采样率为 44.1KHz，看起来频率分辨率将在 20000 / 256 bins 左右，即 78.125。所以这就是为什么我不能得到一个特定的频率，只有那个 78 范围内的一个。还有其他方法可以更具体地了解正弦音吗？

如果您列出相关数字会有所帮助，例如：buffer_size、分辨率等。

我猜我的缓冲区大小目前是 512。我将在哪里/如何增加它？

您的频率分辨率将受到采样频率和缓冲区大小的限制。

你好 Codeman，你能把修改后的增加缓冲区大小的代码发给我吗（我面临着类似的问题，我已经研究了一个月了）

【参考方案1】：

您可以通过将抛物线曲线拟合到峰值幅度 bin 周围的 3 个 FFT bin 幅度，然后找到该抛物线的极值来获得粗略的频率估计。

可以通过将 FFT 窗口的变换用作插值内核并进行逐次逼近来优化对插值点的最大值的估计来创建更好的估计。 （零填充和使用更长的 FFT 将为您提供类似类型的插值估计。）

稳定信号的简单方法是，如果可能的话，使用更长的 FFT 以及跨越更长时间间隔的更多样本。

【讨论】：

如何增加缓冲区大小？那是我可以在代码中手动更改的东西吗？

@codeman：是的，FFT 对它接受的点数没有实际限制。如果你想通过 262.144 分而不是 512，它仍然可以工作。但是，频率分辨率 * 时间分辨率 = 2，因此填充该缓冲区需要 512 倍的时间。

【参考方案2】：

这里有很多问题：

1) 您的频率轴间距为 f_max/N，或大约 80Hz，因此您不会获得比这更好的分辨率。

2) 你的信号非常接近Nyquist frequency（即，20KHz/44.1KHz 几乎是 0.5），当你接近奈奎斯特极限时，如果你想要准确，你需要非常小心结果。 （也就是说，在 20KHz 时，每个完整的振荡周期只记录大约两个数据点。）

3) 由于 20KHz 处于人类听觉的边缘（对于大多数人来说更高），因此许多麦克风并不真正担心它。 Here's iPhone 的测量值。

【参考方案3】：

也许你的采样频率不够高？

【讨论】：

正弦音为 44.1 kHz，据我了解，这适用于我生成的频率，因为其中一半比我生成的声音多。

给定无限时间就足够了。如果您需要在固定时间段内获得更好的频率分辨率，则确实需要提高采样频率。

【参考方案4】：

如果您对输入一无所知，FFT 是一种非常好的获取频谱的方法。如果您知道输入是纯正弦波，您可以做得更好。首先计算 FFT 以大致了解正弦的位置。获取最小值和最大值以估计幅度[或从 FFT 中获取 - 对所有输入求平方，将它们相加，取平方根]，在给定估计的频率和幅度的情况下获取开始和结束时的相位。

通常，您会发现相位不匹配。这是因为最后的相位偏离了 2*Δf * N。f - Δf 是对频率的更好估计。请记住，这种方法对噪声非常敏感。该方法之所以有效，是因为输入是纯正弦波，而噪声就是除此之外的一切。使用这种方法迭代快速炸毁；你甚至会遇到舍入错误（也不是正弦）

另一个类似的技巧是减去估计波。两个正弦之间的差异是两个正弦的乘积，一个增加了频率（在您的情况下为 ±38.5 kHz），另一个减少了频率（Δ_f_，小于 100 Hz）。另见Heterodyne detection