CN103440873A - 一种基于相似性的音乐推荐方法 - Google Patents

Abstract

本发明公开了一种基于混合特征和高斯混合模型的音乐相似度检测方法，基本思路如下：用伽马通倒谱系数进行相似度检测，并将多种特征的加权相似度作为最终的检测结果；提出一种基于帧轴的调制频谱特征，用该特征表示音乐的长时特征，并且将长时特征与短时特征的组合作为下一步建模的输入；使用高斯混合模型对各个音乐特征建模，首先利用动态K均值方法对模型进行初始化，接着用期望最大化算法进行模型训练，获得准确的模型参数，最后使用对数似然比算法获得音乐之间的相似度。本发明对音乐特征的获取更为充分和深入，提高了音乐推荐的准确程度。本发明可降低特征矢量维数，降低音乐数据库信息存储量，提高音乐推荐的准确程度。

Description

一种基于相似性的音乐推荐方法

技术领域

本发明涉及一种互联网的音乐检索技术，特别是一种基于相似性的音乐推荐方法。

背景技术

伴随着网络技术的快速发展，大量音乐在互联网中被分享，在线音乐曲库更新速度很快；传统的音乐检索方式，往往是通过输入歌曲或歌手名字的方式进行，但是在音乐数据如此庞大的信息时代，很难发现歌曲及歌手信息未知的音乐；音乐推荐则使人们能够更加快捷地获得所需音乐。

近年来，许多音乐网站开始提供音乐推荐功能，但是推荐效果尚无法令用户满意，很多推荐仅是对热门曲目的罗列，没有根据用户兴趣进行个性化推荐；针对当前音乐推荐方法尚不完善的现状，迫切地需求一种能深度分析音乐信号，模拟人耳认知的音乐相似性检测算法，以提高音乐推荐的精准度。

音乐推荐是一种根据用户信息分析用户兴趣，产生推荐列表的过程；由于人类对于音乐的认知会受到诸如年龄、背景、情绪等多种因素的影响，对于音乐推荐系统而言，不仅包含了音乐风格的识别，还包括乐曲之间的相似度检测问题；也就是说推荐系统需要模拟人耳的感知特性，根据用户常听曲目，推荐类似音乐给用户。

2006年，在公开号为US20070174274A1的美国专利（Kim H G,Eom K W,Kim J Y,et al.Method and apparatus for searching similar music:U.S.PatentApplication11/487,327.2006）中，申请人公开了一种音乐相似度检测的方法；其基本思路是：(1)提取歌曲的风格特征，称为第一特征，这些特征主要包括改进离散余弦变换（Modified Discrete Cosine Transform，MDCT）域的谱中心、带宽、滚动等信息；(2)获得第一特征的最值（最大值、最小值）、均值和方差，称为第二特征；(3)利用第二特征计算歌曲之间的相似度。该方法的缺点如下：(1)特征提取过程只是简单地获取了频域的几个特征，没有充分地挖掘音乐信号的潜在特点；(2)通过获取最值、均值和方差的方法对特征进行降维，但由于信号的时变性，会导致特征信息的丢失。

2007年，在文献（Signal-based timbre similarity measures for automatic musicrecommendation.USA:The Cooper Union for the Advancement of Science and Art,2007）中，Terence L Magno提出了一种音乐相似度的检测方法；其基本思路是：(1)提取音乐信号的梅尔频率倒谱系数（Mel-Frequency Cepstrum Coefficients，MFCC）特征来表征音乐的音调信息；(2)对特征矢量使用高斯混合模型（Gaussian mixture model，GMM）方法进行建模，获得歌曲的高斯混合模型参数；(3)使用“推土机”距离（Earth Mover’s Distance，EMD）方法，计算歌曲模型之间的距离，进而获得歌曲之间的相似度。该方法的缺点如下：(1)使用梅尔频率倒谱系数来表征音乐信号的音色特征，由于梅尔频率倒谱系数多用于语音信号的处理中，对于背景较为复杂的音乐信号，该特征的应用有一定的限制；(2)“推土机”距离过于依赖模型的精确度，无法反映信号的时变性，从而影响音乐相似性的检测结果。

发明内容

为解决现有技术存在的上述问题，本发明设计了一种能充分挖掘音乐信号的潜在特点和特征，并能反映信号时变性的基于相似性的音乐推荐方法。

为了实现上述目的，本发明提出一种基于混合特征和高斯混合模型的音乐相似度检测方法，基本思路如下：(1)用伽马通倒谱系数（Gamma tone-FrequencyCepstrum Coefficients，GFCC）进行相似度检测，并将多种特征的加权相似度作为最终的检测结果；(2)提出一种基于帧轴的调制频谱特征，用该特征表示音乐的长时特征，并且将长时特征与短时特征的组合作为下一步建模的输入；(3)使用高斯混合模型(Gaussian mixed model,GMM)对各个音乐特征建模，首先利用动态K均值方法对模型进行初始化，接着用期望最大化（Expectation Maximization，EM）算法进行模型训练，获得准确的模型参数，最后使用对数似然比算法获得音乐之间的相似度。

为了检测音乐之间的相似性，需要获得能被计算机识别的反映音乐特点的信息，即特征提取；由于特征的维数一般较高，无法直接计算相似性，因此，需要通过机器学习方法对特征矢量进行建模；最后，通过计算模型之间的相似度，产生推荐列表。

一种基于相似性的音乐推荐方法，包括以下步骤：

A、建立数据库

建立数据库的流程包括特征提取、建立模型和将获得的GMM模型参数保存在模型数据库中，具体步骤如下：

A1、特征提取：以帧为单位提取音乐信号的伽马通倒谱系数特征、情绪特征和八音频谱对比度特征；

A11、预处理

A111、预加重

预处理模块输入为采样率44.1KHz的单声道脉冲编码调制音乐文件，文件中各个数据即是音乐波形的采样值，也是待处理的信号x(n)，这里n表示采样点的序号，定义y(n)为预加重后的输出信号，则

y(n)=x(n)-μx(n-1)

其中，μ称为预加重因子，其取值范围是0<μ<1；

A112、加窗分帧

定义一帧音乐信号的长度为N_f，N_f取值范围是512≤N_f≤8192，帧间交叠的长度为N₀，N₀取值范围是0.25N_f≤N₀≤0.75N_f，对预加重的输出y(n)加汉明窗w(n)=0.54-0.46cos[(2n+1)π/N_f]进行分帧处理，处理后的输出音频信号为s_w(n)，n=0、1、…、N_f-1：

s_w(n)=y(n)·w(n)

A12、提取三类短时特征

A121、提取伽马通倒谱系数

伽马通滤波器是一种基于人耳耳蜗原理的滤波器。伽马通滤波器组中第i个滤波器的时域表示形式为：

g_i(m)＝Bⁿm^n-1exp(-2πBm)cos(2πf_im+φ)·U(m),0≤m＜N_f,0≤i＜M_g    (1)

式中，φ为初始相位，其取值范围是0≤φ<2π，它不影响滤波器组的滤波性能；M_g为滤波器组中的滤波器个数，f_i为伽马通滤波器组中第i个滤波器的中心频率，Bⁿ为第i个滤波器增益，n为第i个滤波器的阶数；B为衰减因子，B的大小与第i个滤波器等效带宽ERB(f_i)有关，即B=1.019ERB(f_i)，第i个滤波器等效带宽ERB(f_i)表示为：

ERB(f_i)＝24.7(4.37f_i1000+1)

对式(1)做快速傅里叶变换，得到伽马通滤波器组中第i个滤波器的频域表达式为：

