基于时频精细分析的乐器分类方法

1.本发明涉及信号处理技术领域，具体涉及一种基于时频精细分析的乐器分类方法。


背景技术：

2.通过解析音频文件，获取音频的乐器特征，乐器识别可以快速精准地识别出其中所含的乐器种类，从而有利于将其应用到音乐智能归类、乐器教学、音乐创作等多个领域。
3.在过去的研究中，学者们针对乐器识别所提出的时频特征取得了一定的分类效果，但是这些特征的尺度往往很大，描述的特征往往是整个时域或者频域的特征。比如振幅包络描述的是一段时间范围内音频幅度的总体变化，但是它却没办法准确描绘单个时域基频周期之间的精细变化。又比如频谱质心这一特征，只是从整个频谱图的宏观维度刻画了频域能量的重心，却无法精细刻画频域里不同次数谐波之间以及某次谐波和附近非谐波之间的大小关系，换句话说，频谱质心并没有办法准确传递频谱结构的信息。然而，频谱结构对乐器音色往往起着决定性的作用。因此，本发明认为过去的许多特征研究忽视了音乐的时频特征的精细化表征，亟需解决。
4.时频精细分析的基础在于基频的精准估计，然而基频的精确估计和谐波的精确标记往往是个复杂的过程，既要结合复杂的理论，又要考虑现实音频数据的异常情况，如：频谱混叠、噪声叠加、频谱偏移等，对基频估计算法和谐波标记算法的适应性和鲁棒性有极高的要求。因此，用代码实现适应各种现实异常情况的基频精准估计和谐波精确标记从而为后面的研究打下平台基础是有一定挑战性的。


技术实现要素：

