基频检测算法详细总结

基频估计算法 F0 estimate methods

基音频率，简称基音或基频，是语音信号最重要的韵律特征参数之一，基频的复杂之处在于不同的人说同一个字词，发出的语音基频可能不同，同一个人在不同的时间说相同的字词时基频也可能不同。通常，基频与发音人声带的长短、薄厚、韧性、劲度和发音习惯等有关系，在很大程度上代表了个人的特征。此外，基频还随着人的性别、年龄不同而有所不同。基频也受说话人情绪，习惯，环境等的影响。

一般来说，男性说话者的基音频率较低，而女性说话者和小孩的基音频率相对较高。可见，基频是随发音人、发音内容、发音时间等因素动态变化的。在汉语语音中，基频的变化主要包含了声调信息，声调具有辅助辨义的作用。另外，普通话中存在着变调和协同发音现象，同一个字在不同的语流中具有不同的声调。因此准确可靠地进行基频检测对汉语语音信号的处理显得尤为重要。

从上世纪末到现在，基频提取方法一共出现了几十乃至上百种。但目前为止，尚未找到一个完善的基频检测方法可以对各类人群（男女老少），不同语种，各类应用领域和各种环境条件都能通用的检测算法。

估计F0的方法可以分为两大类，传统信号处理方法和机器学习方法。传统信号处理方法又可分为：基于时域、基于频域、基于倒频域，或混合方法。但是如何评估比较这些算法还没有一个客观公正的标准。根据大量的应用实例的需求，对一个基频检测与提取算法的评价主要包括如下几个方面：

• 稳定性：无论在高信噪比还是低信噪比的情况下，算法都能检测和提取到语音信号的基频。

• 精准性：检测和提取到的语音信号的基频比较精准。

• 实时性：时延小，算法能够几乎实时进行处理。

• 复杂性：算法容易在不同的计算平台上实现，特别是在计算能力较弱的移动平台和嵌入式系统上。另外，算法中应该尽量少的包含各种变换和逆变换。

考虑到基频并不是每一帧都会有，在提取基频之前、之后或同时，往往要判断一下基频的有无，这称为清浊音判别（unvoiced / voiced decision，简称 U/V decision）。逐帧提取的基频常常含有错误，其中最常见的错误是「倍频错误」和「半频错误」，即提取出的基频是真实基频的两倍或者一半。为了纠正这些错误，通常会再对结果进行平滑（smoothing）处理，得到光滑的基音轨迹（pitch contour）。

基于规则的方法（rule-based）

也就是纯信号处理的方法。

既可以在频域上实现，也可以在时域上实现，还有倒谱域等。频域上的算法一般稳定性较好，但精准性较差(一般误差在10%左右)；时域算法与频域算法正好相反。

算法都包含如下三个主要步骤：
1.预处理：滤波，加窗分帧等
2.搜寻：可能的基频值F0（候选值）
3.跟踪：选择最可能的F0轨迹(这很重要，因为每个时刻我们都有几个候选轨迹)

基于时域

时域法以声音的波形为输入，其基本原理是寻找波形的最小正周期。当然，实际信号的周期性只能是近似的。以时域的自相关函数(ACF)法和平均幅度差函数(AMDF)法最为经典。

在应用时域方法之前，对信号进行滤波以只留下低频。设置阈值:最小频率和最大频率，例如75 ~ 500hz。F0估计仅对谐波。那些有停顿或噪声斑点的部分不仅没有观察意义，还可能改变相邻帧，并在应用插值或平滑时导致错误。所选帧的长度使其至少包含三个基频周期。

1. 自相关（auto-correlation function, ACF）
主要的算法就是自相关，随后又出现了一系列在此基础上的改进算法。自相关方法非常简单：计算自相关函数，然后定义它的第一个最大值，该最大值将反映信号中最显著的频率成分。那么，在使用自相关的问题是什么，为什么第一个最大值不能总是对应所需的频率?即使在高质量录音的近乎完美的条件下，由于信号的复杂结构，算法也容易出现错误。在接近真实的情况下，错误的数量急剧增加。此外，在最初质量较差和有噪声的录音中，我们可能会面临缺乏所期望的峰值等问题。

计算自相关函数前，为了有效地去除共振峰结构对短时自相关函数具体形状的影响，我们使用两个预处理过程。首先，使用60-90Hz的低通滤波器BPF，其截止频率在1700Hz，其次，在经过两个非线性处理单元NL1和NL2（可以相同或不同），这里用的是电平削波函数。

