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

我正试图从麦克风输入中获得音高。 首先,我通过FFT将信号从时域分解到频域。 在执行FFT之前,我已将Hamming窗口应用于信号。 然后我得到了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  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的谐波。

非常感谢帮助。

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

这是我将谐波产品频谱应用于完整数据集时所看到的,汉明窗口,具有5个下采样 – 乘法阶段:

完整数据

这只是底部千赫兹,但频谱几乎已经超过1 KHz。

如果我将长音频片段分成8192个样本块(4096样本50%重叠)和Hamming-window每个块并在其上运行HPS,这就是HPS的矩阵。 这是一部关于整个数据集的HPS频谱的电影。 基本频率似乎相当稳定。

0-500赫兹

完整的源代码在这里 – 有很多代码可以帮助分块数据并可视化运行在块上的HPS的输出,但核心HPSfunction,从def 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的结果。)


原来我把它写成评论但是不得不继续修改它因为它令人困惑所以这里它是一个迷你答案。

根据我对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 cX * cX + cY * cY; } // in HarmonicProductSpectrum array[i] = sqrt(mag2(data[i]) * mag2(hps2[i]) * mag2(hps3[i]) * mag2(hps4[i]) * mag2(hps5[i])); 

您在下面的评论中提出了两种方法的代数缺陷,但上述内容应该是正确的。 当你将一个Complex分配给一个浮点数时,我不确定C#会做什么 – 也许它使用了真实的组件? 我曾经认为这是一个编译器错误,但是使用上面的代码,你正在使用复杂的数据做正确的事情,并且只为array[i]分配一个float

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

补充:您还应该对幅度(abs())求和,而不是取复数和的幅度。

但是谐波乘积谱算法(HPS),特别是当仅使用整数比率的下采样时,通常不提供更好的音调估计分辨率。 相反,它提供了更强大的粗略音调估计(不太可能被谐波愚弄),而不是使用单个裸FFT幅度峰值来获得具有弱或基本频谱内容缺失的连续泛音丰富音色。

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