$H_i\left(k\right)=\left|\underset{m=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{g_i\left(m\right)e}^{(-j2\mathrm{\&pi;km}/N_f)}\right|,0\&le;k<N_f,0\&le;i<M_g---\left(2\right)$

将伽马通滤波器与倒谱系数相结合，得到伽马通倒谱系数。其提取过程为：加窗分帧的音频信号，经快速傅里叶变换后输入至伽马通滤波器组滤波，伽马通滤波器组滤波的输出经对数压缩及离散余弦变换的结果即为所求的伽马通倒谱系数。伽马通倒谱系数的具体提取步骤如下：

A1211、对预处理后的信号s_w(n)进行快速傅里叶变换(FFT)，将时域的音乐信号转换到频域，并计算其能量|S_w(k)|²：

${\left|S_w\left(k\right)\right|}^2={\left|\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{s_w\left(n\right)e}^{-j2\mathrm{\&pi;kn}/N_f}\right|}^2,0\&le;k<N_f---\left(3\right)$

A1212、将式(3)计算出的音乐信号能量通过伽马通滤波器组中第i个滤波器，i=0、1、…、M_g-1，并对滤波器的输出进行对数压缩，得到：

$m_i=\ln \left[\underset{k=0}{\overset{N-1}{\mathrm{\&Sigma;}}}\right|S_w^2\left(k\right)\left|H_i\left(k\right)\right],0\&le;i<M_g,0\&le;k<N_f---\left(4\right)$

A1213、将式(4)表示的对数压缩结果进行离散余弦变换，得到伽马通倒谱系数为：

$\mathrm{GFCC}\left(i\right)=\sqrt{\frac{2}{N}}\underset{g=0}{\overset{M_{\mathrm{gf}}-1}{\mathrm{\&Sigma;}}}{m}_g\cos [\mathrm{\&pi;i}(g+0.5)/M_g],i=\mathrm{0,1},\&CenterDot;\&CenterDot;\&CenterDot;,L-1$

式中，L为伽马通倒谱系数的个数，0<L<M_g；第t帧的伽马通倒谱系数矢量为

$x_t^{\mathrm{GFCC}}={[{\mathrm{GFCC}}_t\left(0\right),{\mathrm{GFCC}}_t\left(1\right),\&CenterDot;\&CenterDot;\&CenterDot;,{\mathrm{GFCC}}_t(L-1)]}^T$

A122、提取情绪特征矢量

音乐中情绪特征包括：时域均值TA、频域均值FA、时域方差TV、频域方差FV、频域最大值FM、时域重心TC、频域重心FC、时域带宽TB、时域滚动TR、频域滚动FR、时域流量TF和频域流量FF；定义x_t(n)表示第t帧音乐信号的离散采样值，n=0、1、…、N_f-1；表示x_t(n)经过快速傅里叶变换后的频域幅值，k=0,1,…,N_f-1，X_t(k)；则时域与频域的各个特征参数的计算公式如下：

A1221、计算时域、频域均值

${\mathrm{TA}}_t=\frac{1}{N_f}\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}{x}_t\left(j\right)$

${\mathrm{FA}}_t=\frac{1}{N_f}\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}{X}_t\left(j\right)$

式中，TA_t为第t帧音乐信号的时域均值，FA_t为第t帧音乐信号的频域均值；

A1222、计算时域、频域方差

${\mathrm{TV}}_t=\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}{(x_t\left(j\right)-{\mathrm{TA}}_t)}^2$

${\mathrm{FV}}_t=\stackrel{N_f-1}{\mathrm{\&Sigma;}}{(X_t\left(j\right)-{\mathrm{FA}}_t)}^2$

式中，TV_t为第t帧音乐信号的时域方差，FV_t为第t帧音乐信号的频域方差；

A1223、计算频域最大值

FM_t＝max{X_t(j)},0≤j＜N_f

式中，FM_t为第t帧音乐信号的频域最大值；

A1224、计算时域、频域重心

$T{C}_t=\frac{\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}j\&CenterDot;x_t\left(j\right)}{\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}{x}_t\left(j\right)}$

${\mathrm{FC}}_t=\frac{\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}j\&CenterDot;X_t\left(j\right)}{\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}{X}_t\left(j\right)}$

式中，TC_t为第t帧音乐信号的时域重心，FC_t为第t帧音乐信号的频域重心；

A1225、计算时域带宽

${\mathrm{TB}}_t=\sqrt{\frac{\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}[{\left|x_t\left(j\right)\right|}^2\&times;(j-{\mathrm{TC}}_t)]}{\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}{\left|x_t\left(j\right)\right|}^2}}$

式中，TB_t为第t帧音乐信号的时域带宽。

A1226、计算时域、频域滚动

$\underset{j=0}{\overset{{\mathrm{TR}}_t}{\mathrm{\&Sigma;}}}{\left|x_t\left(j\right)\right|}^2=0.85\&times;\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}{\left|x_t\left(j\right)\right|}^2$

$\underset{j=0}{\overset{{\mathrm{FR}}_t}{\mathrm{\&Sigma;}}}{\left|X_t\left(j\right)\right|}^2=0.85\&times;\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}{\left|X_t\left(j\right)\right|}^2$

式中，TR_t为第t帧音乐信号的时域滚动，FR_t为第t帧音乐信号的频域滚动；

A1227、计算时域、频域流量

${\mathrm{TF}}_t=\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}{\left[\right|x_t\left(j\right)|-x_{t-1}\left(j\right)]}^2$

${\mathrm{FF}}_t=\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}{\left[\right|X_t\left(j\right)|-X_{t-1}\left(j\right)]}^2$

式中，TF_t为第t帧音乐信号的时域流量，FF_t为第t帧音乐信号的频域流量。

影响情绪模型的时域特征有6种，包括时域方差、时域重心、时域均值、时域带宽、时域滚动和时域流量；频域特征6种，包括频域均值、频域方差、频域重心、频域滚动、频域最大值和频域流量；因此最后获得的情绪特征矢量表示为：

$x_t^{\mathrm{MODD}}={[{\mathrm{TV}}_t,{\mathrm{TC}}_t,{\mathrm{TA}}_t,{\mathrm{TB}}_t,{\mathrm{TR}}_t,{\mathrm{TF}}_t{,\mathrm{FA}}_t,{\mathrm{FV}}_t,{\mathrm{FC}}_t,{\mathrm{FR}}_t,{\mathrm{FM}}_t,{\mathrm{FF}}_t]}^T$

$={[{\mathrm{MOOD}}_t\left(0\right),{\mathrm{MOOD}}_t\left(1\right),\&CenterDot;\&CenterDot;\&CenterDot;,{\mathrm{MOOD}}_t(M-1)]}^T$

式中，M为情绪模型特征矢量维数，其取值范围是3<M<17；

A123、提取八音频谱对比度矢量

帧长为N_f的音乐信号经过快速Fourier变换后，其频域值为

将变换后的结果按频率划分子带，将第b个子频带的向量记为

这里N_b表示该子频带内的点的个数；对于采样率为44.1KHz的音乐信号，每个带通滤波器的频率范围如表1所示：

表1八音频谱滤波器带宽

滤波器编号	物理频率域内带宽范围(Hz)	离散频域内带宽范围(k)
			1	[0,200)	[0,9)
2	[200,400)	[9,19)

3	[400,800)	[19,37)
			4	[800,1600)	[37,74)
5	[1600,3200)	[74,149)
			6	[3200,6400)	[149,297)
7	[6400,12800)	[297,594)
			8	[12800,22050)	[594,1024)

表1中第二列是物理频率域的子频带划分，第三列是离散频域内的子频带划分；将这些点按照降序重新排列，得到{P′_b,1,P′_b,2,…,P′_b,Nb}，式中，P′_b,1≥P′_b,2≥…≥P′_b,Nb，则该子频带的频谱峰值和谷值为：

$\mathrm{Peak}\left(b\right)=10\&CenterDot;{\log }_{10}\left(\frac{1}{{\mathrm{\&alpha;N}}_b}\underset{i=1}{\overset{{\mathrm{\&alpha;N}}_b}{\mathrm{\&Sigma;}}}{P}_{b,i}^{\&prime;}\right)$

$\mathrm{Valley}\left(b\right)=10\&CenterDot;{\log }_{10}\left(\frac{1}{{\mathrm{\&alpha;N}}_b}\underset{i=1}{\overset{{\mathrm{\&alpha;N}}_b}{\mathrm{\&Sigma;}}}{P}_{b,N_b-i+1}^{\&prime;}\right)$

式中，系数α取值范围在0.02～0.2之间；每个子频带内的频谱峰值与谷值之差，称为频谱对比度SC(b)

SC(b)=Peak(b)-Vallay(b)

因此，一帧音乐信号的八音频谱对比度矢量能用各个子频带频谱对比度和子频带频谱谷值来表示，

X_OSC＝[Vallay(0),…,Valley(B-1),SC(0),…,SC(B-1)]

式中，B表示子频带的个数，其取值范围是3<M<33，第t帧的特征矢量为：

$x_t^{\mathrm{OSC}}={[{\mathrm{OSC}}_t\left(0\right),{\mathrm{OSC}}_t\left(1\right),\&CenterDot;\&CenterDot;\&CenterDot;,{\mathrm{OSC}}_t(2B-1)]}^T$

$=[{\mathrm{Vallay}}_t\left(0\right),\&CenterDot;\&CenterDot;\&CenterDot;,{\mathrm{Valley}}_t(B-1),{\mathrm{SC}}_t\left(0\right),\&CenterDot;\&CenterDot;\&CenterDot;{,\mathrm{SC}}_t(B-1)]$

A13、计算长时频谱对比度、长时频谱谷值、长时频谱能量、长时频谱中心和长时频谱平坦度

以帧为单位提取的特征矢量有伽马通倒谱系数、情绪特征矢量和八音频谱对比度矢量，对特征矢量沿帧轴进行离散Fourier变换，得到信号的长时特征；

定义x_p表示伽马通倒谱系数、情绪特征矢量和八音频谱对比度矢量中任意一个特征矢量，x_p=[x_p(0),x_p(1),…,x_p(D-1)]^T表示第p帧提取的特征矢量，D为对应特征矢量的长度。对其沿帧轴进行长度为W的离散傅里叶变换，W取值范围256≤W≤1024，相邻窗之间重叠50%，得到长时频谱M_t(m,d)，m是频域内的索引号，0≤m<W；