CL是削波的域值，通过滤波和削波可以有效消除共振峰和一些能量较小噪声的影响，使得自相关波形的峰值更加的突出。这样在进行峰值搜索的时候，误差可以减少很多。
最后计算自相关函数‚确定自相管函数的最大峰位置，得到基音频率：

2. 短时平均幅度差（average magnitude difference function, AMDF）
即把平移后的信号与原信号相减。对于相减的结果，可以求绝对值再积分。

3. YIN算法
以上两种方法的改进。求平方再积分，YIN 的论文把这样做的结果简称为差函数（difference function）：

差函数的谷值对应的平移量τ，可以代表信号在t处的周期。
YIN 的论文指出，做减法优于做乘法。这是因为，如果信号幅度在逐渐增大，那么把信号往后平移得越多，与原信号相乘再积分的结果就会越大。这将造成估计出的周期过大，基频过低，产生「半频错误」。做减法则没有这个问题。
YIN这个名字本身指的是自相关的便利和不准确之间的平衡：YIN这个名字来自东方哲学的阴YIN和阳YANG，暗指它所涉及的自相关和抵消之间的相互作用。它表明算法的核心思想是在差函数上寻找「谷值」。之所以想到用「阴阳相生」来形容差函数与自相关函数的关系，是因为差函数可以用自相关函数表示出来：

下图：一段周期约 150 个采样点的语音信号（上）及其差函数（中）与 CMNDF（下），取自 YIN 论文中的图 1、3

一个典型信号的差函数如图 3（中）所示，我们要找的是它尽量靠左，又比较深的谷。然而，差函数在原点（平移量为 0）时会取得最小值 0，这个谷会产生干扰，应该避开它。一种方法是设定周期的最小允许值，只在这个阈值右边寻找谷。然而问题是，原点处的谷的宽度并不固定，导致阈值难以设置。YIN 提出的解决办法是，在差函数的基础上，定义一种累积均值归一化差函数（cumulative mean normalized difference function, CMNDF）：

这相当于说，把差函数在τ处的值，用差函数在τ左边的平均值去归一化。CMNDF有两个好处：

• 上图可以看到，CMNDF 去除了差函数在原点处的谷，避免了它的干扰；
• CMNDF 的纵坐标是归一化的，不再依赖于信号本身的幅度。

若在 CMNDF 上直接求最小值，容易产生一个问题：由于信号的周期性并不理想，最小值可能并不是最左边的深谷，而是第二个甚至第三个。这会导致「半频错误」。YIN 的解决办法是，设置一个谷深阈值（如 0.1，图（下）里的虚线），取比此阈值更深的第一个谷的位置作为周期。CMNDF 「归一化」的特点使得选取谷深阈值成为可能。谷深阈值的选取带有任意性，不过 YIN 的作者用实验表明，YIN 的性能对此阈值并不敏感。若没有比阈值更深的谷，则选择最深的谷的位置作为基频，也可以认为这一帧没有基频。

以上就是 YIN 的核心方法。在此基础上还有一些非核心的改进，比如在 CMNDF 的每个谷处进行二次插值，以确定谷底的精确位置。
　　
代码实现：
我们可以通过使用python中的librosa库的yin函数来进行YIN算法提取音频。关于yin函数详细信息可见：https://librosa.org/doc/latest/generated/librosa.yin.html?highlight=yin

import matplotlib.pyplot as plt
import librosa
import numpy as np

for i in range(1, 11):
    filename = "Flowtron_sigma1.0/times" + str(i) +  "_sid0_sigma1.0.wav"
    y, sr = librosa.load(filename, sr=None)
    f0 = librosa.yin(y, fmin=80, fmax=400)
    f0[np.isnan(f0)] = 0
    times = librosa.times_like(f0)
    plt.plot(times, f0, "_", linewidth=1)
    plt.xlabel("Time(s)")
    plt.ylabel("F0")

基于频域

频域法则会先对信号做傅里叶变换，得到频谱（仅取幅度谱，舍弃相位谱）。频谱上在基频的整数倍处会有尖峰，频域法的基本原理就是要求出这些尖峰频率的最大公约数。

当我们谈到频域时，最突出的方面似乎是信号的谐波结构。换句话说，谱峰的频率可以被F0整除。在倒谱分析的帮助下，可以将这种周期性模式转换为一个明显的峰值。倒谱是估计功率谱对数的傅里叶变换(FFT)。倒谱峰对应于频谱中最具周期性的分量。

1. SWIPE算法