5.本发明的目的是为了解决现有技术中的上述缺陷，提供一种基于时频精细分析的乐器分类方法，首先对乐器音频作谐波精确估计并标记出其频谱各谐波精确位置以及计算出时域单周期基频变化序列，然后提取乐器音频的时频精细特征融合传统音频特征进行乐器分类。
6.本发明的目的可以通过采取如下技术方案达到：
7.一种基于时频精细分析的乐器分类方法，所述乐器分类方法包括以下步骤：
8.s1、输入乐器单音音频，精准估计音频基频；
9.s2、根据音频的基频大小进行fft频域谐波标记；
10.s3、根据音频的基频大小进行时域单周期基频序列估计；
11.s4、提取时频精细特征和梅尔频率倒谱系数特征；
12.s5、基于提取的视频精细特征和梅尔频率倒谱系数特征，进行特征相关性和变异系数评估，建立随机森林模型对乐器进行分类。
13.进一步地，所述步骤s1过程如下：
14.s101、采用短时能量法进行乐器单音音频端点检测，截取有效音频，短时能量的计
算公式如下：
[0015][0016]
式中s[n]是输入的乐器单音音频信号，w[n]是窗函数，n是s[n]和w[n]的时域索引，nf是信号帧的帧长，z是窗函数的平移量，ez是平移量为z时对应的短时能量；不同的平移量对应的ez形成短时能量序列；为乐器单音音频有效音频起止端点的短时能量设置阈值，当前帧的短时能量超过起始端点的阈值时记录当前帧的中点为有效音频的起始端点，当前帧的短时能量低于终止端点的阈值时记录当前帧的中点为有效音频的终止端点；
[0017]
s102、通过时域自相关粗略估计出有效音频x[n]的基频f0，离散信号自相关公式如下：
[0018][0019]
式子中τ代表滞后值，即平移量，将x[n]向左平移τ个单位的距离后与原信号相乘，最后对乘积结果进行累加求和即当前滞后值下的自相关值；当滞后值为零时自相关值最大；根据自相关值和不同滞后值的关系后找到除滞后值为零之外的最大极大值点，该最大极大值点与滞后值为零的点的滞后值之差即为基波周期的大致估计，从而估计出基频f0；
[0020]
s103、对音频信号作快速傅里叶变换，快速傅里叶变换简称fft，并记录当前fft最大值对应的频率f
max
；根据时域自相关结果判断f
max
是否是干扰成分；如果是则从fft频谱去除该成分再重新获得fft最大值对应的频率f
max
，重复上述过程，直至根据时域自相关结果认为f
max
不是干扰成分；
[0021]
s104、采用窄带谱能量频率估计计算快速傅里叶变换的最大峰值对应频率值的去噪频率值fc，排除噪声干扰；在fft频谱数组的索引从1开始的情况下使用窄带谱能量频率估计计算fft频谱的去噪频率值fc的公式如下：
[0022][0023]
其中，k0是快速傅里叶变换的最大峰值对应的fft索引，k是fft频谱的索引，|xk|是fft频谱中第k个点对应的fft模值，fs是音频信号的采样率，n是fft点数；
[0024]
s105、根据去噪频率值fc和自相关估算的基频f0，计算出最大峰值大约位于nm次谐波处；其中：
[0025][0026]
得到基频估计值为f1:
[0027][0028]
s106、根据基频估计值为f1对音频时域信号进行重采样，新的采样率fs的计算公式：
[0029]fs
＝f1×2m
ꢀꢀꢀꢀ
公式(6)
[0030]
其中m的取值要满足：
[0031][0032]
保证新采样率不至于过大而导致采样点数n无法覆盖足够的基波周期；将重采样后的时域信号重新从步骤s102开始执行，直至基频估计值f1和上一次循环的基频估计值相比变化幅度小于0.001或者循环次数达到50次限制，输出最后一次的基频估计值f1。
[0033]
进一步地，所述步骤s2中，根据音频的基频大小进行fft频域谐波标记过程如下：
[0034]
s201、采用短时能量法进行乐器音频端点检测，截取有效音频部分；
[0035]
s202、根据步骤s1的基频估计值f1对音频时域信号进行重采样，使音频的采样率为f1的幂次方，减小频谱泄露，并把fft频谱转换为对数谱；
[0036]
s203、确定下一个谐波峰搜索范围在fft频谱上的左端点l和右端点r，如果是前3次谐波，则计算公式如下：
[0037][0038][0039]
其中p是上一次谐波峰在fft频谱中的位置索引，如果标记第一次谐波峰的位置，则p＝1；
[0040]
如果是标记第4次谐波及以后的谐波，则左端点l和右端点r的计算公式如下：
[0041]
l＝p+bt-2
ꢀꢀꢀꢀ
公式(10)
[0042]
r＝p+1.09t1ꢀꢀꢀꢀ
公式(11)
[0043]
其中，bt是fft频谱上第三次谐波峰和第一次谐波峰的距离的一半；t1是fft频谱离当前标记谐波最近的两次谐波fft点数间隔；
[0044]
s204、在搜索范围内求出所有极大值点，并求出所有极大值点的最大值点作为当前谐波的所求峰值点；如果搜索范围内不存在极大值点，则当标记前四次谐波时扩展搜索范围的左右边界，标记其他次数谐波时只扩展右边界；扩展前两次谐波的搜索边界时，左端点l和右端点r要满足以下条件：
[0045][0046][0047]
扩展第3和4次谐波的搜索边界时，左端点l和右端点r要满足以下条件：
[0048]
l-p≥0.7t1ꢀꢀꢀꢀ
公式(14)
[0049]
r-p≤1.3t1ꢀꢀꢀꢀ
公式(15)
[0050]
扩展第5次及以上次数谐波的搜索边界时，右端点r要满足以下条件：
[0051]
r-p≤1.3t1ꢀꢀꢀꢀ
公式(16)
[0052]
然后重新在新搜索范围内求出所有极大值点的最大值点作为当前谐波的所求峰值点，记录下以最大值点为中心的0.1倍谐波间距内的所有fft对数幅度值和最大值点的索
引值，分别存入矩阵hwb和数列hwi，hwb的每一行存放每一次谐波峰对数值及其附近的频谱对数值；如果仍然找不到极大值点则继续扩大搜索范围，循环往复。
[0053]
进一步地，所述步骤s3中，根据音频的基频大小进行时域单周期基频序列估计的过程如下：
[0054]
s301、采用短时能量法进行乐器音频端点检测，截取有效音频部分；
[0055]
s302、使用高通滤波器滤除低频干扰；
[0056]
s303、使用互相关运算标记单基波周期的分界点，以音频信号中间某一点o为原点，向两边等距延伸，延伸距离均为上一次运算得到的单基波周期长度，延伸范围内的信号为移位向量，移位向量的长度为上一次运算得到的单基波周期长度的两倍，因此有自适应调整的特点；如果是第一次互相关运算，则单周期长度参考s2步骤的基频估计结果；以o点为原点，向左或向右取四个单基波周期长度的信号作为固定向量，其中的单基波周期长度也是上一次运算得到的单基波周期长度；把移位向量当作滑动窗口和固定向量做互相关运算，得到互相关值和相对位移的对应关系；从对应关系中记录互相关值的极值点序列和对应的相对位移序列；从相对位移序列挑选出合理值，以o点为基准计算下一个单基波周期分隔点，分隔点和o的距离即为单基波周期长度；循环往复，从而得到时域单周期基频序列，互相关r[τ]的计算公式如下：
[0057][0058]
式子中τ代表滞后值，即平移量，g[n]和v[n]是两个不同的信号，将v[n]向左平移τ个单位的距离后，与g[n]相乘，对乘积结果进行累加求和即当前滞后值下的相似度。
[0059]
进一步地，所述步骤s4中，提取时频精细特征和梅尔频率倒谱系数特征的过程如下：
[0060]
s401、提取频域精细特征，包括：对数值大于零谐波能量质心、谐波峰毛刺平均衰弱比、对数值大于零最大谐波次数、两次谐波间斜率带宽比、谐波峰谷比、谐波能量下降率、前6次谐波奇次谐波与偶次谐波能量比，一次谐波峰谷比，谐波间斜率带宽比下降率，频域对数值大于零谐波能量与总能量之比，谐波峰谷差比半谐波间距，一次谐波峰谷差比半谐波间距；
[0061]
从步骤s2获得的hwb中可以获得各次谐波峰最大值点的对数幅度数列p[nh]，并将对数幅度数组中的小于等于零的数置零，索引从1开始计算。对数值大于零谐波能量质心c的计算公式如下：
[0062][0063]
其中l
p
是p[nh]的长度，也是最大谐波次数，nh是当前谐波次数，也是p[nh]的索引；
[0064]
对数值大于零最大谐波次数就是对数幅度数组p[nh]中的小于等于零的数置零后非零值对应的最大索引，索引从1开始；
[0065]
谐波峰毛刺平均衰弱比的计算需要在步骤s2中fft频谱的对数频谱大于零的情况下，获得以谐波峰最大值为中心的一个谐波间距范围内的离谐波峰最大值最近的二十个非谐波峰极大值a[n
p
]，n
p
是a[n
p
]的索引，如果范围内极大值个数不够二十个则取全部；如果
极大值个数小于等于10则直接取a[n
p
]的均值作为谐波峰周围毛刺的平均幅度；如果极大值个数大于10则在a[n
p
]里寻找长度为10连续的极大值组，该极大值组的fft索引跨度包含谐波峰最大值所在的索引或者与谐波峰最大值所在的索引之间不存在其他极大值并且总和最大；求得该极大值组的均值作为谐波峰周围毛刺的平均幅度；令谐波峰周围毛刺的平均幅度为a，单个谐波的谐波峰毛刺衰弱比pdr公式如下：
[0066][0067]
其中p为当前谐波峰的最大幅度，得到多个谐波的谐波峰毛刺衰弱比后取平均即为谐波峰毛刺平均衰弱比；
[0068]
令步骤s2的fft频谱的对数谱为mag[k]，步骤s2所求hwi为mag[k]中谐波峰最大值所在的索引序列，简称为hi[nh]，数列的索引全部都从1开始，则第3和4次谐波间斜率带宽比sfwr计算公式为：
[0069][0070]
计算谐波间斜率带宽比下降率时，先计算对数谱中直流分量和一次谐波间斜率带宽比为r1，再计算第5次谐波到第8次谐波之间各个谐波间距的谐波间斜率带宽比并取均值为r2，如果不够8次谐波，则取至最高次谐波；谐波间斜率带宽比下降率sfwrd的计算公式如下：
[0071][0072]
对谐波峰谷比基于fft对数频谱中对数值大于零的部分计算，令mag[k]中小于等于零的元素均为零，得到一个新数列mag[k]；对于一个谐波峰，其有左右两个波谷，左边波谷值为v
p1
，其值为mag[k]在左边谐波间距的波谷范围内的均值；右边波谷值为v
p2
，其值为mag[k]在右边谐波间距的波谷范围内的均值；在m次谐波峰谷比的计算中v
p1
的计算公式如下：
[0073][0074]
其中l是波谷值计算的左边界，r是波谷值计算的右边界；
[0075][0076][0077]
在m次谐波峰谷比的计算中v2的计算公式如下：
[0078][0079][0080][0081]
已知各次谐波峰最大值点的对数幅度数列为p[nh]，第m次谐波峰谷比mr的计算公式如下：
[0082][0083]
其中p[m]表示第m次谐波峰最大值点的对数幅度。计算多次谐波的谐波峰谷比的均值作为谐波峰谷比特征参数；
[0084]
一次谐波峰谷比是计算第m次谐波峰谷比的特殊情况，其步骤类似；在一次谐波峰谷比的计算中v
p1
的计算公式如下：
[0085][0086][0087][0088]
在一次谐波峰谷比的计算中v
p2
的计算公式如下：
[0089][0090][0091][0092]
一次谐波峰谷比fhpvr的计算公式如下：
[0093][0094]
谐波能量下降率的计算分为对数值大于零最大谐波次数num大于四和小于等于四两种情况；第一种情况下，令head为前三次谐波的第二大对数幅度，其谐波次数为s，令tail为num-2、num-1、num次谐波中对数幅度大于零的谐波的均值，这种情况下谐波能量下降率hdr的计算公式为：
[0095][0096]
如果对数值大于零最大谐波次数num小于等于四，则令head为前num-1次谐波中的最大对数幅度值，且对应谐波次数为s，tail为num次谐波对数幅度值；这种情况下谐波能量下降率hdr的计算公式为：
[0097][0098]
前6次谐波奇次谐波与偶次谐波能量比的计算中，令矩阵hwb为hb
ij
,前6次谐波奇次谐波与偶次谐波能量比r
oe6
的计算公式如下：
[0099][0100]
其中矩阵的行列索引都是由1开始，wid是hb
ij
一行的长度；
[0101]
计算频域对数值大于零谐波能量与总能量之比时，令矩阵hwb中小于等于零的元素为零，得到新矩阵uphb
ij
，则频域对数值大于零谐波能量与总能量之比uhtr的计算公式如
下：
[0102][0103]
式中row和col分别是矩阵uphb
ij
的行数和列数；
[0104]
谐波峰谷差比频率间距的计算需要针对m次谐波计算出谐波谷的值v和谐波峰最大对数值p[m]的差的绝对值，差的绝对值与谐波谷和峰值的频率间距相比为m次谐波峰谷差比频率间距，针对多个谐波求此特征值取平均数即为最终的谐波峰谷差比频率间距值；一般对对数频谱中第五至八次谐波中谐波峰对数值大于零的谐波进行计算并取平均值；对于一个谐波峰，其有左右两个波谷，左边波谷值为v
p1
，右边波谷值为v
p2
；v
p1
和v
p2
的计算和m次谐波峰谷比中的算法一致；则m次谐波峰谷差比频率间距dfr的计算公式如下：
[0105][0106]
对多个谐波计算dfr并取均值即为谐波峰谷差比频率间距值；
[0107]
一次谐波峰谷差比频率间距fdfr是m次谐波峰谷差比频率间距的特殊情况，左右波谷值v
p1
和v
p2
计算和一次谐波峰谷比一致；其计算公式为：
[0108][0109]
第1和2次谐波能量比则是一次谐波最大对数幅度与二次谐波最大对数幅度之比；
[0110]
s402、提取时域精细特征，包括单周期基频序列标准差、单周期基频序列均值、单周期基频序列中位数、单周期基频序列均值与中位数距离、单周期基频序列单位时间过均值次数；令步骤s3所求得时域单周期基频序列为f[n
t
]，其长度为len，n
t
表示时域单周期的序号。单周期基频序列均值f
mean
的计算公式为：
[0111][0112]
单周期基频序列标准差f
std
的计算公式为：
[0113][0114]
对f[n
t
]进行排序，则单周期基频序列中位数f
median
的计算公式为：
[0115][0116]
单周期基频序列均值与中位数距离f
mm
的计算公式为：
[0117]fmm
＝|f
mean-f
median
|
ꢀꢀꢀꢀ
公式(45)
[0118]
单周期基频序列单位时间过均值次数nc的计算公式为：
[0119][0120]
其中t1是时域音频信号的持续时间，u(t)的表达式如下：
[0121][0122]
s403、提取13维梅尔频率倒谱系数静态特征：首先对输入的音频数据进行预处理，其次对音频作fft变换并对fft模值求平方值得到音频的能量谱图，然后利用梅尔滤波器组进行滤波得到每一个滤波器的能量输出，最后对滤波器输出作对数运算并通过离散余弦变换dct得到梅尔频率倒谱系数特征。
[0123]
进一步地，所述步骤s5中，进行特征相关性和变异系数评估，建立随机森林模型对乐器进行分类的过程如下：
[0124]
s501、计算特征向量各元素之间的皮尔逊相关系数并评估特征有效性；
[0125]
s502、计算特征向量各元素的变异系数并评估特征有效性；
[0126]
s503、将时频精细特征和梅尔频率倒谱系数特征融合成新特征向量，将融合特征送入随机森林模型中进行训练，训练时采用五折交叉验证；
[0127]
随机森林模型的处理过程如下：利用自助法对训练集进行多次、有放回地随机采样，采样得到的训练样本即为θ；
[0128]
对于每个训练样本θ，形成对应的决策树；在构造的每棵决策树中，从每个训练样本θ的q个特征中，随机选取q个特征进行分类，其中，q≤q，并且表示对向下取整操作；计算选取出的q个特征的基尼不纯度，并认为基尼不纯度最小的特征为最优的分类特征；根据所选的最优分类特征，将最优分类特征所在节点进行二分，从剩余特征中挑选出次优分类特征；重复以上步骤，直到θ剩余特征数量为0或者当前决策树的分类效果最佳；
[0129]
将得到的所有30棵决策树进行组合，得到随机森林模型，并通过投票法得到最终的分类结果。
[0130]
随机森林模型能够处理具有高维特征的输入样本，而且不需要降维；随机森林模型的基本单位是决策树，它通过组合多个弱分类器，最终结果通过投票或取均值，使得整体模型的结果具有较高的精确度和泛化性能。
[0131]
本发明相对于现有技术具有如下的优点及效果：
[0132]
(1)本发明提出多个频域精细特征，相比传统的频域特征的大尺度描述，专注于描述频谱精细结构，提出了多个指标表征频谱谐波峰、毛刺、谐波谷等的精细特征，更能精确表征决定音色的频谱结构。频域精细特征更能表征乐器演奏起振特性、乐器本身共鸣腔特性等引起的频谱的细微差别。比如小号和小提琴对数频谱图在谐波峰最大值周围存在大小与谐波峰最大值相比比例较大的毛刺，而钢琴的对数频谱图中谐波峰十分尖锐，并不存在这种情况。频域精细特征就可以比传统频域特征更准确地描述这种频谱细微差别。
[0133]
(2)本发明提出的时域精细特征相比传统的时域特征，精准标记了时域单基波周期的间隔并记录基频变化序列，着重描述基频的精细变化。对于不同的乐器，其单音音频的
基频稳定性存在差异性，比如靠唇部振动控制基频的小号，其基频稳定性比击弦后自由振动的钢琴弦要弱，本发明提出的时域精细特征就可以准确表征这种差别，而传统时域特征并不关注此特点。
[0134]
(3)本发明得到的输入融合特征的随机森林模型的乐器分类误差率比只用13维梅尔频率倒谱系数静态特征的随机森林模型降低5.7个百分点。
[0135]
(4)本发明设计了基于自相关算法和窄带谱能量频率估计的循环迭代基频估计算法进行基频估计，实现了在低信噪比的情况下精确估计基频的效果。根据所估基频本发明设计了基于自适应窗口搜索的谐波标记算法和时域单周期基频序列估计算法精确标记频谱的谐波位置和获取单周期基频序列，一个方面克服了谐波标记时谐波偏移的问题，实现谐波的精确标记，另一方面单周期基频序列基于连续的时域单周期，更能表征时域基频的连续变化。
附图说明
[0136]
此处所说明的附图用来提供对本发明的进一步理解，构成本技术的一部分，本发明的示意性实施例及其说明用于解释本发明，并不构成对本发明的不当限定。在附图中：
[0137]
图1是本发明实施例中的一种基于时频精细分析的乐器分类方法的流程图；
[0138]
图2是本发明实施例中基于自相关算法和窄带谱能量频率估计算法的循环迭代基频估计算法流程图；
[0139]
图3是本发明实施例中基于自适应窗口搜索的谐波标记算法的流程图；
[0140]
图4是本发明实施例中时域单周期基频序列估计算法的流程图。
具体实施方式
[0141]
为使本发明实施例的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0142]
实施例1
[0143]
图1是本发明实施例提供的一种基于时频精细分析的乐器分类方法的流程图。
[0144]
s1、输入乐器单音音频，精准估计音频基频；
[0145]
如图2是本发明实施例中基于自相关算法和窄带谱能量频率估计算法的循环迭代基频估计算法流程图；
[0146]
s101、采用短时能量法进行乐器单音音频端点检测，截取有效音频，短时能量的计算公式如下：
[0147][0148]
式中s[n]是输入的乐器单音音频信号，w[n]是窗函数，n是s[n]和w[n]的时域索引，nf是信号帧的帧长，z是窗函数的平移量，ez是平移量为z时对应的短时能量；不同的平移量对应的ez形成短时能量序列；为乐器单音音频有效音频起止端点的短时能量设置阈值，
当前帧的短时能量超过起始端点的阈值时记录当前帧的中点为有效音频的起始端点，当前帧的短时能量低于终止端点的阈值时记录当前帧的中点为有效音频的终止端点；
[0149]
在本实例中主要是进行单音音频信号起始点的检测，采用的短时能量阈值为所有音频帧的短时能量的最大值的0.3倍；
[0150]
s102、通过时域自相关粗略估计出有效音频x[n]的基频f0，离散信号自相关公式如下：
[0151][0152]
式子中τ代表滞后值，即平移量，将x[n]向左平移τ个单位的距离后与原信号相乘，最后对乘积结果进行累加求和即当前滞后值下的自相关值；当滞后值为零时自相关值最大。根据自相关值和不同滞后值的关系后找到除滞后值为零之外的最大极大值点，其与滞后值为零的点的滞后值之差即为基波周期的大致估计，从而估计出基频f0；
[0153]
s103、对音频信号作快速傅里叶变换(fast fourier transform，fft)，并记录当前fft最大值对应的频率f
max
；根据时域自相关结果判断f
max
是否是干扰成分；如果是则从fft频谱去除该成分再重新获得fft最大值对应的频率f
max
，重复上述过程，直至认为f
max
不是干扰成分；
[0154]
在本实例中自相关估计的基频f0小于400赫兹时fft点数是4096个点，大于等于400赫兹时fft点数是2048个点；
[0155]
s104、采用窄带谱能量频率估计计算快速傅里叶变换的最大峰值对应频率值的去噪频率值fc，排除噪声干扰；在fft频谱数组的索引从1开始的情况下使用窄带谱能量频率估计计算fft频谱的去噪频率值fc的公式如下：
[0156][0157]
其中，k0是快速傅里叶变换的最大峰值对应的fft索引，k是fft频谱的索引，|xk|是fft频谱中第k个点对应的fft模值，fs是音频信号的采样率，n是fft点数；
[0158]
s105、根据去噪频率值fc和自相关估算的基频f0，计算出最大峰值大约位于nm次谐波处；其中：
[0159][0160]
得到基频估计值为f1:
[0161][0162]
s106、根据基频估计值为f1对音频时域信号进行重采样，新的采样率fs的计算公式：
[0163]fs
＝f1×2m
ꢀꢀꢀꢀ
公式(6)
[0164]
其中m的取值要满足：
[0165][0166]
保证新采样率不至于过大而导致采样点数n无法覆盖足够的基波周期；将重采样后的时域信号重新从步骤s102开始执行，直至基频估计值f1和上一次循环的基频估计值相比变化幅度小于0.001或者循环次数达到50次限制，输出最后一次的基频估计值f1。
[0167]
s2、根据音频的基频大小进行fft频域谐波标记；
[0168]
图3是本发明实施例中基于自适应窗口搜索的谐波标记算法的流程；
[0169]
s201、采用短时能量法进行乐器音频端点检测，截取有效音频部分；在本实例中主要是进行单音音频信号起始点的检测，采用的短时能量阈值为所有音频帧的短时能量的最大值的0.3倍；
[0170]
s202、根据s1步骤的基频估计值f1对音频时域信号进行重采样，使音频的采样率为f1的幂次方，减小频谱泄露，并把fft频谱转换为对数谱；在本实施例中，由于需要精细化频谱分析，因此谐波标记时fft点数为16384个点。
[0171]
s203、确定下一个谐波峰搜索范围在fft频谱上的左端点l和右端点r，如果是前3次谐波，则计算公式如下：
[0172][0173][0174]
其中p是上一次谐波峰在fft频谱中的位置索引，如果标记第一次谐波峰的位置，则p＝1；
[0175]
如果是标记第4次谐波及以后的谐波，则左端点l和右端点r的计算公式如下：
[0176]
l＝p+bt-2
ꢀꢀꢀꢀ
公式(10)
[0177]
r＝p+1.09t1ꢀꢀꢀꢀ
公式(11)
[0178]
其中，bt是fft频谱上第三次谐波峰和第一次谐波峰的距离的一半；t1是fft频谱离当前标记谐波最近的两次谐波fft点数间隔；
[0179]
s204、在搜索范围内求出所有极大值点，并求出所有极大值点的最大值点作为当前谐波的所求峰值点；如果搜索范围内不存在极大值点，则当标记前四次谐波时扩展搜索范围的左右边界，标记其他次数谐波时只扩展右边界；扩展前两次谐波的搜索边界时，左端点l和右端点r要满足以下条件：
[0180][0181][0182]
扩展第3和4次谐波的搜索边界时，左端点l和右端点r要满足以下条件：
[0183]
l-p≥0.7t1ꢀꢀꢀꢀ
公式(14)
[0184]
r-p≤1.3t1ꢀꢀꢀꢀ
公式(15)
[0185]
扩展第5次及以上次数谐波的搜索边界时，右端点r要满足以下条件：
[0186]
r-p≤1.3t1ꢀꢀꢀꢀ
公式(16)
[0187]
然后重新在新搜索范围内求出所有极大值点的最大值点作为当前谐波的所求峰值点，记录下以最大值点为中心的0.1倍谐波间距内的所有fft对数幅度值和最大值点的索引值，分别存入矩阵hwb和数列hwi，hwb的每一行存放每一次谐波峰对数值及其附近的频谱对数值；如果仍然找不到极大值点则继续扩大搜索范围，循环往复。
[0188]
s3、根据音频的基频大小进行时域单周期基频序列估计的过程如下：
[0189]
图4是本发明实施例中时域单周期基频序列估计算法的流程图；
[0190]
s301、采用短时能量法进行乐器音频端点检测，截取有效音频部分；在本实例中主要是进行单音音频信号起始点的检测，采用的短时能量阈值为所有音频帧的短时能量的最大值的0.3倍；
[0191]
s302、使用高通滤波器滤除低频干扰。在本实例中截止频率设置为步骤s1的基频估计值f1的0.3倍。
[0192]
s303、使用互相关运算标记单基波周期的分界点，以音频信号中间某一点o为原点，向两边等距延伸，延伸距离均为上一次运算得到的单基波周期长度，延伸范围内的信号为移位向量，移位向量的长度为上一次运算得到的单基波周期长度的两倍，因此有自适应调整的特点；如果是第一次互相关运算，则单周期长度参考s2步骤的基频估计结果；以o点为原点，向左或向右取四个单基波周期长度的信号作为固定向量，其中的单基波周期长度也是上一次运算得到的单基波周期长度；把移位向量当做滑动窗口和固定向量做互相关运算，得到互相关值和相对位移的对应关系；从对应关系中记录互相关值的极值点序列和对应的相对位移序列；从相对位移序列挑选出合理值，以o点为基准计算下一个单基波周期分隔点，分隔点和o的距离即为单基波周期长度；循环往复，从而得到时域单周期基频序列，互相关r[τ]的计算公式如下：
[0193][0194]
式子中τ代表滞后值，即平移量，g[n]和v[n]是两个不同的信号，将v[n]向左平移τ个单位的距离后，与g[n]相乘，对乘积结果进行累加求和即当前滞后值下的相似度。
[0195]
s4、提取时频精细特征和梅尔频率倒谱系数特征的过程如下：
[0196]
s401、提取频域精细特征，包括：对数值大于零谐波能量质心、谐波峰毛刺平均衰弱比、对数值大于零最大谐波次数、两次谐波间斜率带宽比、谐波峰谷比、谐波能量下降率、前6次谐波奇次谐波与偶次谐波能量比，一次谐波峰谷比，谐波间斜率带宽比下降率，频域对数值大于零谐波能量与总能量之比，谐波峰谷差比半谐波间距，一次谐波峰谷差比半谐波间距；
[0197]
从步骤s2获得的hwb中可以获得各次谐波峰最大值点的对数幅度数列p[nh]，并将对数幅度数组中的小于等于零的数置零，索引从1开始计算。对数值大于零谐波能量质心c的计算公式如下：
[0198][0199]
其中l
p
是p[nh]的长度，也是最大谐波次数，nh是当前谐波次数，也是p[nh]的索引；
[0200]
对数值大于零最大谐波次数就是对数幅度数组p[nh]中的小于等于零的数置零后
非零值对应的最大索引，索引从1开始；
[0201]
谐波峰毛刺平均衰弱比的计算需要在步骤s2中fft频谱的对数频谱大于零的情况下，获得以谐波峰最大值为中心的一个谐波间距范围内的离谐波峰最大值最近的二十个非谐波峰极大值a[n
p
]，n
p
是a[n
p
]的索引，如果范围内极大值个数不够二十个则取全部；如果极大值个数小于等于10则直接取a[n
p
]的均值作为谐波峰周围毛刺的平均幅度；如果极大值个数大于10则在a[n
p
]里寻找长度为10连续的极大值组，该极大值组的fft索引跨度包含谐波峰最大值所在的索引或者与谐波峰最大值所在的索引之间不存在其他极大值并且总和最大；求得该极大值组的均值作为谐波峰周围毛刺的平均幅度；令谐波峰周围毛刺的平均幅度为a，单个谐波的谐波峰毛刺衰弱比pdr公式如下：
[0202][0203]
其中p为当前谐波峰的最大幅度，得到多个谐波的谐波峰毛刺衰弱比后取平均即为谐波峰毛刺平均衰弱比；
[0204]
令步骤s2的fft频谱的对数谱为mag[k]，步骤s2所求hwi为mag[k]中谐波峰最大值所在的索引序列，简称为hi[nh]，数列的索引全部都从1开始，则第3和4次谐波间斜率带宽比sfwr计算公式为：
[0205][0206]
计算谐波间斜率带宽比下降率时，先计算对数谱中直流分量和一次谐波间斜率带宽比为r1，再计算第5次谐波到第8次谐波之间各个谐波间距的谐波间斜率带宽比并取均值为r2，如果不够8次谐波，则取至最高次谐波；谐波间斜率带宽比下降率sfwrd的计算公式如下：
[0207][0208]
对谐波峰谷比基于fft对数频谱中对数值大于零的部分计算，令mag[k]中小于等于零的元素均为零，得到一个新数列mag[k]；对于一个谐波峰，其有左右两个波谷，左边波谷值为v
p1
，其值为mag[k]在左边谐波间距的波谷范围内的均值；右边波谷值为v
p2
，其值为mag[k]在右边谐波间距的波谷范围内的均值；在m次谐波峰谷比的计算中v
p1
的计算公式如下：
[0209][0210]
其中l是波谷值计算的左边界，r是波谷值计算的右边界；
[0211][0212][0213]
在m次谐波峰谷比的计算中v2的计算公式如下：
[0214]
[0215][0216][0217]
已知各次谐波峰最大值点的对数幅度数列为p[nh]，第m次谐波峰谷比mr的计算公式如下：
[0218][0219]
其中p[m]表示第m次谐波峰最大值点的对数幅度。计算多次谐波的谐波峰谷比的均值作为谐波峰谷比特征参数；
[0220]
一次谐波峰谷比是计算第m次谐波峰谷比的特殊情况，其步骤类似；在一次谐波峰谷比的计算中v
p1
的计算公式如下：
[0221][0222][0223][0224]
在一次谐波峰谷比的计算中v
p2
的计算公式如下：
[0225][0226][0227][0228]
一次谐波峰谷比fhpvr的计算公式如下：
[0229][0230]
谐波能量下降率的计算分为对数值大于零最大谐波次数num大于四和小于等于四两种情况；第一种情况下，令head为前三次谐波的第二大对数幅度，其谐波次数为s，令tail为num-2、num-1、num次谐波中对数幅度大于零的谐波的均值，这种情况下谐波能量下降率hdr的计算公式为：
[0231][0232]
如果对数值大于零最大谐波次数num小于等于四，则令head为前num-1次谐波中的最大对数幅度值，且对应谐波次数为s，tail为num次谐波对数幅度值；这种情况下谐波能量下降率hdr的计算公式为：
[0233][0234]
前6次谐波奇次谐波与偶次谐波能量比的计算中，令矩阵hwb为hb
ij
,前6次谐波奇次谐波与偶次谐波能量比r
oe6
的计算公式如下：
[0235][0236]
其中矩阵的行列索引都是由1开始，wid是hb
ij
一行的长度；
[0237]
计算频域对数值大于零谐波能量与总能量之比时，令矩阵hwb中小于等于零的元素为零，得到新矩阵uphb
ij
，则频域对数值大于零谐波能量与总能量之比uhtr的计算公式如下：
[0238][0239]
式中row和col分别是矩阵uphb
ij
的行数和列数；
[0240]
谐波峰谷差比频率间距的计算需要针对m次谐波计算出谐波谷的值v和谐波峰最大对数值p[m]的差的绝对值，差的绝对值与谐波谷和峰值的频率间距相比为m次谐波峰谷差比频率间距，针对多个谐波求此特征值取平均数即为最终的谐波峰谷差比频率间距值；一般对对数频谱中第5至8次谐波中谐波峰对数值大于零的谐波进行计算并取平均值；对于一个谐波峰，其有左右两个波谷，左边波谷值为v
p1
，右边波谷值为v
p2
；v
p1
和v
p2
的计算和m次谐波峰谷比中的算法一致；则m次谐波峰谷差比频率间距dfr的计算公式如下：
[0241][0242]
对多个谐波计算dfr并取均值即为谐波峰谷差比频率间距值；
[0243]
一次谐波峰谷差比频率间距fdfr是m次谐波峰谷差比频率间距的特殊情况，左右波谷值v
p1
和v
p2
计算和一次谐波峰谷比一致；其计算公式为：
[0244][0245]
第1和2次谐波能量比则是一次谐波最大对数幅度与二次谐波最大对数幅度之比；
[0246]
s402、提取时域精细特征，包括单周期基频序列标准差、单周期基频序列均值、单周期基频序列中位数、单周期基频序列均值与中位数距离、单周期基频序列单位时间过均值次数；令步骤s3所求得时域单周期基频序列为f[n
t
]，其长度为len，n
t
表示时域单周期的序号。单周期基频序列均值f
mean
的计算公式为：
[0247][0248]
单周期基频序列标准差f
std
的计算公式为：
[0249][0250]
对f[n
t
]进行排序，则单周期基频序列中位数f
median
的计算公式为：
[0251][0252]
单周期基频序列均值与中位数距离f
mm
的计算公式为：
[0253]fmm
＝|f
mean-f
median
|
ꢀꢀꢀꢀ
公式(45)
[0254]
单周期基频序列单位时间过均值次数nc的计算公式为：
[0255][0256]
其中t1是时域音频信号的持续时间，u(t)的表达式如下：
[0257][0258]
s403、提取13维梅尔频率倒谱系数静态特征：首先对输入的音频数据进行预处理，其次对音频作fft变换并对fft模值求平方值得到音频的能量谱图，然后利用梅尔滤波器组进行滤波得到每一个滤波器的能量输出，最后对滤波器输出作对数运算并通过离散余弦变换换(discrete cosine transform,dct)得到梅尔频率倒谱系数特征。
[0259]
s5、基于提取的视频精细特征和梅尔频率倒谱系数特征，进行特征相关性和变异系数评估，建立随机森林模型对乐器进行分类的过程如下：
[0260]
s501、计算特征向量各元素之间的皮尔逊相关系数并评估特征有效性；
[0261]
s502、计算特征向量各元素的变异系数并评估特征有效性；
[0262]
s503、将时频精细特征和梅尔频率倒谱系数特征融合成新特征向量，在本实施例中将5种乐器的822个音频的融合特征输入随机森林模型中进行训练，训练时采用五折交叉验证；再将其中的13维梅尔频率倒谱系数静态特征单独输入随机森林模型中进行训练，训练时采用五折交叉验证；随机森林能够处理具有高维特征的输入样本，而且不需要降维；随机森林的基本单位是决策树，它通过组合多个弱分类器，最终结果通过投票或取均值，使得整体模型的结果具有较高的精确度和泛化性能；本实例最终得到的输入融合特征的随机森林模型的乐器分类误差率比只用13维梅尔频率倒谱系数静态特征的随机森林模型降低6％，如表1所示。
[0263]
表1.特征种类及分类准确率表
[0264][0265][0266]
实施例2
[0267]
基于上述实施例1公开的一种基于时频精细分析的乐器分类方法，本实施例继续给出基于时频精细分析的乐器分类方法的分类实施过程，具体过程如下：
[0268]
s1、参照实施例1中对应步骤，获得输入的小提琴单音音频的基频f1；
[0269]
s2、记录下以各次谐波最大值点为中心的0.1倍谐波间距内的所有fft对数幅度值和最大值点的索引值，分别存入矩阵hwb和数列hwi，hwb的每一行存放每一次谐波峰对数值及其附近的频谱对数值；
[0270]
s3、参照实施例1中对应步骤，此处不再赘述；
[0271]
s4、提取时频精细特征和13维梅尔频率倒谱系数静态特征，采用matlab代码对各个特征进行自动计算，算得该音频的各个特征值如下表所示。
[0272]
表2.对应音频的时频精细特征和13维梅尔频率倒谱系数静态特征表
[0273]
[0274][0275]
s5、基于实例一得到的输入融合特征的随机森林模型，将以上特征输入该模型得到的分类结果是小提琴，分类正确。
[0276]
上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