$M_t(m,d)=\underset{p=0}{\overset{W-1}{\mathrm{\&Sigma;}}}{x}_{(t\&times;W/2)+p}\left(d\right)e^{-j2\mathrm{\&pi;}\frac{p}{W}m},0\&le;m<W,0\&le;d<D$

将长时频谱M_t(m,d)，0≤m<W，划分为J个调制子频带，J的取值范围6≤J≤32，根据离散傅里叶变换的对称性，W取512时，只需对前256个进行划分即可，因此得到长时频谱子频带范围依次为[0,3)、[3,6)、[6,12)、[12,24)、[24,48)、[48,96)、[48,96)和[192,256)；

在各个子频带内获取长时频谱对比度、长时频谱谷值、长时频谱能量、长时频谱中心和长时频谱平坦度，具体计算如下；

A131、计算长时频谱对比度

长时频谱峰值(MSP)、谷值(MSV)做差即得到长时频谱对比度；

$\mathrm{MSP}(j,d)=\underset{{\mathrm{\&phi;}}_{j,l}\&le;m\&le;{\mathrm{\&phi;}}_{j,h}}{\max }\left[M(m,d)\right]$

$\mathrm{MSV}(j,d)=\underset{{\mathrm{\&phi;}}_{j,l}\&le;m\&le;{\mathrm{\&phi;}}_{j,h}}{\min }\left[M(m,d)\right]$

MSC(j,d)=MSP(j,d)-MSV(j,d)

式中，φ_j,l与φ_j,h对应第j个子频带的低频索引号与高频索引号，与表2对应，式中0≤j<J；

表2调制子频带范围

调制子频带号	频谱索引范围
		1	[0,3)
2	[3,6)
		3	[6,12)
4	[12,24)
		5	[24,48)
6	[48,96)
		7	[96,192)
8	[192,256)

A132、长时频谱能量

长时频谱能量反映了长时频谱频带内的能量分布情况，其计算公式为：