用频域法提取基频的一般原理为：对每一个候选频率 f，设计一个函数来表达 f 是基频的可能性，称为显著度函数（salience function）；使得显著度函数取得最大值的频率 f 就当作基频。

考虑到信号的幅度谱|x(f)|在基频 f 的整数倍处会有峰值，于是可以用幅度谱在 f 的各个倍数处的值的总和来代表 f 的显著度：

这称为谐波相加法（harmonic summation），它构成了大部分频域法的骨架。在具体的算法中，求和可能只取固定的项数，也可能一直取到信号所能表示的频率上限（即采样率的一半，称为奈奎斯特频率）；求和的各项也可能乘以不同的权重wk ，一般k越小，权重越大。

• 半频错误：若采用不加权的谐波相加法，并在求和时取到频率上限，就容易产生半频错误：这是因为f0的整数倍，也都是f0/2的整数倍，这会使得f0/2取得更大的显著度。
• 倍频错误：若是在求和时施加随k递减的权重，则又容易产生倍频错误，因为f0的 k 次谐波对于2f0来说也是k/2次谐波，权重更大，这会使得2f0取得更大的显著度。

SWIPE 的论文提出了如下核心论点：一个频率 f 若是好的基频候选值，那么幅度谱不仅要在 f 的整数倍处取得峰值，还要在 f 的半整数倍处取得谷值。SWIPE 用如下的函数衡量幅度谱在 f 处的峰值相对于紧邻它的两个谷值有多显著，称为峰谷距（peak-to-valley distance，PVD）：

中间峰值，减去两边谷值，则该式值最大

各次谐波的平均峰谷距（average PVD, APVD），就可以代表 f 的显著度：

平均峰谷距还可以表达成幅度谱与一个核函数的内积，这个核函数在 f 的整数倍处有向上的冲激，在 f 的半整数倍处有向下的冲激，如下图所示。为了表达「平均」峰谷距而不是「总」峰谷距，核函数需要经过某种归一化。

信号的幅度谱（虚线）与核函数（实线）做内积，取自 SWIPE 论文中的图 1：

平均峰谷距对于对抗倍频错误有奇效：当候选基频 f 等于真实基频 f0 的两倍时，f0 的奇次谐波将与核函数中向下的冲激对齐，这会减小显著度。它也能对抗半频错误：如果 f=f0/2 ，那么核函数中第奇数个向上的冲激将无法对应到频谱上的峰值，从而被「浪费」掉，这也会减小显著度。

SWIPE 算法在平均峰谷距的基础上，考虑了许多细节，让算法的性能精益求精。这些细节包括：

在典型信号中，锯齿波（sawtooth wave）的幅度谱与元音一样，都呈现 1/k 的衰减规律，这正是 SWIPE 算法全名（sawtooth wave inspired pitch estimator）的来源。不过，与其说是受到了「锯齿波」的启发，不如说算法是针对「语音中的元音」而设计的。

考虑了上述 6 条细节的 SWIPE 算法，可以表示成下面图 5 中的幅度谱与核函数的内积：
SWIPE 算法中实际使用的、经过弯折的幅度谱（虚线）与核函数（实线），取自 SWIPE 论文的图 7。

为了进一步减少「半频错误」，SWIPE 算法的作者又提出了这样一条改进（原文 II.I 节）：在计算平均峰谷距中的求和时，只考虑候选频率及其素数倍频率。改进后的算法称为 SWIPE’，撇号在英语中读作 prime，正是「素数」的意思。SWIPE’ 算法相当于把核函数改成了下面图的样子。这样，如果候选频率 是真实基频的整数分之一，那么幅度谱上的各次谐波的峰（基波除外）都将位于 的合数倍处，与核函数的峰完美错开，只能得到很小的显著度。

SWIPE算法的核函数（实线），只保留了第一个和第素数个峰。取自 SWIPE 论文中的图 9

混合法

有一个相当有说服力的名字YAAPT Yet Another algorithm of Pitch Tracking。它被归类为混合，因为它同时使用时间和频率数据。

YAAPT包括几个阶段，首先预处理开始。在这一阶段，初始信号的值被平方。这一步与YIN中的累积平均归一化差函数的目标相同:放大和恢复自相关的阻塞峰值。两个版本的信号都经过滤波，通常在50－1500 Hz或50－ 900 Hz的频谱范围内。
然后，根据变换后信号的频谱，计算出F0的基本轨迹。利用谱谐波相关函数(Spectral Harmonics Correlation，SHC)确定F0的候选值。

