摘要
副歌就是我们日常说的高潮。
副歌分析属于MIR领域的一个小分支，MIR（Music information retrieval)是从音乐中检索信息的跨学科科学，该领域需要数学、乐理、信号处理、机器学习、概率、算法等学科的背景知识。
业界类似方案
- 1：根据能量值、音高区间差异定位
- 2：根据动态歌词定位
- 3：谱聚类&能量值
第一种方式太简单暴力，很容易命中一些间奏，且起始点不准。
第二种依赖精准的歌词，对纯音乐无歌词无能为力。
本文使用的是第三种方式，是作者结合业界谱聚类算法以及自己对歌曲的理解整合的一个方案。下面一点点剖析，本文实验歌曲是著名实力加偶像派歌手不罕见的《你好再见》，不妨打开虾米听着歌曲看文章，更容易理解。
提取大步骤
- 1：解析频谱
- 2：谱聚类分段分类
- 3：分析副歌段落
基础知识
在开始进入正题之前，很有必要花1分钟了解下音频基础知识，老司机略过
- 响度[必须了解]：人耳感受到的声音强弱，它是人对声音大小的一个主观感觉量，响度的大小取决于音强、音高、音色、音长等条件
- 频率[必须了解]：是单位时间内完成周期性变化的次数，是描述周期运动频繁程度的量，例如1000hz/s。在音乐里可以简单类比大家说的音高。所以很多KTV 遇到高音就大声喊是没用的，那是响度，只会吓唬人，而频率没上去，俗称跑调。另外一个很重要的点就是*人耳能听到的平均频率范围，大概是20Hz-20000Hz
- 采样率[必须了解] ：每秒从连续信号中提取并组成离散信号的采样个数，它用赫兹（Hz）来表示。平时我们听得大部分高品质或无损歌曲都是44.1KHz采样频率，那么为什么大部分录音采样率是44.1KHz采样频率，因为根据人耳能听到的频率段来设定的，至于为什么乘以2了可以搜下Nyquist采样定理
- 码率：音频的码率可以简单理解每秒中包含的数据量，例如192bps，意思就是每秒钟存储了192b的数据
- 采样位深度：也叫采样位深，音频的位深度决定动态范围。它关乎到一个采样点能多精确地描述被采样的数据，也会影响到听歌时的细节。同比可以看图片里一个像素存储的RGBAlpha四种数据，对应位深度可能是32bit。歌曲大部分都是16bit，但目前高质量的数字音频系统已经使用24-32bit的量化精度。而有些对音质要求较低的场合，比如网络电话，也可能使用8bit
- 所以由以上几点可以推导出 ：音频大小b ≈采样率Hz×采样位大小bit×声道数x时间s
- wav：是微软公司开发的一种声音文件格式，是一种无损的音频文件格式
- mp3：有损压缩，牺牲部分高频（12KHz到16KHz）。
其实平时大家理解的无损仅仅是相对的，很多高逼格耳机音箱，如果音频源不好，也只能装逼用；
- 共振峰[必须了解]：共振峰是指在声音的频谱中能量相对集中的一些区域，共振峰不但是音色的决定因素，而且反映了声道（共振腔）的物理特征。例如刘德华和赵本山唱歌声音不一样，很大程度由共振峰决定。共振峰值携带了声音的辨识属性，如同人的身份证
- 频谱包络[必须了解] ：频谱包络是将不同频率的振幅最高点连结起来形成的曲线，就叫频谱包络线，是语音合成技术的重要组成部分
1：解析频谱
y, sr = librosa.load('audio/你好再见.mp3')平时我们听到的所有声音都可以理解为一种波，例如一个音叉发出频率为440HZ 的声波，他的波形可能是这样的：
傅里叶变换FFT
然鹅，一首歌曲里有太多各种各种的频率的波交叠在一起，就像如下这个gif 图
那么如果我们直接去分析这个交叠后的波形会发现很多有价值的数据都没有了，比如频率、相位。
那么，怎么把揉在一起的波形再拆开呢？这个时候傅里叶变换FFT就出现了，他可以把叠在一起的不同频率波形分开，就如同一束光透过棱镜后被分出各种颜色的光。图
傅里叶变换FFT是将按时间或空间采样的信号与按频率采样的相同信号进行关联的数学公式。在信号处理中，傅里叶变换可以揭示信号的重要特征（即其频率分量）
它是个伟大的公式，不仅仅用在音乐领域，还试图预测即将到来的地震，识别距离遥远的星系的组成部分，寻找热量大爆炸残余物中的新物理成分，从x射线衍射模式揭示蛋白质的结构，为NASA分析数字信号，研究乐器的声学原理，改进水循环的模型，寻找脉冲星（自转的中子星），用核磁共振研究分子结构。傅里叶变换已经被用于通过破译油画中的化学物质，来识别假冒的绘画。
下图黑色波形就是被混合后的，其他则是分解之后的。
STFT -短时傅里叶变换
而FFT 只能分离瞬间一个点的信号，只变换出了频域，但声音毕竟是连续的有时域概念的，且歌曲并不是平稳信号，如果说通过对信号加窗处理，信号在一窗内可看成平稳信号（可以简单理解频率不变），然后再用傅里叶变换。这样就能得到既有频域又有时域的信号，所以又有了 STFT -短时傅里叶变换
经过STFT 变换后的信号可以看下这个图：
到波形层面如下
常数Q转换
好，我们是不是可以立即去对歌曲做短时傅里叶变换拿到数据进行分析了？因为在频域一个点上的频率可以从0到20000，这个维度去计算一来代价很大，二来也不太符合人耳听到后的感受。我们听歌听到的音符也就是音阶，这些都可以对应到钢琴上的某个键。所以我们不妨把STFT变换后的数据在对频域进行log 压缩（为什么是log ,因为频率和音阶的关系也是这样的），把频域数据归属到若干个音阶上。这样一首歌曲经过这样变化后就可以得到如下这种图
刚才的变换学术上就叫常数Q转换（Constant Q transform），与短时距傅立叶转换一样为重要时频分析工具，其中特别适用于音乐信号的分析，这个转换产生的频谱最大的特色是在频率轴为对数标度而不是线性标度，且窗口长度会随着频率而改变。
经过CQT后的数据里表示能量值的维度，还需要转换为DB，因为这样可以使分布跨度很大的能量值相对集中到一个小区间内，不至于使太低的能量被忽略掉。
C = librosa.amplitude_to_db(np.abs(librosa.cqt(y=y, sr=sr,bins_per_octave=BINS_PER_OCTAVE,n_bins=N_OCTAVES * BINS_PER_OCTAVE)),ref=np.max)MFCC（Mel频率倒谱系数）
mfcc是从频谱提取的一种特征，这个特征在自动语音和说话人识别中广泛的使用一首歌曲我们可以得到它的频谱包络（连接所有共振峰值点的平滑曲线，共振峰值携带了声音的辨识属性，如同人的身份证），但是对于人类来说，人类听觉的感知只会聚焦在某些特定的区域而不是整个频谱包络，所以为了最大程度让数据符合人耳听到的频率变化，MFCC 就出现了。
MFCC 的具体提取原理网上有太多介绍，简单说它是先做信号预处理(预加重、分帧、加窗)、FFT、对数、DCT。
我们提取一首歌的mfcc
mfcc = librosa.feature.mfcc(y=y, sr=sr)
Msync = librosa.util.sync(mfcc, beats)降维之节拍分隔取均值
刚才利用CQT 算是对频域维度很好的进行了降维，但是时域维度上看，假如一首歌有4分钟，采样率为44100HZ/S,那么意味着时域上有46044100个采集点，可以想象下这个矩阵的维度，这样计算代价是非常大的，回到我们副歌提取的目标上看，完全没必要使用这么多的数据，怎么可以更低维的数据呢？
一种办法是用较低的采样率重新采样一次，但即便这样，数据维度仍旧很大，如果强力降低采样率，会直接影响分析结果，因为降低到一定程度很可能连人也听不出来了，何况机器。
另一种比较取巧的方法是获取一定时间段的平均值，无外乎就平均算法和中位数算法，考虑到歌曲的动态幅度，采用中位数比较稳妥。
那是简单段怎么切分比较好呢？节拍，对，没有比它更适合的了。
接下来就是计算歌曲的节拍分布（这个技术是另外一个话题，本文不展开），然后每个拍之间求频谱中位数(频率维度)
左图是没有节拍分隔降维的，右图是降维后的，可以发现从图纹理看，差别不大，也就是信息还原度比较高。当然，数据量已经减少了很多很多了
tempo, beats = librosa.beat.beat_track(y=y, sr=sr, trim=False)
Csync = librosa.util.sync(C, beats, aggregate=np.median)2：谱聚类分段分类
这一步骤很关键，也较多点，我们必须先了解下什么是谱聚类。
谱聚类（spectral clustering）是描述图的一种矩阵，在降维，分类，聚类等领域有广泛的运用。谱聚类是一种基于图论的聚类方法。将样本看成顶点，样本的相似度看作带权边。这样，把样本集划分成K个簇的过程就等同于一个图的分割问题。
如上图所示，每个顶点是一个样本。边的权值就是样本之间的相似度。然后我们需要分割，分割之后要使得连接不同组之间的边的权重尽可能低（组间相似度小），组内部的边的权重尽可能高（组内相似度高）。
是不是看着挺简单的？其实做起来并不简单，它要比传统的Kmean 聚类的复杂些，基本步骤如下：
- 1 根据数据构造一个Graph，Graph的每一个节点对应一个数据点，将各个点连接起来（随后将那些已经被连接起来但并不怎么相似的点，通过cut/RatioCut/NCut 的方式剪开），并且边的权重用于表示数据之间的相似度。把这个Graph用邻接矩阵的形式表示出来，记为 W。
- 2 把W的每一列元素加起来得到N个数，把它们放在对角线上（其他地方都是零），组成一个N*N的对角矩阵，记为度矩阵D，并把 W - D 的结果记为拉普拉斯矩阵 L = D - W。
- 3 求出L的前个k特征值（前k个指按照特征值的大小从小到大排序得到），以及对应的特征向量。
- 4 把这k个特征（列）向量排列在一起组成一个N*k的矩阵，将其中每一行看作k维空间中的一个向量，并使用 K-means 算法进行聚类。聚类的结果中每一行所属的类别就是原来 Graph 中的节点亦即最初的N个数据点分别所属的类别。
或许你已经看出来，前面的拉普拉斯矩阵的构造是为了将数据之间的关系反映到矩阵中，那么计算特征值以及特征向量从而达到将维度从N维降到k维，然后维度降下来了，我们就可以轻松的使用kmeans聚类，聚类完成后再将数据投影到原始数据上，这样原始数据就成功聚类啦
- 谱聚类算法的主要优点有：
1）谱聚类只需要数据之间的相似度矩阵，因此对于处理稀疏数据的聚类很有效。这点传统聚类算法比如K-Means很难做到
2）由于使用了降维，因此在处理高维数据聚类时的复杂度比传统聚类算法好。
- 谱聚类算法的主要缺点有：
1）如果最终聚类的维度非常高，则由于降维的幅度不够，谱聚类的运行速度和最后的聚类效果均不好。
2) 聚类效果依赖于相似矩阵，不同的相似矩阵得到的最终聚类效果可能很不同。
类似于谱聚类的算法除了kmean ,还有 kmedoids DBSCAN 等（注意KNN 不是聚类算法，它是分类算法，和PCA 支持向量机等都属于机器学习分类领域）
到这里谱聚类原理就介绍完了，你可以疑惑谱聚类不就是求一个把一个大图分隔出多个子图的RatioCut的解吗？是的，RatioCut 经过推导后，神奇的把一个NP难度的问题转换成拉普拉斯矩阵特征值（向量）的问题，数学就是这么美。
2.1：频谱自相关邻接矩阵
通过上面的介绍，我们已经了解了什么是谱聚类，那么接下来我们要构造关键的拉普拉斯矩阵 L = D - W中的W；
我们的目标是找到副歌段落，那么首先要把歌曲分成不同的段落，而歌曲往往是段落间有重复规律的。利用刚才得到的CQT数据Csync，我们来做下自相关，自相关可以理解为自己和自己逐帧对比，得出帧之间的相关度（底层使用的knn算法），生成邻接矩阵 R
R = librosa.segment.recurrence_matrix(Csync, width=3, mode='affinity',sym=True)
2.2 加强对角线
为了得到更连续平滑的重复片段，可以采用加窗的方法，求窗内的中位数，这样可以使断掉的片段连续光滑起来，还可以减低噪点的影响。这个窗的长度参数很关键，不能太大或太小，可以多次尝试取个合适的值。
df = librosa.segment.timelag_filter(scipy.ndimage.median_filter)
Rf = df(R, size=(1, 7))2.3 相邻关系邻接矩阵
歌曲段落除了在音调上有规律外，还有一个就是音色，通俗的说就是不同乐器或人的声音特征，往往一段内是同一种特征, 而提到这个就用到了上文提到的MFCC。
这个mfcc数据不适合再去做自相关，它表现出的不同音色的区别，更适用于分段，那么我们可以计算前一帧和当前帧的相似度(这里使用高斯核函数)，然后做矩阵变换构造出邻接矩阵R_path。
mfcc = librosa.feature.mfcc(y=y, sr=sr)
Msync = librosa.util.sync(mfcc, beats)
path_distance = np.sum(np.diff(Msync, axis=1)**2, axis=0)
sigma = np.median(path_distance)
path_sim = np.exp(-path_distance / sigma)
R_path = np.diag(path_sim, k=1) + np.diag(path_sim, k=-1)
2.4组合邻接矩阵
自相关邻接矩阵R和临近帧距离邻接矩阵R_path都已经有了，那么毕竟求拉普拉斯矩阵时只需要一个邻接矩阵，所以我们要把两者合一，做归一化后进行融合。这个里融合其实也很关键，两者权重非常影响结果，这一步我还有些疑问。
deg_path = np.sum(R_path, axis=1)
deg_rec = np.sum(Rf, axis=1)
mu = deg_path.dot(deg_path + deg_rec) / np.sum((deg_path + deg_rec)**2)
A = mu * Rf + (1 - mu) * R_path
2.5 求拉普拉斯矩阵
L = scipy.sparse.csgraph.laplacian(A, normed=True)2.6 求特征向量
evals, evecs = scipy.linalg.eigh(L)
evecs = scipy.ndimage.median_filter(evecs, size=(9, 1))
Cnorm = np.cumsum(evecs**2, axis=1)**0.52.7 利用kmean对特征向量聚类
//这个参数很关键，它决定了要分几种类别，太大太小都不好
k = 5
X = evecs[:, :k] / Cnorm[:, k-1:k]
KM = sklearn.cluster.KMeans(n_clusters=k)
seg_ids = KM.fit_predict(X)谱聚类后画出分段的情况如下图，这个就是歌曲《你好再见》的频谱图被聚类分段后的样子
3：分析副歌片段
经过以上步骤已经把歌曲分为了几个片段，并且片段之间进行了分类。那么到底哪一个分类或者哪一个片段属于副歌呢？这个还没有很管用成熟的算法可以使用，但是可以根据副歌一些特点来进行甄别。最明显的副歌片段在整首歌里的能量值是偏高的，顺着这个思路我们可以往下走。
3.1 剔除过短时间片段，例如可以认为10秒以内的都是无效片段，极少有10s 以内的副歌的
3.2 计算每拍DB均值，然后计算每段的能量(每段有n拍)
meanCsync = np.mean(Csync,0)
#根据节拍算均值
for bf ,bt , label in zip(zip(bound_beats, bound_beats[1:]), zip(bound_times, bound_times[1:]),bound_segs):
......
energyListSorted = sorted(timeEnergyMap.values())3.3 找到最大能量对应的label 值(也就是谱聚类出来的那个分类的label)
targetLabel = energyLabelMap.get(energyListSorted[-1])3.4 找到这个label 里按时间维度最后的那个片段，副歌大多在后面，这样命中的概率更高，副歌质量更高
fugeSeg = labelTimelistMap.get(targetLabel)[-1]
// the fuge time is 165.326077098-214.157641723
总结展望
到这里一首歌曲的副歌片段就定位出来了，我通过对大量歌曲的人工核对，准确率在80%以上，bad case 中有很多是很怪的歌曲，人工也听不出副歌在哪里；理论上正常歌曲的准确率在90%以上。
后面为了提高准确率，还可以加上帧人声识别和歌词定位，通过这两种方式可以更大概率的定位到副歌片段。
当然，我个人认为终极方案还是要结合机器学习和音频算法相结合的方式。非常感谢这期间@屠零 博士给予音频信号领域的指导答疑~
参考论文 http://bmcfee.github.io/papers/ismir2014_spectral.pdf
常用工具包librosa sklearn msaf Essntia等
