如何使用谐波乘积谱获取基频？

Question

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

    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

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

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

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

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

Answer 2

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

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

这只是底部的千赫兹，但频谱在 1 KHz 以上几乎是死的。

如果我将长音频剪辑分成 8192 个样本块（4096 个样本 50% 重叠）和 Hamming-window 每个块和运行 HPS，这是HPS 矩阵。这是整个数据集上 HPS 光谱的电影。基频好像挺稳定的。

full source code is here—有很多代码可以帮助对数据进行分块并可视化 HPS 运行ning 在分块上的输出，但是核心 HPS 函数从 def hps(…，是短。但它有几个技巧。

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

但假设您从纯正频率频谱开始，适当地获取其幅度，看起来您应该会得到合理的结果。尝试保存你的 HarmonicProductSpectrum 的输出，看看它是否像上面那样。

同样，完整的源代码位于 https://gist.github.com/fasiha/957035272009eb1c9eb370936a6af2eb。（我在那里尝试了另外几个频谱估计器，来自 Scipy 的 Welch 方法和我的 Blackman-Tukey 频谱估计器的端口。我不确定你是否打算实施 HPS 或者你是否会考虑其他音高估计器，所以我将 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 应该只包含真实的 float，您可以应用 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]));

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

如何使用谐波乘积谱获取基频？

How to get the fundamental frequency using Harmonic Product Spectrum?

c#

signal-processing

fft

pitch

audio-processing