$\mathrm{MSE}(j,d)=10\&CenterDot;{\log }_{10}\left\{(1+\underset{m={\mathrm{\&phi;}}_{j,l}}{\overset{{\mathrm{\&phi;}}_{j,h}}{\mathrm{\&Sigma;}}}{\left[M(m,d)\right]}^2\right\}$

A133、长时频谱中心

长时频谱中心，反映了每个长时频谱子频带的分布情况；

$\mathrm{MSCEN}(j,d)=\frac{\underset{m={\mathrm{\&phi;}}_{j,l}}{\overset{{\mathrm{\&phi;}}_{j,h}}{\mathrm{\&Sigma;}}}M(m,d)\&times;m}{\underset{m={\mathrm{\&phi;}}_{j,l}}{\overset{{\mathrm{\&phi;}}_{j,h}}{\mathrm{\&Sigma;}}}M(m,d)}$

A134、长时频谱平坦度

长时频谱平坦度反映了长时频谱各频带的频谱分布情况，长时频谱平坦度的值较大，说明长时频谱分布较为均匀，反之，说明长时频谱集中分布于某几个频带内；长时频谱的几何均值与代数均值的商定义为

$\mathrm{MSF}(j,d)=\frac{{\mathrm{\&phi;}}_{j,h}-{\mathrm{\&phi;}}_{j,l}+1\sqrt{\underset{m={\mathrm{\&phi;}}_{j,l}}{\overset{{\mathrm{\&phi;}}_{j,h}}{\mathrm{\&Pi;}}}M(m,d)}}{\frac{1}{{\mathrm{\&phi;}}_{j,h}-{\mathrm{\&phi;}}_{j,l}+1}\underset{m={\mathrm{\&phi;}}_{j,l}}{\overset{{\mathrm{\&phi;}}_{j,h}}{\mathrm{\&Sigma;}}}M(m,d)}$

A2、建立模型

用高斯混合模型对调制的伽马通倒谱系数、情绪特征矢量和八音频谱对比度矢量分别进行建模，并使用K均值聚类算法对模型参数进行粗略估计，然后用期望最大化估计算法获得精确的模型参数；最后，将获得的模型参数保存在模型数据库中；

A21、估计高斯混合模型参数

高斯混合模型的数目直接影响建模效果，采取动态K均值的方式调整K均值聚类的个数，进而调整混合模型个数；具体做法是：设置某一阈值，判断聚类的半径是否小于该阈值，若不成立，则增加聚类数目，直到聚类半径小于该阈值为止；具体步骤如下：

A211、动态K均值初始化高斯混合模型

歌曲的总帧数为N_frames，某种特征矢量的维数为D，则每帧的D维特征矢量表示为c_j=[c_j,1,c_j,2,…,c_j,D]，j=0,1,…,N_frames-1；高斯混合模型定义为M个单高斯概率密度函数的加权和，其表达式为：

$p\left(c_j\right)=\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_i{N}_i(c_j;{\mathrm{\&mu;}}_j,{\mathrm{\&Sigma;}}_j)$

式中，i是高斯概率密度函数索引号，πi是相应的权值，满足

N_i(c_j;μ_j,Σ_j)为相应的高斯概率密度函数的表达式，其均值为μ_j＝[μ_j,1,μ_j,2,…,μ_j,N]，方差Σ_j是一个N×N的矩阵，N_i(c_j;μ_j,Σ_j)的表达式为：

$N_i(c_j;{\mathrm{\&mu;}}_j,{\mathrm{\&Sigma;}}_j)=\frac{1}{{\left(2\mathrm{\&pi;}\right)}^{\frac{d}{2}}{\left|{\mathrm{\&Sigma;}}_j\right|}^{\frac{1}{2}}}\exp [-\frac{1}{2}{(c_j-{\mathrm{\&mu;}}_j)}^T{\mathrm{\&Sigma;}}_j^{-1}(c_j-{\mathrm{\&mu;}}_j)]$

则单个高斯概率密度函数的模型Θ_j为：

Θ_j＝{π_j,μ_j,Σ_j},j＝1,2,…,M

高斯混合模型就表示为：

Θ＝{Θ₁,Θ₂,…,Θ_M}

利用一首歌曲的特征向量获得高斯混合模型时，高斯混合模型的个数直接影响建模参数的准确度，进而影响歌曲间的相似度度量结果；采取动态K均值方法先确定聚类的个数，进而确定高斯混合模型的个数；

在聚类时，首先聚类数目初始化为1，所有的特征矢量被聚成一类，聚类半径大于阈值半径；增大聚类数目为2，此时有一个类内的半径大于阈值半径；继续增大聚类数目为3，此时各个聚类的半径均小于阈值半径，则确定3为最终的聚类个数；

聚类半径定义为聚类内的特征到聚类中心的最大距离，即

$R_{\max }=\underset{c_j\&Element;Z_i}{\max }\left(\right|c_j-{\mathrm{\&mu;}}_i\left|\right)---\left(44\right)$

式中，i=1,2,…,K表示聚类后的索引号，Z_i是相应的类，N_p,i是类Z_i内的特征数目，c_j是属于类Z_i中的特征矢量，μ_i是Z_i类的均值矢量，表示为：

${\mathrm{\&mu;}}_i=\frac{1}{N_i}\underset{c_j\&Element;Z_i}{\mathrm{\&Sigma;}}{c}_j---\left(45\right)$

应用动态K均值方法确定聚类数目后，需要计算各个聚类的均值、方差和权值来初始化高斯混合模型；均值由式(49)求得；权值用聚类内的特征向量个数除以特征矢量总数来计算，其表达式为：

${\mathrm{\&pi;}}_i=\frac{N_{p,i}}{N_{\mathrm{frames}}}---\left(46\right)$

为了获得各个聚类的方差，为每个聚类构建N_p,i×N的矩阵，即：

C_i＝[c₁ c₂ … c_p,i]^T    (47)

式中，

是矩阵的第k列，是特征矢量的第k维，则方差矩阵的第m行n列元素的求法为：

${\mathrm{\&Sigma;}}_i^{(m,n)}=\mathrm{cov}(C_i^{(:,m)},C_i^{(:,n)})$

                            (48)

$=E\left[(C_i^{(:,m)}-M_m)\left((C_i^{(:,n)}-M_n\right)\right]$

式中，m,n=1,2,…,N，M_k是每个元素均为m_k的向量，式中m_k是聚类Z_i中第k维特征的均值，表示为：

$m_k=\frac{1}{N_{p,i}}\underset{j=1}{\overset{N_{p,i}}{\mathrm{\&Sigma;}}}{C}_i^{(j,k)}---\left(49\right)$

A22、精确估计高斯混合模型参数

高斯混合模型参数估计就是依据某种参数估计准则确定模型参数的过程，用极大似然准则进行估计；一首歌曲的特征矩阵表示为这里N_frames表示帧的总数目，

是每帧的特征矢量，假设它们彼此独立，则高斯混合模型的似然函数定义为，在数理统计学中，似然函数是一种关于统计模型中的参数的函数,表示模型参数中的似然性：

$L(\mathrm{\&Theta;}/C)=p(C/\mathrm{\&Theta;})=\underset{i=1}{\overset{N_{\mathrm{frames}}}{\mathrm{\&Pi;}}}p(c_i/\mathrm{\&Theta;})---\left(50\right)$

将似然函数L(Θ/C)看成Θ的非线性函数，参数估计的目的是要找到使似然函数L(Θ/C)最大的模型Θ，即：

${\mathrm{\&Theta;}}^{*}=\arg \underset{\mathrm{\&Theta;}}{\max }p(C/\mathrm{\&Theta;})---\left(51\right)$

对式(51)取对数，获得对数似然函数

$J\left(\mathrm{\&Theta;}\right)=\log L(\mathrm{\&Theta;}/C)=\log \underset{i=1}{\overset{N_{\mathrm{frames}}}{\mathrm{\&Pi;}}}p(c_i/\mathrm{\&Theta;})---\left(52\right)$

实际中常用期望最大化算法来估计模型参数，于是得到模型权值为：

${\mathrm{\&pi;}}_k=\frac{1}{N_{\mathrm{frames}}}\underset{i=1}{\overset{N_{\mathrm{frames}}}{\mathrm{\&Sigma;}}}p\left(k\right|c_i,{\mathrm{\&Theta;}}^{\mathrm{old}})---\left(53\right)$

均值为：

${\mathrm{\&mu;}}_k=\frac{\underset{i=1}{\overset{N_{\mathrm{frames}}}{\mathrm{\&Sigma;}}}{c}_{i}p\left(k\right|c_i,{\mathrm{\&Theta;}}^{\mathrm{old}})}{\underset{i=1}{\overset{N_{\mathrm{frames}}}{\mathrm{\&Sigma;}}}p\left(k\right|c_i,{\mathrm{\&Theta;}}^{\mathrm{old}})}---\left(54\right)$

方差矩阵为：

${\mathrm{\&Sigma;}}_k=\frac{\underset{i=1}{\overset{N_{\mathrm{frames}}}{\mathrm{\&Sigma;}}}p\left(k\right|c_i,{\mathrm{\&Theta;}}^{\mathrm{old}})(c_i-{\mathrm{\&mu;}}_k){(c_i-{\mathrm{\&mu;}}_k)}^T}{\underset{i=1}{\overset{N_{\mathrm{frames}}}{\mathrm{\&Sigma;}}}p\left(k\right|c,{\mathrm{\&Theta;}}^{\mathrm{old}})}---\left(55\right)$

式中，p(k|c_i,Θ^old)称为后验概率；根据贝叶斯准则，该后验概率为：

$p\left(k\right|c_i,{\mathrm{\&Theta;}}^{\mathrm{old}})=\frac{p(k,c_i|{\mathrm{\&Theta;}}^{\mathrm{old}})}{p\left(c_i\right|{\mathrm{\&Theta;}}^{\mathrm{old}})}=\frac{{\mathrm{\&pi;}}_k^{\mathrm{old}}N(c_i;{\mathrm{\&mu;}}_k^{\mathrm{old}},{\mathrm{\&Sigma;}}_k^{\mathrm{old}})}{\underset{l=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_l^{\mathrm{old}}N(c_i;{\mathrm{\&mu;}}_l^{\mathrm{old}},{\mathrm{\&Sigma;}}_l^{\mathrm{old}})}---\left(56\right)$

期望最大化估计的步骤为：

A221、根据式(41)和(56)得到后验概率；

A222、根据式(53)～(55)计算新的模型参数，以此更新上次的模型参数，当新旧模型参数几乎一致时，迭代停止，否则返回(a)步骤继续迭代计算；

B、生成推荐列表

生成推荐列表的流程包括提取输入歌曲的特征矢量、分别计算该歌曲与数据库中歌曲的加权相似度和对获得的相似度按降序排列产生推荐列表，具体步骤如下：

B1、提取音乐信号的特征矢量，提取步骤与步骤A1所述的步骤完全一致，只是输入音乐信号不同，生成推荐列表的输入音乐信号为用户所试听的一支曲目，经步骤A1处理的输出即为产生类似该曲目风格的音乐，而步骤A1中的输入音乐信号由建立数据库的许多音乐信号组成；

B2、分别计算该歌曲与数据库中歌曲的加权相似度

得到特征矢量的统计模型后，通过模型之间的相似度来判定歌曲的相似度；若M₁维的高斯混合模型模型A为：

${\mathrm{\&Theta;}}_A=\{({\mathrm{\&mu;}}_{1,A},{\mathrm{\&Sigma;}}_{1,A},{\mathrm{\&pi;}}_{1,A}),\&CenterDot;\&CenterDot;\&CenterDot;,({\mathrm{\&mu;}}_{M,A},{\mathrm{\&Sigma;}}_{M,A},{\mathrm{\&pi;}}_{M_1,A})\}---\left(57\right)$

对应的特征矩阵为：

$C_A={[c_{1,A},c_{2,A},\&CenterDot;\&CenterDot;\&CenterDot;,c_{N_{\mathrm{Aframes}},A}]}^T---\left(58\right)$

式中，N_Aframes为模型A特征矢量的个数；

另一M₂维的高斯混合模型模型B为：

${\mathrm{\&Theta;}}_B=\{({\mathrm{\&mu;}}_{1,B},{\mathrm{\&Sigma;}}_{1,B},{\mathrm{\&pi;}}_{1,B}),\&CenterDot;\&CenterDot;\&CenterDot;,({\mathrm{\&mu;}}_{M,B},{\mathrm{\&Sigma;}}_{M,B},{\mathrm{\&pi;}}_{M_2,B})\}---\left(59\right)$

对应的特征矢量为：

$C_B={[c_{1,B},c_{2,B},\&CenterDot;\&CenterDot;\&CenterDot;,c_{N_{\mathrm{Bframes}},B}]}^T---\left(60\right)$

式中，N_Bframes为模型B特征矢量的个数；

相似性的计算公式由下式获得：

r(A,B)＝logL(Θ_A/C_A)+logL(Θ_B/C_B)-logL(Θ_A/C_B)-logL(Θ_B/C_A)    (61)

式(61)的最大值为：

r_max(A,B)＝logL(Θ_A/C_A)+logL(Θ_B/C_B)    (62)

则模型A与模型B的相似度定义为：

$\mathrm{sim}(A,B)=\frac{r_{\max }(A,B)-r(A,B)}{r_{\max }(A,B)}=\frac{\log L({\mathrm{\&Theta;}}_A/C_B)+\log L({\mathrm{\&Theta;}}_B/C_A)}{\log L({\mathrm{\&Theta;}}_A/C_A)+\log L({\mathrm{\&Theta;}}_B/C_B)}---\left(63\right)$

为了计算歌曲m与歌曲n的相似度，利用式(63)分别计算上述三类特征的相似度，即sim_gfcc(m,n)、sim_osc(m,n)和sim_mood(m,n)，则两首歌曲的总相似度表示为：

sim(m,n)＝w_gfcc·sim_gfcc(m,n)+w_osc·sim_osc(m,n)+w_mood·sim_mood(m,n)    (64)

式中，w_gfcc、w_osc与w_mood为三类特征相似度的加权值，取值范围是0<w_gfcc<1、0<w_mood<1、0<w_osc<1，且w_gfcc+w_mood+w_osc=1；

B3、对获得的相似度降序排列，产生推荐列表。

本发明所述的预加重因子μ的最佳值为0.97；所述的一帧音乐信号的长度N_f的最佳值为2048；所述的帧间交叠的长度的最佳值为N₀=0.5N_f=1024；所述的初始相位φ最佳值为0；所述的伽马通倒谱系数的个数最佳值为L=26；所述的情绪模型特征矢量维数最佳值为M=12；所述的子频带的个数B的最佳值为8；所述的W最佳值为512帧；所述的三类特征相似度的加权值的最佳值分别为w_gfcc=0.6、w_mood=0.2和w_osc=0.2；所述的系数α的最佳值为0.2。

与现有技术相比，本发明具有以下有益效果：

1、本发明用Gamma tone倒谱系数代替传统的梅尔倒谱系数，进行音色特征的提取，对音乐特征的获取更为充分和深入，提高了音乐推荐的准确程度。

2、本发明通过时间轴的调制技术，可降低特征矢量维数，降低音乐数据库信息存储量；长时特征与短时特征结合，充分获取了音乐的静态与动态特征，较多地保留了音乐信号特征，提高了音乐推荐的准确程度。

3、本发明使用EMD距离计算相似度，与音乐信号的动态特征无关，保证了推荐结果的稳定性。客观测试结果显示，本发明的推荐准确度在86%以上。主观测试结果表明，本发明与人的主观感知近似。与其他算法相比，本发明算法的推荐效果好于现有算法。

附图说明

本发明共有附图6张，其中：

图1是音乐模型数据库建立流程图。

图2是音乐推荐列表生成流程图。

图3是伽马通倒谱系数提取过程。

图4是动态K均值聚类数目为1的示意图。

图5是动态K均值聚类数目为2的示意图。

图6是动态K均值聚类数目为3的示意图。

具体实施方式

下面结合附图对本发明进行进一步地描述。图1所示为本发明步骤A的流程图，图2所示为本发明步骤B的流程图，图3所示为本发明步骤A121的流程图，图4-6所示为本发明步骤A211的动态K均值方法的聚类过程示意图。

为验证本发明技术的有效性，按照图1-6所示的流程图，进行了如下的客观测试，具体步骤如下：

1、建立数据库

建立5个测试数据库，每个曲库中含有200首歌曲，每首歌曲为44.1kHz的PCM文件；其中每个曲库中曲目组成如下：(1)20首歌曲作为测试曲目，此处称为“种子歌曲”，它们来自同一歌手；(2)80首歌曲与“种子歌曲”属于同一种风格；(3)100首歌曲与“种子歌曲”既不是同一歌手，也不是同一风格；

每个测试数据库的“种子歌曲”风格不同，以验证方法的推荐性能；从曲库的组成可以看出，对于任意一个种子歌曲，在曲库内，有另外百分之十的曲目跟它很相似(同一歌手)，百分之四十的曲目为基本相似(同种风格)，百分之五十的曲目不相似(不同歌手、不同风格)；

2、确定评价指标

(1)客观评价标准

为了衡量推荐结果的好坏，进而评价本发明方法的优劣，使用推荐准确度作为客观评价指标；对于任意一个“种子歌曲”，将推荐列表中同一歌手和同一风格歌曲数目所占的比例定义为客观推荐准确度，即：

由此可知，对于该曲库的整体推荐精确度，为多次更换“种子歌曲”后推荐准确度的均值，即：

其中，推荐列表歌曲数目一般设置为5或10，当列表数目过大时，会导致列表中排序靠后的曲目与“种子歌曲”的相似度较低，影响推荐指标的评价；对于某一次产生10列表的推荐中，如果推荐列表中有5首歌曲来自同一歌手、2首歌曲来自同一风格、3首歌曲既不是同一歌手也不是同一风格，则根据上述定义，此时的客观推荐准确度为70%；

歌曲的相似与否是人耳的一种主观感受，为了衡量本发明方法的相似度与人耳真实感受的差距，引入主观评价指标；为了获得人的主观相似度，针对曲库中的“种子歌曲”，评价另外199首歌曲与该歌曲的相似程度，并且进行评分，评分标准如表3所示；

表3主观相似度评分标准

分数	含义	描述
			4	很相似	如果一个人喜欢其中的一首歌曲，则一定会喜欢另一首歌
3	相似	如果一个人喜欢其中的一首歌曲，很可能也喜欢另一首歌
			2	不相似	一个人是否喜欢其中的一首歌曲，不影响其喜欢另一首歌
1	完全不同	一个人不可能在同一时间，同时喜欢这两首歌

表3给出了几种主观感受对应的打分值，打分的范围限制在0～5分之间，分数不限整数，可给予该范围内的小数分值，如2.6等；若听众无法确定分数时，此处分值随机生成；评分时，听众完全看不到音乐的任何信息(音乐名字、风格等)，分数完全凭借自己的主观感受给出；这样，对于每一个种子歌曲，可以获得一个1×200的主观评分矩阵；

歌曲的相似性不仅是一种主观的听觉感受，还因人而异，因此，引入相似系数来度量不同人之间打分矩阵的相近程度；若两个不同人的打分矩阵分别为X与Y，则二者的相似系数为：

${\mathrm{\&rho;}}_{\mathrm{XY}}=\frac{\mathrm{cov}(X,Y)}{{\mathrm{\&sigma;}}_X{\mathrm{\&sigma;}}_Y}---\left(67\right)$

其中，cov(X,Y)是二者的协方差，σ_X与σ_Y是打分矩阵的方差；研究表明，不同歌曲之间的相似度矩阵的相关系数上限为0.613，Cohen对于相关系数的解释见表4；因此，可以通过计算相似度矩阵与人主观相似度矩阵的相关系数的方法，来评价本发明技术的性能；

表4相关系数含义解释

相关性	负相关	正相关
			小	[-0.29,-0.10]	[0.10,0.29]
中	[-0.49,-0.30]	[0.30,0.49]
			大	[-1.00,-0.50]	[0.50,1.00]

2.3.3测试结果

表5给出了本发明与已有的几种推荐系统的结果比较；从表5可以看出，本发明的推荐精度高于其他已有算法；这是因为本发明在特征提取处，充分融合了多种特征；在建模时，动态初始化模型参数，使得模型更加精准；在最后的相似性计算中，结合对数似然原理，获得更加接近人耳的相似性结果；

表5本发明方法与已有其它方法的性能比较

方法	准确度	相似度系数
			动态K均值方法	80%	无
Magno方法	无	0.547
			SSPK2方法	77.05%	无
本发明方法	86%	0.59

表5中，动态K均值方法来自于D.M.Kim,K.S.Kim,K.H.Park.著《Amusicrecommendation system with a dynamic k-means clustering algorithm[C].International Conference on Machine Learning and Applications》Cincinnati,OH,USA,2007:399-403。Magno方法来自于T.L.Magno.著《Signal-based timbresimilarity measures for automatic music recommendation[D].USA:The CooperUnion for the Advancement of Science and Art》2007。SSPK2方法来自于陈捷,许洁萍,刘璇著《基于内容的音乐相似计算研究》，第七届和谐人机环境联合学术会议论文集[C],北京,2011。

Claims (2)

1.一种基于相似性的音乐推荐方法，其特征在于：包括以下步骤：

A、建立数据库

建立数据库的流程包括特征提取、建立模型和将获得的GMM模型参数保存在模型数据库中，具体步骤如下：

A1、特征提取：以帧为单位提取音乐信号的伽马通倒谱系数特征、情绪特征和八音频谱对比度特征；

A11、预处理

A111、预加重

预处理模块输入为采样率44.1KHz的单声道脉冲编码调制音乐文件，文件中各个数据即是音乐波形的采样值，也是待处理的信号x(n)，这里n表示采样点的序号，定义y(n)为预加重后的输出信号，则

y(n)=x(n)-μx(n-1)

其中，μ称为预加重因子，其取值范围是0<μ<1；

A112、加窗分帧

定义一帧音乐信号的长度为N_f，N_f取值范围是512≤N_f≤8192，帧间交叠的长度为N₀，N₀取值范围是0.25N_f≤N₀≤0.75N_f，对预加重的输出y(n)加汉明窗w(n)=0.54-0.46cos[(2n+1)π/N_f]进行分帧处理，处理后的输出音频信号为s_w(n)，n=0、1、…、N_f-1：

s_w(n)=y(n)·w(n)

A12、提取三类短时特征

A121、提取伽马通倒谱系数

伽马通滤波器是一种基于人耳耳蜗原理的滤波器；伽马通滤波器组中第i个滤波器的时域表示形式为：

g_i(m)＝Bⁿm^n-1exp(-2πBm)cos(2πf_im+φ)·U(m),0≤m＜N_f,0≤i＜M_g    (1)

式中，φ为初始相位，其取值范围是0≤φ<2π，它不影响滤波器组的滤波性能；M_g为滤波器组中的滤波器个数，f_i为伽马通滤波器组中第i个滤波器的中心频率，Bⁿ为第i个滤波器增益，n为第i个滤波器的阶数；B为衰减因子，B的大小与第i个滤波器等效带宽ERB(f_i)有关，即B=1.019ERB(f_i)，第i个滤波器等效带宽ERB(f_i)表示为：

ERB(f_i)＝24.7(4.37f_i/1000+1)

对式(1)做快速傅里叶变换，得到伽马通滤波器组中第i个滤波器的频域表达式为：

$H_i\left(k\right)=\left|\underset{m=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{g_i\left(m\right)e}^{(-j2\mathrm{\&pi;km}/N_f)}\right|,0\&le;k<N_f,0\&le;i<M_g---\left(2\right)$

将伽马通滤波器与倒谱系数相结合，得到伽马通倒谱系数；其提取过程为：加窗分帧的音频信号，经快速傅里叶变换后输入至伽马通滤波器组滤波，伽马通滤波器组滤波的输出经对数压缩及离散余弦变换的结果即为所求的伽马通倒谱系数；伽马通倒谱系数的具体提取步骤如下：

A1211、对预处理后的信号s_w(n)进行快速傅里叶变换(FFT)，将时域的音乐信号转换到频域，并计算其能量|S_w(k)|²：

${\left|S_w\left(k\right)\right|}^2={\left|\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}{s_w\left(n\right)e}^{-j2\mathrm{\&pi;kn}/N_f}\right|}^2,0\&le;k<N_f---\left(3\right)$

A1212、将式(3)计算出的音乐信号能量通过伽马通滤波器组中第i个滤波器，i=0、1、…、M_g-1，并对滤波器的输出进行对数压缩，得到：

$m_i=\ln \left[\underset{k=0}{\overset{N-1}{\mathrm{\&Sigma;}}}\right|S_w^2\left(k\right)\left|H_i\left(k\right)\right],0\&le;i<M_g,0\&le;k<N_f---\left(4\right)$

A1213、将式(4)表示的对数压缩结果进行离散余弦变换，得到伽马通倒谱系数为：

$\mathrm{GFCC}\left(i\right)=\sqrt{\frac{2}{N}}\underset{g=0}{\overset{M_{\mathrm{gf}}-1}{\mathrm{\&Sigma;}}}{m}_g\cos [\mathrm{\&pi;i}(g+0.5)/M_g],i=\mathrm{0,1},\&CenterDot;\&CenterDot;\&CenterDot;,L-1$

式中，L为伽马通倒谱系数的个数，0<L<M_g；第t帧的伽马通倒谱系数矢量为

$x_t^{\mathrm{GFCC}}={[{\mathrm{GFCC}}_t\left(0\right),{\mathrm{GFCC}}_t\left(1\right),\&CenterDot;\&CenterDot;\&CenterDot;,{\mathrm{GFCC}}_t(L-1)]}^T$

A122、提取情绪特征矢量

音乐中情绪特征包括：时域均值TA、频域均值FA、时域方差TV、频域方差FV、频域最大值FM、时域重心TC、频域重心FC、时域带宽TB、时域滚动TR、频域滚动FR、时域流量TF和频域流量FF；定义x_t(n)表示第t帧音乐信号的离散采样值，n=0、1、…、N_f-1；表示x_t(n)经过快速傅里叶变换后的频域幅值，k=0,1,…,N_f-1，X_t(k)；则时域与频域的各个特征参数的计算公式如下：

A1221、计算时域、频域均值

${\mathrm{TA}}_t=\frac{1}{N_f}\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}{x}_t\left(j\right)$

${\mathrm{FA}}_t=\frac{1}{N_f}\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}{X}_t\left(j\right)$

式中，TA_t为第t帧音乐信号的时域均值，FA_t为第t帧音乐信号的频域均值；

A1222、计算时域、频域方差

${\mathrm{TV}}_t=\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}{(x_t\left(j\right)-{\mathrm{TA}}_t)}^2$

