项目笔记（四） 实验——通过对midi自相似矩阵的高斯核对角卷积分割乐段

零、写在前面

要做seq2seq交响乐翻译，这大概是最后一块拼图了。读了些文献，参考了里面的方法，复现了一下，效果挺不错的。

一、原理

1.piano_roll与帧向量v

这是一个二维矩阵，两个轴分别为时间和音高，对应的数值表示力度。它把各个乐器的音都合成在一起了。可以调整的参数是fs：扫描频率。默认是fs=100，表示对于midi音频每秒扫描100次。为了在计算后面的自相似矩阵时减少计算量和内存，在不损失精度的情况下，fs可以降至fs=10，即每秒扫描十次。最后每0.1秒一帧，每一帧是一个128维（一共有128种音高）向量，每维的数值代表音符的力度（自然，没有音符的音高数值为0），第i帧的向量记为 vi v i ，总帧数记为Tx。

事实上，文献中的方法都是针对于音频而非midi文件：对于原始音频做常数Q变换或MFCC，得到频谱图后每一帧提取特征向量。这个特征向量就等价于上段中的帧向量。

2.自相似矩阵S

这是一个二维矩阵，两个轴的大小等于piano_roll扫描后的总帧数。对于自相似矩阵S的每一点，我们这样来计算：

S(i,j)=distance(vi,vj) S ( i , j ) = d i s t a n c e ( v i , v j )

这里引入了距离函数distance()，可选的包括欧氏距离，L2范数，余弦距离。尝试了前两者，几乎没有差别。显然地，两个相同的向量之间的距离（无论如何计算距离）应该恒等于常数0，所以，自相似矩阵的对角线应该为常数0，将自相似矩阵以灰度图的形式绘制出来时，可以能清楚地看到对角线的常数0。

作为例子，下面是土耳其进行曲的自相似矩阵灰度图：

3.高斯核对角卷积

这里的高斯核是一个二维矩阵，以灰度图的形式绘制（颜色越浅表示数值越高）：

它的大小是超参数，在本实验中使用 filter_size = 256，即对于一个时间点的计算，要考虑前后各12.8秒的上下文。它的标准差同样是一个超参数，这里设定其值为128。

对角卷积即为将高斯核的中心置于自相似矩阵的对角线上，从S(filter_size,filter_size)移动至S(总帧数Tx-filter_size, 总帧数Tx-filter_size) 进行卷积操作（这是因为考虑到乐段分割点不会位于离头尾很近的地方，所以没有进行补零填充）。卷积后得到的数值组成了一个Tx维的向量，称作新奇度得分N。对它进行滤波后就可以用作乐段分割了：极大值点对应分割点。

二、代码

首先导入所需要的库：

import numpy as np
import pretty_midi
import matplotlib.pyplot as plt
import keras
import scipy.signal
from scipy.ndimage import filters
import cv2
           

超参数设定等：

filter_size = 
standard_deviation = 
p = print
           

打开midi文件，获得piano_roll :

file_name = 'x.mid'
pm = pretty_midi.PrettyMIDI(file_name)
piano_roll = pm.get_piano_roll(fs=10)
Tx = len(piano_roll[0])

x = np.zeros(shape=(128,Tx))
x[:,:]=piano_roll
x=np.transpose(x)
           

计算自相似矩阵，这一步比较耗时，所以计算完后将矩阵导出，之后就不用重新计算了：

S=np.zeros((Tx,Tx))
p('正在计算自相似矩阵...')
for i in range(Tx):
    S[i]=(np.sum(np.sqrt(np.square(x-x[i])),axis=)) #欧氏距离

p(S)
np.save('S of '+ file_name + '.npy',S)
p('已储存相似度矩阵。')
           

绘制自相似矩阵：

plt.figure()
plt.imshow(S,cmap = plt.cm.gray)
plt.title('Similarity matrix')
plt.show()
           

下面计算高斯卷积核并绘图检查：

def gaussian_kernel_2d_opencv(filter_size,standard_deviation):
    kx = cv2.getGaussianKernel(filter_size,standard_deviation)
    ky = cv2.getGaussianKernel(filter_size,standard_deviation)
    return np.multiply(kx,np.transpose(ky))
filter =gaussian_kernel_2d_opencv(filter_size,standard_deviation)

ones = np.ones((filter_size,filter_size))
for i in range (filter_size):
    if i < filter_size/:
        ones[i,:int(filter_size/)]= -
    else:
        ones[i,int(filter_size/):]=-
filter=np.multiply(filter,ones)

plt.figure()
plt.imshow(filter,cmap = plt.cm.gray)
plt.show()
           

卷积操作：

N = np.zeros(Tx)
p('正在卷积计算新奇度函数...')
for i in range (Tx):
    if i>filter_size/ and (Tx-i)>filter_size/:
        slice_x = S[:, i-int(filter_size/) : int(i+filter_size/)]
        slice_xy = slice_x[i-int(filter_size/) : int(i+filter_size/)]
        N[i]=np.sum(np.multiply(slice_xy,filter))
p(N)
           

最后，我们对新奇度函数进行滤波：

这里滤波器的带宽仍然是一个超参数，设定为10时效果可以接受，曲线较为平滑，也没有很大的延迟。

p('正在对新奇度函数进行卷积滤波...')
filter_size_1d=
#filter_1d = cv2.getGaussianKernel(filter_size_1d,standard_deviation)
filter_1d = np.ones(filter_size_1d)
H = scipy.signal.convolve(N,np.reshape(filter_1d,(filter_size_1d)))
           

最后得到的H的极值就可以作为midi的乐段分割点了。

三、结果

我们以手工标注的贝多芬第七交响曲第二乐章做测试。

自相似矩阵：

我们将新奇度函数N、输出滤波后的新奇度函数H（橙色）和手工标注的分割点（绿色）绘至同一张图上，结果如下：

可以看到，这个模型的计算结果准确率尚可接受。