如何使用谐波产品频谱获得基频？

Question

我正在尝试从麦克风输入中获取音高。首先，我通过 FFT 将信号从时域分解到频域。在执行 FFT 之前，我已将汉明窗应用于信号。然后我得到了 FFT 的复杂结果。然后我将结果传递给谐波乘积频谱，在其中对结果进行下采样，然后将下采样的峰值相乘，并给出一个复数值。那我应该怎么做才能得到基频呢？

    public float[] HarmonicProductSpectrum(Complex[] data)
    {
        Complex[] hps2 = Downsample(data, 2);
        Complex[] hps3 = Downsample(data, 3);
        Complex[] hps4 = Downsample(data, 4);
        Complex[] hps5 = Downsample(data, 5);
        float[] array = new float[hps5.Length];

        for (int i = 0; i < array.Length; i++)
        {
            checked
            {
                array[i] = data[i].X * hps2[i].X * hps3[i].X * hps4[i].X * hps5[i].X;
            }
        }
        return array;
    }

    public Complex[] Downsample(Complex[] data, int n)
    {
        Complex[] array = new Complex[Convert.ToInt32(Math.Ceiling(data.Length * 1.0 / n))];
        for (int i = 0; i < array.Length; i++)
        {
            array[i].X = data[i * n].X;
        }
        return array;
    } 

我已经尝试使用，

    magnitude[i] = (float)Math.Sqrt(array[i] * array[i] + (data[i].Y * data[i].Y));  

在 HarmonicProductSpectrum 方法的 for 循环中。然后尝试使用最大的 bin，

        float max_mag = float.MinValue;
        float max_index = -1;

        for (int i = 0; i < array.Length / 2; i++)
            if (magnitude[i] > max_mag)
            {
                max_mag = magnitude[i];
                max_index = i;
            }

然后我尝试使用频率来获取频率，

    var frequency = max_index * 44100 / 1024;

但是对于 A4 音符（440 Hz），我得到了像 1248.926、1205,859、2454.785 这样的垃圾值，而这些值看起来不像 A4 的谐波。

非常感谢您的帮助。

Answer 1

我在 Python 中实现了谐波积谱，以确保您的数据和算法运行良好。

以下是我在将谐波积谱应用于完整数据集时看到的结果，汉明窗，具有 5 个下采样-乘法阶段：

这只是最底层的千赫兹，但在 1 KHz 以上的频谱几乎已经死了。

如果我将长音频片段分成 8192 个样本块（4096 个样本 50% 重叠），并对每个块进行汉明窗口并在其上运行 HPS，这就是 HPS 的矩阵。这是整个数据集上 HPS 频谱的电影。基频似乎相当稳定。

full source code is here——有很多代码可以帮助对数据进行分块并可视化在块上运行的 HPS 的输出，但核心 HPS 功能（从def hps(… 开始）很短。但它有一些技巧。

鉴于您发现峰值的奇怪频率，可能是您在从 0 到 44.1 KHz 的整个频谱上运行？您只想保留“正”频率，即从 0 到 22.05 KHz，并对其应用 HPS 算法（下采样-乘法）。

但是假设您从纯正频率频谱开始，正确地获取它的幅度，看起来您应该得到合理的结果。尝试保存您的HarmonicProductSpectrum 的输出，看看它是否与上述类似。

同样，完整的源代码位于https://gist.github.com/fasiha/957035272009eb1c9eb370936a6af2eb。 （在那里我尝试了另外几个谱估计器，来自 Scipy 的 Welch 方法和我的 Blackman-Tukey 谱估计器端口。我不确定您是否准备实施 HPS 或者您是否会考虑其他音高估计器，所以我m 将 Welch/Blackman-Tukey 结果留在那里。）

---

原创我写了这篇评论，但不得不不断修改它，因为它令人困惑，所以这里是一个迷你答案。

根据我对 this intro to HPS 的简短阅读，我认为在找到四个抽取的响应后，您的量级并没有正确。

你想要：

array[i] = sqrt(data[i] * Complex.conjugate(data[i]) *
                hps2[i] * Complex.conjugate(hps2[i]) *
                hps3[i] * Complex.conjugate(hps3[i]) *
                hps4[i] * Complex.conjugate(hps4[i]) *
                hps5[i] * Complex.conjugate(hps5[i])).X;

这使用sqrt(x * Complex.conjugate(x)) 技巧找到x 的震级，然后将所有 5 个震级相乘。

（实际上，它会将sqrt 移到产品之外，所以你只做一个sqrt，节省一些时间，但结果相同。所以也许这是另一个技巧。）

最后一招：将结果取为实数部分，因为有时由于浮点精度问题，一个微小的虚数部分（如 1e-15）仍然存在。

完成此操作后，array 应该只包含真实的floats，您可以应用 max-bin-finding。

---

如果没有Conjugate 方法，那么老式的方法应该可以工作：

public float mag2(Complex c) { return c.X * c.X + c.Y * c.Y; }

// in HarmonicProductSpectrum 
array[i] = sqrt(mag2(data[i]) * mag2(hps2[i]) * mag2(hps3[i]) * mag2(hps4[i]) * mag2(hps5[i]));

您在下面的 cmets 中建议的两种方法存在代数缺陷，但上面应该是正确的。我不确定当您将 Complex 分配给浮点数时 C# 会做什么——也许它使用了真正的组件？我原以为那是编译器错误，但使用上面的代码，您对复杂数据做正确的事情，并且只将 float 分配给 array[i]。

Answer 2

要获得音高估计值，您必须将总的 bin 频率估计值除以用于该总和的下采样率。

补充：您还应该对幅度求和（abs()），而不是取复数和的幅度。

但是谐波积谱算法 (HPS)，尤其是在仅使用整数比率的下采样时，通常不能提供更好的音高估计分辨率。相反，它提供了比使用单个裸 FFT 幅度峰值更强大的粗略音高估计（不太可能被谐波欺骗），用于具有弱或缺少基本频谱内容的连续泛音丰富的音色。

如果您知道如何通过分数比（使用插值等）对频谱进行下采样，则可以尝试更细粒度的下采样，以便从 HPS 中获得更好的音高估计。或者，您可以使用 HPS 结果通知您使用另一种音高或频率估计方法搜索的较窄频率范围。