${\mathrm{FV}}_t=\stackrel{N_f-1}{\mathrm{\&Sigma;}}{(X_t\left(j\right)-{\mathrm{FA}}_t)}^2$

式中，TV_t为第t帧音乐信号的时域方差，FV_t为第t帧音乐信号的频域方差；

A1223、计算频域最大值

FM_t＝max{X_t(j)},0≤j＜N_f

式中，FM_t为第t帧音乐信号的频域最大值；

A1224、计算时域、频域重心

${\mathrm{TC}}_t=\frac{\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}{x}_t\left(j\right)\&times;j}{\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}{x}_t\left(j\right)}$

${\mathrm{FC}}_t=\frac{\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}j\&CenterDot;X_t\left(j\right)}{\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}{X}_t\left(j\right)}$

式中，TC_t为第t帧音乐信号的时域重心，FC_t为第t帧音乐信号的频域重心；

A1225、计算时域带宽

${\mathrm{TB}}_t=\sqrt{\frac{\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}[{\left|x_t\left(j\right)\right|}^2\&times;(j-{\mathrm{TC}}_t)]}{\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}{\left|x_t\left(j\right)\right|}^2}}$

式中，TB_t为第t帧音乐信号的时域带宽；

A1226、计算时域、频域滚动

$\underset{j=0}{\overset{{\mathrm{TR}}_t}{\mathrm{\&Sigma;}}}{\left|x_t\left(j\right)\right|}^2=0.85\&times;\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}{\left|x_t\left(j\right)\right|}^2$