其中S(t,f)是第ｔ帧和频率ｆ的幅度频谱，WL是窗口的长度，单位为Hz, NH是谐波的数量(本文作者建议使用前三个谐波)。
基于功率谱对有音-无音帧进行了定义。然后，考虑音高加倍/音高减半的可能性，寻找最优轨迹。
然后，对于初始信号和转换信号，用归一化交叉相关(Normalized Cross Correlation，NCCF)来确定F0的候选信号，而不是自相关。

下一步是评估所有可能的候选值并计算他们的权重。候选的权重不仅取决于NCCF的振幅峰值，还取决于它们与F0轨迹的接近程度，这也是由频谱决定的。因此，频域被认为是相当稳定的。
对于所有剩余的候选值对，转移成本矩阵计算转移成本，以找到最佳的可能轨迹。

例子
我们选取了几段男女声音的片段，既有中性的声音，也有带有情感色彩的声音，并将它们组合在一起。让我们看一下信号的频谱图，它的强度(橙色)和F0(蓝色)。你可以在Praat用Ctrl+O(打开从文件中读取)打开它，然后按查看&编辑。

从录音中可以清楚地看出，男性和女性在带感情说话的音调都更高。此外，男性情感言语的F0与女性的F0相似。

跟踪Tracking

在菜单中选择Analyze periodicity to Pitch (ac)，通过自相关估计F0。设置窗口允许为F0候选估计设置3个参数，以及为在所有候选中构建最可能的F0轨迹的寻径算法设置6个参数。
路径查找器结果可以通过单击OK和View &编辑Pitch文件。从图中可以看出，除了选定的轨道之外，还有其他频率较低的候选轨道。

让我们以两个提供音高跟踪的库为例。首先是aubio，其中默认算法是YIN。并利用YAAPT算法进行amfm_分解。将Praat中的F0值插入到一个单独的文件(PraatPitch.txt)中(您可以手动操作:选择音频文件，单击“查看和编辑”，选择整个文件，然后在上面的菜单中选择“Pitch-Pitch列表”)。

现在比较这三种算法(YIN, YAAPT, Praat)的结果。

import amfm_decompy.basic_tools as basic
import amfm_decompy.pYAAPT as pYAAPT
import matplotlib.pyplot as plt
import numpy as np
import sys from aubio
import source, pitch 
# load audio
signal = basic.SignalObj('/home/eva/Documents/papers/habr/media/audio.wav')
filename = '/home/eva/Documents/papers/habr/media/audio.wav'
# YAAPT pitches 
pitchY = pYAAPT.yaapt(signal, frame_length=40, tda_frame_length=40, f0_min=75, f0_max=600)
# YIN pitches
downsample = 1
samplerate = 0
win_s = 1764 // downsample # fft size
hop_s = 441 // downsample # hop size
s = source(filename, samplerate, hop_s)
samplerate = s.samplerate
tolerance = 0.8
pitch_o = pitch("yin", win_s, hop_s, samplerate) pitch_o.set_unit("midi")
pitch_o.set_tolerance(tolerance)
pitchesYIN = []
confidences = [] 
total_frames = 0
while True:
     samples, read = s()
     pitch = pitch_o(samples)[0]
     pitch = int(round(pitch))
     confidence = pitch_o.get_confidence()
     pitchesYIN += [pitch]
     confidences += [confidence]
     total_frames += read
     if read < hop_s:
         break          
# load PRAAT pitches
praat = np.genfromtxt('/home/eva/Documents/papers/habr/PraatPitch.txt', filling_values=0)
praat = praat[:,1]
# plot
fig, (ax1,ax2,ax3) = plt.subplots(3, 1, sharex=True, sharey=True, figsize=(12, 8))
ax1.plot(np.asarray(pitchesYIN), label='YIN', color='green')
ax1.legend(loc="upper right")
ax2.plot(pitchY.samp_values, label='YAAPT', color='blue')
ax2.legend(loc="upper right")
ax3.plot(praat, label='Praat', color='red')
ax3.legend(loc="upper right") plt.show()

正如我们所看到的，在默认参数下，YIN不如Praat，显示出一个非常平坦的轨迹，相对于Praat的值更低。此外，它完全失去了男声和女声之间的过渡，以及情感和非情感话语之间的过渡。
YAAPT在女性情感言语的高音测试中失败，但总体上表现出较好的效果。我们不能肯定地说为什么它的效果更好，但我们可以假设它与从三个来源提取候选值以及比YIN更准确地计算他们的权重有关。
将Praat作为语音处理的标准算法(因为它在研究人员中被广泛使用)，我们可以得出结论，YAAPT比YIN更可靠和准确，尽管我们的测试对它来说相当困难。