技术特征：
1.一种基于时频精细分析的乐器分类方法，其特征在于，所述乐器分类方法包括以下步骤：s1、输入乐器单音音频，精准估计音频基频；s2、根据音频的基频大小进行fft频域谐波标记；s3、根据音频的基频大小进行时域单周期基频序列估计；s4、提取时频精细特征和梅尔频率倒谱系数特征；s5、基于提取的视频精细特征和梅尔频率倒谱系数特征，进行特征相关性和变异系数评估，建立随机森林模型对乐器进行分类。2.根据权利要求1所述的一种基于时频精细分析的乐器分类方法，其特征在于，所述步骤s1过程如下：s101、采用短时能量法进行乐器单音音频端点检测，截取有效音频，短时能量的计算公式如下：式中s[n]是输入的乐器单音音频信号，w[n]是窗函数，n是s[n]和w[n]的时域索引，n
f
是信号帧的帧长，z是窗函数的平移量，e
z
是平移量为z时对应的短时能量；不同的平移量对应的e
z
形成短时能量序列；为乐器单音音频有效音频起止端点的短时能量设置阈值，当前帧的短时能量超过起始端点的阈值时记录当前帧的中点为有效音频的起始端点，当前帧的短时能量低于终止端点的阈值时记录当前帧的中点为有效音频的终止端点；s102、通过时域自相关粗略估计出有效音频x[n]的基频f0，离散信号自相关公式如下：式子中τ代表滞后值，即平移量，将x[n]向左平移τ个单位的距离后与原信号相乘，最后对乘积结果进行累加求和即当前滞后值下的自相关值；当滞后值为零时自相关值最大；根据自相关值和不同滞后值的关系后找到除滞后值为零之外的最大极大值点，该最大极大值点与滞后值为零的点的滞后值之差即为基波周期的大致估计，从而估计出基频f0；s103、对音频信号作快速傅里叶变换，快速傅里叶变换简称fft，并记录当前fft最大值对应的频率f
max
；根据时域自相关结果判断f
max
是否是干扰成分；如果是则从fft频谱去除该成分再重新获得fft最大值对应的频率f
max
，重复上述过程，直至根据时域自相关结果认为f
max
不是干扰成分；s104、采用窄带谱能量频率估计计算快速傅里叶变换的最大峰值对应频率值的去噪频率值f
c
，排除噪声干扰；在fft频谱数组的索引从1开始的情况下使用窄带谱能量频率估计计算fft频谱的去噪频率值f
c
的公式如下：其中，k0是快速傅里叶变换的最大峰值对应的fft索引，k是fft频谱的索引，|x
k
|是fft
频谱中第k个点对应的fft模值，f
s
是音频信号的采样率，n是fft点数；s105、根据去噪频率值f
c
和自相关估算的基频f0，计算出最大峰值大约位于n
m
次谐波处；其中：得到基频估计值为f1:s106、根据基频估计值为f1对音频时域信号进行重采样，新的采样率f
s
的计算公式：f
s
＝f1×2m
ꢀꢀꢀ
公式(6)其中m的取值要满足：将重采样后的时域信号重新从步骤s102开始执行，直至基频估计值f1和上一次循环的基频估计值相比变化幅度小于0.001或者循环次数达到50次限制，输出最后一次的基频估计值f1。3.根据权利要求1所述的一种基于时频精细分析的乐器分类方法，其特征在于，所述步骤s2中，根据音频的基频大小进行fft频域谐波标记过程如下：s201、采用短时能量法进行乐器音频端点检测，截取有效音频部分；s202、根据步骤s1的基频估计值f1对音频时域信号进行重采样，使音频的采样率为f1的幂次方，并把fft频谱转换为对数谱；s203、确定下一个谐波峰搜索范围在fft频谱上的左端点l和右端点r，如果是前3次谐波，则计算公式如下：波，则计算公式如下：其中p是上一次谐波峰在fft频谱中的位置索引，如果标记第一次谐波峰的位置，则p＝1；如果是标记第4次谐波及以后的谐波，则左端点l和右端点r的计算公式如下：l＝p+bt-2
ꢀꢀꢀ
公式(10)r＝p+1.09t1ꢀꢀ
公式(11)其中，bt是fft频谱上第3次谐波峰和第1次谐波峰的距离的一半；t1是fft频谱离当前标记谐波最近的两次谐波fft点数间隔；s204、在搜索范围内求出所有极大值点，并求出所有极大值点的最大值点作为当前谐波的所求峰值点；如果搜索范围内不存在极大值点，则当标记前4次谐波时扩展搜索范围的左右边界，标记其他次数谐波时只扩展右边界；扩展前两次谐波的搜索边界时，左端点l和右端点r要满足以下条件：
扩展第3和4次谐波的搜索边界时，左端点l和右端点r要满足以下条件：l-p≥0.7t1ꢀꢀ
公式(14)r-p≤1.3t1ꢀꢀ
公式(15)扩展第5次及以上次数谐波的搜索边界时，右端点r要满足以下条件：r-p≤1.3t1ꢀꢀ
公式(16)然后重新在新搜索范围内求出所有极大值点的最大值点作为当前谐波的所求峰值点，记录下以最大值点为中心的0.1倍谐波间距内的所有fft对数幅度值和最大值点的索引值，分别存入矩阵hwb和数列hwi，hwb的每一行存放每一次谐波峰对数值及其附近的频谱对数值；如果仍然找不到极大值点则继续扩大搜索范围，循环往复直到找到极大值点。4.根据权利要求1所述的一种基于时频精细分析的乐器分类方法，其特征在于，所述步骤s3中，根据音频的基频大小进行时域单周期基频序列估计的过程如下：s301、采用短时能量法进行乐器音频端点检测，截取有效音频部分；s302、使用高通滤波器滤除低频干扰；s303、使用互相关运算标记单基波周期的分界点，以音频信号中间某一点o为原点，向两边等距延伸，延伸距离均为上一次运算得到的单基波周期长度，延伸范围内的信号为移位向量，移位向量的长度为上一次运算得到的单基波周期长度的两倍，因此有自适应调整的特点；如果是第一次互相关运算，则单周期长度取步骤s2的基频估计结果；以o点为原点，向左或向右取四个单基波周期长度的信号作为固定向量，其中的单基波周期长度也是上一次运算得到的单基波周期长度；把移位向量当作滑动窗口和固定向量做互相关运算，得到互相关值和相对位移的对应关系；从对应关系中记录互相关值的极值点序列和对应的相对位移序列；从相对位移序列挑选出合理值，以o点为基准计算下一个单基波周期分隔点，分隔点和o的距离即为单基波周期长度；循环往复，从而得到时域单周期基频序列，互相关r[τ]的计算公式如下：式子中τ代表滞后值，即平移量，g[n]和v[n]是两个不同的信号，将v[n]向左平移τ个单位的距离后，与g[n]相乘，对乘积结果进行累加求和即当前滞后值下的相似度。5.根据权利要求4所述的一种基于时频精细分析的乐器分类方法，其特征在于，所述步骤s4中，提取时频精细特征和梅尔频率倒谱系数特征的过程如下：s401、提取频域精细特征，包括：对数值大于零谐波能量质心、谐波峰毛刺平均衰弱比、对数值大于零最大谐波次数、两次谐波间斜率带宽比、谐波峰谷比、谐波能量下降率、前6次谐波奇次谐波与偶次谐波能量比，一次谐波峰谷比，谐波间斜率带宽比下降率，频域对数值大于零谐波能量与总能量之比，谐波峰谷差比半谐波间距，一次谐波峰谷差比半谐波间距；从步骤s2获得的hwb中可以获得各次谐波峰最大值点的对数幅度数列p[n
h
]，并将对数
幅度数组中的小于等于零的数置零，索引从1开始计算，对数值大于零谐波能量质心c的计算公式如下：其中l
p
是p[n
h
]的长度，也是最大谐波次数，n
h
是当前谐波次数，也是p[n
h
]的索引；对数值大于零最大谐波次数就是对数幅度数组p[n
h
]中的小于等于零的数置零后非零值对应的最大索引，索引从1开始；谐波峰毛刺平均衰弱比的计算需要在步骤s2中fft频谱的对数频谱大于零的情况下，获得以谐波峰最大值为中心的一个谐波间距范围内的离谐波峰最大值最近的20个非谐波峰极大值a[n
p
]，n
p
是a[n
p
]的索引，如果范围内极大值个数不够20则取全部；如果极大值个数小于等于10则直接取a[n
p
]的均值作为谐波峰周围毛刺的平均幅度；如果极大值个数大于10则在a[n
p
]里寻找长度为10连续的极大值组，该极大值组的fft索引跨度包含谐波峰最大值所在的索引或者与谐波峰最大值所在的索引之间不存在其他极大值并且总和最大；求得该极大值组的均值作为谐波峰周围毛刺的平均幅度；令谐波峰周围毛刺的平均幅度为a，单个谐波的谐波峰毛刺衰弱比pdr公式如下：其中p为当前谐波峰的最大幅度，得到多个谐波的谐波峰毛刺衰弱比后取平均即为谐波峰毛刺平均衰弱比；令步骤s2得到的fft频谱的对数谱为mag[k]，步骤s2所求hwi为mag[k]中谐波峰最大值所在的索引序列，简称为hi[n
h
]，数列的索引全部都从1开始，则第3和4次谐波间斜率带宽比sfwr计算公式为：计算谐波间斜率带宽比下降率时，先计算对数谱中直流分量和一次谐波间斜率带宽比为r1，再计算第5次谐波到第8次谐波之间各个谐波间距的谐波间斜率带宽比并取均值为r2，如果不够8次谐波，则取至最高次谐波；谐波间斜率带宽比下降率sfwrd的计算公式如下：对谐波峰谷比基于fft对数频谱中对数值大于零的部分计算，令mag[k]中小于等于零的元素均为零，得到一个新数列mag[k]；对于一个谐波峰，其有左右两个波谷，左边波谷值为v
p1
，其值为mag[k]在左边谐波间距的波谷范围内的均值；右边波谷值为v
p2
，其值为mag[k]在右边谐波间距的波谷范围内的均值；在m次谐波峰谷比的计算中v
p1
的计算公式如下：其中l是波谷值计算的左边界，r是波谷值计算的右边界；
在m次谐波峰谷比的计算中v2的计算公式如下：在m次谐波峰谷比的计算中v2的计算公式如下：在m次谐波峰谷比的计算中v2的计算公式如下：已知各次谐波峰最大值点的对数幅度数列为p[n
h
]，第m次谐波峰谷比mr的计算公式如下：其中p[m]表示第m次谐波峰最大值点的对数幅度。计算多次谐波的谐波峰谷比的均值作为谐波峰谷比特征参数；一次谐波峰谷比是计算第m次谐波峰谷比的特殊情况，步骤类似；在一次谐波峰谷比的计算中v
p1
的计算公式如下：的计算公式如下：的计算公式如下：在一次谐波峰谷比的计算中v
p2
的计算公式如下：的计算公式如下：的计算公式如下：一次谐波峰谷比fhpvr的计算公式如下：谐波能量下降率的计算分为对数值大于零最大谐波次数num大于四和小于等于四两种情况；第一种情况下，令head为前三次谐波的第二大对数幅度，其谐波次数为s，令tail为num-2、num-1、num次谐波中对数幅度大于零的谐波的均值，这种情况下谐波能量下降率hdr的计算公式为：如果对数值大于零最大谐波次数num小于等于四，则令head为前num-1次谐波中的最大对数幅度值，且对应谐波次数为s，tail为num次谐波对数幅度值；这种情况下谐波能量下降率hdr的计算公式为：
前6次谐波奇次谐波与偶次谐波能量比的计算中，令矩阵hwb为hb
ij
,前6次谐波奇次谐波与偶次谐波能量比r
oe6
的计算公式如下：其中矩阵的行列索引都是由1开始，wid是hb
ij
一行的长度；计算频域对数值大于零谐波能量与总能量之比时，令矩阵hwb中小于等于零的元素为零，得到新矩阵uphb
ij
，则频域对数值大于零谐波能量与总能量之比uhtr的计算公式如下：式中row和col分别是矩阵uphb
ij
的行数和列数；谐波峰谷差比频率间距的计算需要针对m次谐波计算出谐波谷的值v和谐波峰最大对数值p[m]的差的绝对值，差的绝对值与谐波谷和峰值的频率间距相比为m次谐波峰谷差比频率间距，针对多个谐波求此特征值取平均数即为最终的谐波峰谷差比频率间距值；一般对对数频谱中第五至八次谐波中谐波峰对数值大于零的谐波进行计算并取平均值；对于一个谐波峰，其有左右两个波谷，左边波谷值为v
p1
，右边波谷值为v
p2
；v
p1
和v
p2
的计算和m次谐波峰谷比中的算法一致；则m次谐波峰谷差比频率间距dfr的计算公式如下：对多个谐波计算dfr并取均值即为谐波峰谷差比频率间距值；一次谐波峰谷差比频率间距fdfr是m次谐波峰谷差比频率间距的特殊情况，左右波谷值v
p1
和v
p2
计算和一次谐波峰谷比一致；其计算公式为：第1和2次谐波能量比则是一次谐波最大对数幅度与二次谐波最大对数幅度之比；s402、提取时域精细特征，包括单周期基频序列标准差、单周期基频序列均值、单周期基频序列中位数、单周期基频序列均值与中位数距离、单周期基频序列单位时间过均值次数；令步骤s3所求得时域单周期基频序列为f[n
t
]，长度为len，n
t
表示时域单周期的序号；单周期基频序列均值f
mean
的计算公式为：单周期基频序列标准差f
std
的计算公式为：
对f[n
t
]进行排序，则单周期基频序列中位数f
median
的计算公式为：单周期基频序列均值与中位数距离f
mm
的计算公式为：f
mm
＝|f
mean-f
median
|
ꢀꢀꢀ
公式(45)单周期基频序列单位时间过均值次数n
c
的计算公式为：其中t1是时域音频信号的持续时间，u(t)的表达式如下：s403、提取13维梅尔频率倒谱系数静态特征：首先对输入的音频数据进行预处理，其次对音频作fft变换并对fft模值求平方值得到音频的能量谱图，然后利用梅尔滤波器组进行滤波得到每一个滤波器的能量输出，最后对滤波器输出作对数运算并通过离散余弦变换dct得到梅尔频率倒谱系数特征。6.根据权利要求1所述的一种基于时频精细分析的乐器分类方法，其特征在于，所述步骤s5中，进行特征相关性和变异系数评估，建立随机森林模型对乐器进行分类的过程如下：s501、计算特征向量各元素之间的皮尔逊相关系数并评估特征有效性；s502、计算特征向量各元素的变异系数并评估特征有效性；s503、将时频精细特征和梅尔频率倒谱系数特征融合成新特征向量，将融合特征送入随机森林模型中进行训练，训练时采用五折交叉验证；其中，随机森林模型的处理过程如下：利用自助法对训练集进行多次、有放回地随机采样，采样得到的训练样本即为θ；对于每个训练样本θ，形成对应的决策树；在构造的每棵决策树中，从每个训练样本θ的q个特征中，随机选取q个特征进行分类，其中，q≤q，并且表示对向下取整操作；计算选取出的q个特征的基尼不纯度，并认为基尼不纯度最小的特征为最优的分类特征；根据所选的最优分类特征，将最优分类特征所在节点进行二分，从剩余特征中挑选出次优分类特征；重复以上步骤，直到θ剩余特征数量为0或者当前决策树的分类效果最佳；将得到的所有30棵决策树进行组合，得到随机森林模型，并通过投票法得到最终的分类结果。

技术总结
本发明公开了一种基于时频精细分析的乐器分类方法，该方法包括：基频估计、频域谐波标记和时域单周期基频序列估计、时频精细特征提取、随机森林模型。本发明设计基于自相关和窄带谱能量频率估计的循环迭代基频估计算法进行基频估计，然后基于自适应窗口搜索的谐波标记算法和时域单周期基频序列估计算法精确标记谐波位置和获取单周期基频序列；通过时频精细特征，表征频谱精细结构和单周期基频变化并提取；最后将时频精细特征与13维梅尔频率倒谱系数特征融合输入到随机森林模型训练，基于随机森林模型的乐器分类误差率比只用梅尔频率倒谱系数特征降低6％。倒谱系数特征降低6％。倒谱系数特征降低6％。


技术研发人员：王一歌 陈浩文 韦岗 曹燕
受保护的技术使用者：华南理工大学
技术研发日：2023.05.23
技术公布日：2023/9/6