$\underset{j=0}{\overset{{\mathrm{FR}}_t}{\mathrm{\&Sigma;}}}{\left|X_t\left(j\right)\right|}^2=0.85\&times;\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}{\left|X_t\left(j\right)\right|}^2$

式中，TR_t为第t帧音乐信号的时域滚动，FR_t为第t帧音乐信号的频域滚动；

A1227、计算时域、频域流量

${\mathrm{TF}}_t=\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}{\left[\right|x_t\left(j\right)|-x_{t-1}\left(j\right)]}^2$

${\mathrm{FF}}_t=\underset{j=0}{\overset{N_f-1}{\mathrm{\&Sigma;}}}{\left[\right|X_t\left(j\right)|-X_{t-1}\left(j\right)]}^2$

式中，TF_t为第t帧音乐信号的时域流量、FF_t为第t帧音乐信号的频域流量；

影响情绪模型的时域特征有6种，包括时域方差、时域重心、时域均值、时域带宽、时域滚动和时域流量；频域特征6种，包括频域均值、频域方差、频域重心、频域滚动、频域最大值和频域流量；因此最后获得的情绪特征矢量表示为：

$x_t^{\mathrm{MODD}}={[{\mathrm{TV}}_t,{\mathrm{TC}}_t,{\mathrm{TA}}_t,{\mathrm{TB}}_t,{\mathrm{TR}}_t,{\mathrm{TF}}_t{,\mathrm{FA}}_t,{\mathrm{FV}}_t,{\mathrm{FC}}_t,{\mathrm{FR}}_t,{\mathrm{FM}}_t,{\mathrm{FF}}_t]}^T$

$={[{\mathrm{MOOD}}_t\left(0\right),{\mathrm{MOOD}}_t\left(1\right),\&CenterDot;\&CenterDot;\&CenterDot;,{\mathrm{MOOD}}_t(M-1)]}^T$

式中，M为情绪模型特征矢量维数，其取值范围是3<M<17；

A123、提取八音频谱对比度矢量

帧长为N_f的音乐信号经过快速Fourier变换后，其频域值为

将变换后的结果按频率划分子带，将第b个子频带的向量记为

这里N_b表示该子频带内的点的个数；对于采样率为44.1KHz的音乐信号，每个带通滤波器的频率范围如表1所示：

表1八音频谱滤波器带宽

滤波器编号 物理频率域内带宽范围(Hz) 离散频域内带宽范围(k) 1 [0,200) [0,9) 2 [200,400) [9,19) 3 [400,800) [19,37) 4 [800,1600) [37,74) 5 [1600,3200) [74,149) 6 [3200,6400) [149,297) 7 [6400,12800) [297,594) 8 [12800,22050) [594,1024)

表1中第二列是物理频率域的子频带划分，第三列是离散频域内的子频带划分；将这些点按照降序重新排列，得到{P′_b,1,P′_b,2,…,P′_b,Nb}，式中，P′_b,1≥P′_b,2≥…≥P′_b,Nb，则该子频带的频谱峰值和谷值为：

$\mathrm{Peak}\left(b\right)=10\&CenterDot;{\log }_{10}\left(\frac{1}{{\mathrm{\&alpha;N}}_b}\underset{i=1}{\overset{{\mathrm{\&alpha;N}}_b}{\mathrm{\&Sigma;}}}{P}_{b,i}^{\&prime;}\right)$

$\mathrm{Valley}\left(b\right)=10\&CenterDot;{\log }_{10}\left(\frac{1}{{\mathrm{\&alpha;N}}_b}\underset{i=1}{\overset{{\mathrm{\&alpha;N}}_b}{\mathrm{\&Sigma;}}}{P}_{b,N_b-i+1}^{\&prime;}\right)$

式中，系数α取值范围在0.02～0.2之间；每个子频带内的频谱峰值与谷值之差，称为频谱对比度SC(b)

SC(b)=Peak(b)-Vallay(b)

因此，一帧音乐信号的八音频谱对比度矢量能用各个子频带频谱对比度和子频带频谱谷值来表示，

X_OSC＝[Vallay(0),…,Valley(B-1),SC(0),…,SC(B-1)]

式中，B表示子频带的个数，其取值范围是3<M<33，第t帧的特征矢量为：

$x_t^{\mathrm{OSC}}={[{\mathrm{OSC}}_t\left(0\right),{\mathrm{OSC}}_t\left(1\right),\&CenterDot;\&CenterDot;\&CenterDot;,{\mathrm{OSC}}_t(2B-1)]}^T$

$=[{\mathrm{Vallay}}_t\left(0\right),\&CenterDot;\&CenterDot;\&CenterDot;,{\mathrm{Valley}}_t(B-1),{\mathrm{SC}}_t\left(0\right),\&CenterDot;\&CenterDot;\&CenterDot;{,\mathrm{SC}}_t(B-1)]$

A13、计算长时频谱对比度、长时频谱谷值、长时频谱能量、长时频谱中心和长时频谱平坦度