用纯信号处理方法来提取基频有两个优点：一是不需要训练数据，二是原理容易解释。但这些方法同样有一些明显的缺点：一是算法中总有许多成分显得很人为，二是实际数据中总会有设计算法时没有考虑到的情况，导致提取结果不准确。近几年出现了几种使用机器学习的基频提取方法，它们牺牲了可解释性，但换取了性能的进一步飞跃。

另一种针对汉语的混合法：

如果在语音序列S中选取了语音帧F：

1. 对语音帧F做快速傅里叶变换，即得到频谱序列P = FFT(F)。如果采样频率是16 kHz，FFT的长度是512，那么FFT的频率分辨率是31.25 Hz。

2. 在频谱序列P上提取两个峰值点PV1，PV2，一个最大值点PK，其中PV1是93 ~ 218 Hz之间的峰值点，PV2是218~ 375 Hz之间的峰值点，PK是93~375 Hz之间的最大值点。注意，PV1，PV2可能不存在，但PK一定存在。

3. 基频F0粗估计。如果PV1存在，则用PV1估计F0；否则，如果PV2存在，则用PV2估计F0；如果PV1，PV2都不存在，则用PK估计F0，得到基频F0的粗估计CF0。

4. 粗估计CF0的调整。针对粗估计CF0是基于峰值点PV1给出的，分为两种情况进行处理：第一种是PV2与PK重合，且PV1不是PV2的半频，则用PV2估计F0。第二种是如果PV2存在，且PV2的峰值比PV1的峰值大，则用PV2估计F0。

5. 对于基频的粗估计CF0，在时域进一步精准估计。首先建立频率Hz与时域的语音序列S的下标之间的对应关系，这种关系是非线性的，为易于实现，用分段线性函数进行了简化，共分为三段，即100 ~ 200 Hz，200~ 300 Hz，300~400 Hz，经过计算，可以得到如下分段函数关系式：
   Y=−0.80X+240
   Y=−0.27X+134
   Y=−0.13X+92
   其中，输入X为基频的粗估计CF0，输出Y为CF0在语音序列 S中对应的下标。如果CF0位于区间100~ 200 Hz，则用函数式(2.1)；如果CF0位于区间200 ~ 300 Hz，则用函数式(2.2)；如果CF0位于区间300~400 Hz，则用函数式(2.3)。

6. 根据基频的粗估计CF0以及上述分段函数关系式，得到CF0在语音序列S中对应的下标SI.在语音序列S中的下标SI附近搜索峰值点SV1，且在下标2 * SI附近搜索峰值点SV2。如果SV1和SV2都存在，则根据SV1和SV2的下标计算出CF0的精确估计F0。

7. 如果粗估计CF0与精确估计F0相差15%以上，直接以粗估计CF0代替精确估计F0。

8. 返回精确估计F0作为本帧语音的基频。

需要说明的是，本算法直接对原始语音进行处理，未做高频滤波。这是考虑到滤波操作使得原始语音波形的峰值点可能发生变化，在时域中提取基频时可能产生偏差。

基于倒谱域

前面博客已经介绍过到倒谱的概念和计算方法。倒谱法可以将包含有基音周期的声门脉冲倒谱和声道响应倒谱分离，所以只要寻找声门脉冲倒谱的最大值，就是该帧的基音周期。

基于线性预测

通过逆滤波法得到预测误差（残差），此时已经去除了共振峰响应，再对残差进行倒谱分析（把声道响应影响减到最小），或者自相关，寻找最大值即可。

基于机器学习方法

用纯信号处理方法，很难同时降低「倍频错误」和「半频错误」的发生率：在算法中努力降低一者，往往会造成另一者升高，呈现「按下葫芦起来瓢」的态势。因此，最近几年出现了一些使用机器学习方法提取基频的研究，它们通过数据驱动（data-driven）的方式，超越了纯信号处理方法的性能。用机器学习来提取基频，同样可以分为时域法和频域法。

部分参考相关资料：基频提取算法综述

疑问：
大家都好奇基于深度学习的基频提取是不有杀鸡用牛刀的感觉，大材小用是否有必要？
没有一种基频提取算法是万能的，不同算法适用于不同信号。哪些算法更适合汉语普通话语音基频的提取呢？
基频提取的准确性对语音合成的质量影响有多大？比如野点过多