以帧为单位提取的特征矢量有伽马通倒谱系数、情绪特征矢量和八音频谱对比度矢量，对特征矢量沿帧轴进行离散Fourier变换，得到信号的长时特征；

定义x_p表示伽马通倒谱系数、情绪特征矢量和八音频谱对比度矢量中任意一个特征矢量，x_p=[x_p(0),x_p(1),…,x_p(D-1)]^T表示第p帧提取的特征矢量，D为对应特征矢量的长度；对其沿帧轴进行长度为W的离散傅里叶变换，W取值范围256≤W≤1024，相邻窗之间重叠50%，得到长时频谱M_t(m,d)，m是频域内的索引号，0≤m<W；

$M_t(m,d)=\underset{p=0}{\overset{W-1}{\mathrm{\&Sigma;}}}{x}_{(t\&times;W/2)+p}\left(d\right)e^{-j2\mathrm{\&pi;}\frac{p}{W}m},0\&le;m<W,0\&le;d<D$

将长时频谱M_t(m,d)，0≤m<W，划分为J个调制子频带，J的取值范围6≤J≤32，根据离散傅里叶变换的对称性，W取512时，只需对前256个进行划分即可，因此得到长时频谱子频带范围依次为[0,3)、[3,6)、[6,12)、[12,24)、[24,48)、[48,96)、[48,96)和[192,256)；

在各个子频带内获取长时频谱对比度、长时频谱谷值、长时频谱能量、长时频谱中心和长时频谱平坦度，具体计算如下；

A131、计算长时频谱对比度

长时频谱峰值(MSP)、谷值(MSV)做差即得到长时频谱对比度；

$\mathrm{MSP}(j,d)=\underset{{\mathrm{\&phi;}}_{j,l}\&le;m\&le;{\mathrm{\&phi;}}_{j,h}}{\max }\left[M(m,d)\right]$

$\mathrm{MSV}(j,d)=\underset{{\mathrm{\&phi;}}_{j,l}\&le;m\&le;{\mathrm{\&phi;}}_{j,h}}{\min }\left[M(m,d)\right]$

MSC(j,d)=MSP(j,d)-MSV(j,d)

式中，φ_j,l与φ_j,h对应第j个子频带的低频索引号与高频索引号，与表2对应，式中0≤j<J；

表2调制子频带范围

调制子频带号 频谱索引范围 1 [0,3) 2 [3,6) 3 [6,12) 4 [12,24) 5 [24,48) 6 [48,96) 7 [96,192) 8 [192,256)

A132、长时频谱能量

长时频谱能量反映了长时频谱频带内的能量分布情况，其计算公式为：

$\mathrm{MSE}(j,d)=10\&CenterDot;{\log }_{10}\left\{(1+\underset{m={\mathrm{\&phi;}}_{j,l}}{\overset{{\mathrm{\&phi;}}_{j,h}}{\mathrm{\&Sigma;}}}{\left[M(m,d)\right]}^2\right\}$

A133、长时频谱中心

长时频谱中心，反映了每个长时频谱子频带的分布情况；

$\mathrm{MSCEN}(j,d)=\frac{\underset{m={\mathrm{\&phi;}}_{j,l}}{\overset{{\mathrm{\&phi;}}_{j,h}}{\mathrm{\&Sigma;}}}M(m,d)\&times;m}{\underset{m={\mathrm{\&phi;}}_{j,l}}{\overset{{\mathrm{\&phi;}}_{j,h}}{\mathrm{\&Sigma;}}}M(m,d)}$

A134、长时频谱平坦度

长时频谱平坦度反映了长时频谱各频带的频谱分布情况，长时频谱平坦度的值较大，说明长时频谱分布较为均匀，反之，说明长时频谱集中分布于某几个频带内；长时频谱的几何均值与代数均值的商定义为

$\mathrm{MSF}(j,d)=\frac{{\mathrm{\&phi;}}_{j,h}-{\mathrm{\&phi;}}_{j,l}+1\sqrt{\underset{m={\mathrm{\&phi;}}_{j,l}}{\overset{{\mathrm{\&phi;}}_{j,h}}{\mathrm{\&Pi;}}}M(m,d)}}{\frac{1}{{\mathrm{\&phi;}}_{j,h}-{\mathrm{\&phi;}}_{j,l}+1}\underset{m={\mathrm{\&phi;}}_{j,l}}{\overset{{\mathrm{\&phi;}}_{j,h}}{\mathrm{\&Sigma;}}}M(m,d)}$

A2、建立模型

用高斯混合模型对调制的伽马通倒谱系数、情绪特征矢量和八音频谱对比度矢量分别进行建模，并使用K均值聚类算法对模型参数进行粗略估计，然后用期望最大化估计算法获得精确的模型参数；最后，将获得的模型参数保存在模型数据库中；

A21、估计高斯混合模型参数

高斯混合模型的数目直接影响建模效果，采取动态K均值的方式调整K均值聚类的个数，进而调整混合模型个数；具体做法是：设置某一阈值，判断聚类的半径是否小于该阈值，若不成立，则增加聚类数目，直到聚类半径小于该阈值为止；具体步骤如下：

A211、动态K均值初始化高斯混合模型

歌曲的总帧数为N_frames，某种特征矢量的维数为D，则每帧的D维特征矢量表示为c_j=[c_j,1,c_j,2,…,c_j,D]，j=0,1,…,N_frames-1；高斯混合模型定义为M个单高斯概率密度函数的加权和，其表达式为：

$p\left(c_j\right)=\underset{i=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_i{N}_i(c_j;{\mathrm{\&mu;}}_j,{\mathrm{\&Sigma;}}_j)$

式中，i是高斯概率密度函数索引号，π_i是相应的权值，满足N_i(c_j;μ_j,Σ_j)为相应的高斯概率密度函数的表达式，其均值为μ_j＝[μ_j,1,μ_j,2,…,μ_j,N]，方差Σ_j是一个N×N的矩阵，N_i(c_j;μ_j,Σ_j)的表达式为：

$N_i(c_j;{\mathrm{\&mu;}}_j,{\mathrm{\&Sigma;}}_j)=\frac{1}{{\left(2\mathrm{\&pi;}\right)}^{\frac{d}{2}}{\left|{\mathrm{\&Sigma;}}_j\right|}^{\frac{1}{2}}}\exp [-\frac{1}{2}{(c_j-{\mathrm{\&mu;}}_j)}^T{\mathrm{\&Sigma;}}_j^{-1}(c_j-{\mathrm{\&mu;}}_j)]$

则单个高斯概率密度函数的模型Θ_j为：

Θ_j＝{π_j,μ_j,Σ_j},j＝1,2,…,M

高斯混合模型就表示为：

Θ＝{Θ₁,Θ₂,…,Θ_M}

利用一首歌曲的特征向量获得高斯混合模型时，高斯混合模型的个数直接影响建模参数的准确度，进而影响歌曲间的相似度度量结果；采取动态K均值方法先确定聚类的个数，进而确定高斯混合模型的个数；

在聚类时，首先聚类数目初始化为1，所有的特征矢量被聚成一类，聚类半径大于阈值半径；增大聚类数目为2，此时有一个类内的半径大于阈值半径；继续增大聚类数目为3，此时各个聚类的半径均小于阈值半径，则确定3为最终的聚类个数；

聚类半径定义为聚类内的特征到聚类中心的最大距离，即

$R_{\max }=\underset{c_j\&Element;Z_i}{\max }\left(\right|c_j-{\mathrm{\&mu;}}_i\left|\right)---\left(44\right)$

式中，i=1,2,…,K表示聚类后的索引号，Z_i是相应的类，N_p,i是类Z_i内的特征数目，c_j是属于类Z_i中的特征矢量，μ_i是Z_i类的均值矢量，表示为：

${\mathrm{\&mu;}}_i=\frac{1}{N_i}\underset{c_j\&Element;Z_i}{\mathrm{\&Sigma;}}{c}_j---\left(45\right)$

应用动态K均值方法确定聚类数目后，需要计算各个聚类的均值、方差和权值来初始化高斯混合模型；均值由式(49)求得；权值用聚类内的特征向量个数除以特征矢量总数来计算，其表达式为：

${\mathrm{\&pi;}}_i=\frac{N_{p,i}}{N_{\mathrm{frames}}}---\left(46\right)$

为了获得各个聚类的方差，为每个聚类构建N_p,i×N的矩阵，即：

C_i＝[c₁ c₂ … c_p,i]^T    (47)

式中，

是矩阵的第k列，是特征矢量的第k维，则方差矩阵的第m行n列元素的求法为：

${\mathrm{\&Sigma;}}_i^{(m,n)}=\mathrm{cov}(C_i^{(:,m)},C_i^{(:,n)})$

                           (48)

$=E\left[(C_i^{(:,m)}-M_m)\left((C_i^{(:,n)}-M_n\right)\right]$

式中，m,n=1,2,…,N，M_k是每个元素均为m_k的向量，式中m_k是聚类Z_i中第k维特征的均值，表示为：

$m_k=\frac{1}{N_{p,i}}\underset{j=1}{\overset{N_{p,i}}{\mathrm{\&Sigma;}}}{C}_i^{(j,k)}---\left(49\right)$

A22、精确估计高斯混合模型参数

高斯混合模型参数估计就是依据某种参数估计准则确定模型参数的过程，用极大似然准则进行估计；一首歌曲的特征矩阵表示为

这里N_frames表示帧的总数目，

是每帧的特征矢量，假设它们彼此独立，则高斯混合模型的似然函数定义为，在数理统计学中，似然函数是一种关于统计模型中的参数的函数,表示模型参数中的似然性：

$L(\mathrm{\&Theta;}/C)=p(C/\mathrm{\&Theta;})=\underset{i=1}{\overset{N_{\mathrm{frames}}}{\mathrm{\&Pi;}}}p(c_i/\mathrm{\&Theta;})---\left(50\right)$

将似然函数L(Θ/C)看成Θ的非线性函数，参数估计的目的是要找到使似然函数L(Θ/C)最大的模型Θ，即：

${\mathrm{\&Theta;}}^{*}=\arg \underset{\mathrm{\&Theta;}}{\max }p(C/\mathrm{\&Theta;})---\left(51\right)$

对式(51)取对数，获得对数似然函数

$J\left(\mathrm{\&Theta;}\right)=\log L(\mathrm{\&Theta;}/C)=\log \underset{i=1}{\overset{N_{\mathrm{frames}}}{\mathrm{\&Pi;}}}p(c_i/\mathrm{\&Theta;})---\left(52\right)$

实际中常用期望最大化算法来估计模型参数，于是得到模型权值为：

${\mathrm{\&pi;}}_k=\frac{1}{N_{\mathrm{frames}}}\underset{i=1}{\overset{N_{\mathrm{frames}}}{\mathrm{\&Sigma;}}}p\left(k\right|c_i,{\mathrm{\&Theta;}}^{\mathrm{old}})---\left(53\right)$

均值为：

${\mathrm{\&mu;}}_k=\frac{\underset{i=1}{\overset{N_{\mathrm{frames}}}{\mathrm{\&Sigma;}}}{c}_{i}p\left(k\right|c_i,{\mathrm{\&Theta;}}^{\mathrm{old}})}{\underset{i=1}{\overset{N_{\mathrm{frames}}}{\mathrm{\&Sigma;}}}p\left(k\right|c_i,{\mathrm{\&Theta;}}^{\mathrm{old}})}---\left(54\right)$

方差矩阵为：

${\mathrm{\&Sigma;}}_k=\frac{\underset{i=1}{\overset{N_{\mathrm{frames}}}{\mathrm{\&Sigma;}}}p\left(k\right|c_i,{\mathrm{\&Theta;}}^{\mathrm{old}})(c_i-{\mathrm{\&mu;}}_k){(c_i-{\mathrm{\&mu;}}_k)}^T}{\underset{i=1}{\overset{N_{\mathrm{frames}}}{\mathrm{\&Sigma;}}}p\left(k\right|c,{\mathrm{\&Theta;}}^{\mathrm{old}})}---\left(55\right)$

式中，p(k|c_i,Θ^old)称为后验概率；根据贝叶斯准则，该后验概率为：

$p\left(k\right|c_i,{\mathrm{\&Theta;}}^{\mathrm{old}})=\frac{p(k,c_i|{\mathrm{\&Theta;}}^{\mathrm{old}})}{p\left(c_i\right|{\mathrm{\&Theta;}}^{\mathrm{old}})}=\frac{{\mathrm{\&pi;}}_k^{\mathrm{old}}N(c_i;{\mathrm{\&mu;}}_k^{\mathrm{old}},{\mathrm{\&Sigma;}}_k^{\mathrm{old}})}{\underset{l=1}{\overset{M}{\mathrm{\&Sigma;}}}{\mathrm{\&pi;}}_l^{\mathrm{old}}N(c_i;{\mathrm{\&mu;}}_l^{\mathrm{old}},{\mathrm{\&Sigma;}}_l^{\mathrm{old}})}---\left(56\right)$

期望最大化估计的步骤为：

A221、根据式(41)和(56)得到后验概率；

A222、根据式(53)～(55)计算新的模型参数，以此更新上次的模型参数，当新旧模型参数几乎一致时，迭代停止，否则返回(a)步骤继续迭代计算；

B、生成推荐列表

生成推荐列表的流程包括提取输入歌曲的特征矢量、分别计算该歌曲与数据库中歌曲的加权相似度和对获得的相似度按降序排列产生推荐列表，具体步骤如下：

B1、提取音乐信号的特征矢量，提取步骤与步骤A1所述的步骤完全一致，只是输入音乐信号不同，生成推荐列表的输入音乐信号为用户所试听的一支曲目，经步骤A1处理的输出即为产生类似该曲目风格的音乐，而步骤A1中的输入音乐信号由建立数据库的许多音乐信号组成；

B2、分别计算该歌曲与数据库中歌曲的加权相似度

得到特征矢量的统计模型后，通过模型之间的相似度来判定歌曲的相似度；若M₁维的高斯混合模型模型A为：

${\mathrm{\&Theta;}}_A=\{({\mathrm{\&mu;}}_{1,A},{\mathrm{\&Sigma;}}_{1,A},{\mathrm{\&pi;}}_{1,A}),\&CenterDot;\&CenterDot;\&CenterDot;,({\mathrm{\&mu;}}_{M,A},{\mathrm{\&Sigma;}}_{M,A},{\mathrm{\&pi;}}_{M_1,A})\}---\left(57\right)$

对应的特征矩阵为：

$C_A={[c_{1,A},c_{2,A},\&CenterDot;\&CenterDot;\&CenterDot;,c_{N_{\mathrm{Aframes}},A}]}^T---\left(58\right)$

式中，N_Aframes为模型A特征矢量的个数；

另一M₂维的高斯混合模型模型B为：

${\mathrm{\&Theta;}}_B=\{({\mathrm{\&mu;}}_{1,B},{\mathrm{\&Sigma;}}_{1,B},{\mathrm{\&pi;}}_{1,B}),\&CenterDot;\&CenterDot;\&CenterDot;,({\mathrm{\&mu;}}_{M,B},{\mathrm{\&Sigma;}}_{M,B},{\mathrm{\&pi;}}_{M_2,B})\}---\left(59\right)$

对应的特征矢量为：

$C_B={[c_{1,B},c_{2,B},\&CenterDot;\&CenterDot;\&CenterDot;,c_{N_{\mathrm{Bframes}},B}]}^T---\left(60\right)$

式中，N_Bframes为模型B特征矢量的个数；

相似性的计算公式由下式获得：

r(A,B)＝logL(Θ_A/C_A)+logL(Θ_B/C_B)-logL(Θ_A/C_B)-logL(Θ_B/C_A)    (61)

式(61)的最大值为：

r_max(A,B)＝logL(Θ_A/C_A)+logL(Θ_B/C_B)    (62)

则模型A与模型B的相似度定义为：

$\mathrm{sim}(A,B)=\frac{r_{\max }(A,B)-r(A,B)}{r_{\max }(A,B)}=\frac{\log L({\mathrm{\&Theta;}}_A/C_B)+\log L({\mathrm{\&Theta;}}_B/C_A)}{\log L({\mathrm{\&Theta;}}_A/C_A)+\log L({\mathrm{\&Theta;}}_B/C_B)}---\left(63\right)$

为了计算歌曲m与歌曲n的相似度，利用式(63)分别计算上述三类特征的相似度，即sim_gfcc(m,n)、sim_osc(m,n)和sim_mood(m,n)，则两首歌曲的总相似度表示为：

sim(m,n)＝w_gfcc·sim_gfcc(m,n)+w_osc·sim_osc(m,n)+w_mood·sim_mood(m,n)    (64)

式中，w_gfcc、w_osc与w_mood为三类特征相似度的加权值，取值范围是0<w_gfcc<1、0<w_mood<1、0<w_osc<1，且w_gfcc+w_mood+w_osc=1；

B3、对获得的相似度降序排列，产生推荐列表。

2.根据权利要求1所述的一种基于相似性的音乐推荐方法，其特征在于：所述的预加重因子μ的最佳值为0.97；所述的一帧音乐信号的长度N_f的最佳值为2048；所述的帧间交叠的长度的最佳值为N₀=0.5N_f=1024；所述的初始相位φ最佳值为0；所述的伽马通倒谱系数的个数最佳值为L=26；所述的情绪模型特征矢量维数最佳值为M=12；所述的子频带的个数B的最佳值为8；所述的W最佳值为512帧；所述的三类特征相似度的加权值的最佳值分别为w_gfcc=0.6、w_mood=0.2和w_osc=0.2；所述的系数α的最佳值为0.2。

Images