实战音乐分类项目 -- 潘登同学的Machine Learning笔记

@

音乐结构解析

• 先欣赏一首歌

• 查看其波形图

查看无条件某一秒内的波形

#%%波形显示
import matplotlib.pyplot as plt
import librosa.display
import numpy as np
from pydub import AudioSegment

# 1秒=1000毫秒
SECOND = 1000
# 音乐文件
AUDIO_PATH = '无条件-陈奕迅(30秒).wav'

def split_music(begin, end, filepath):
 # 导入音乐
    song = AudioSegment.from_wav(filepath)

 # 取begin秒到end秒间的片段
    song = song[begin*SECOND: end*SECOND]

 # 存储为临时文件做备份
    a = filepath.split('.')
    temp_path = a[0] + '(波形显示).' + a[1]
    song.export(temp_path, format='wav')

    return temp_path

music, sr = librosa.load(split_music(0, 1, AUDIO_PATH))

# 宽高比为14:5的图
plt.figure(figsize=(14, 5))
librosa.display.waveplot(music, sr=sr)
plt.show()

• 1秒内的波形如下:

还是太恶心了, 看不出什么所以然, 换成0.1秒内的

接上代码

# 放大
n0 = 9000
n1 = 10000

music = np.array([mic for mic in music])
plt.figure(figsize=(14, 5))
plt.plot(music[n0:n1])
plt.grid()

# 显示图
plt.show()

• 0.1秒内的波形如下:

要不再看看正半轴的图像?

# 只显示正半段的
music = np.array([mic for mic in music if mic > 0])
plt.figure(figsize=(14, 5))
plt.plot(music[n0:n1])
plt.grid()

# 显示图
plt.show()

• 0.1秒内正半轴的波形如下:

傅里叶分析解析音频结构

时域与频域(预备知识)

时域: 就是我们现实生活中随着时间变动的那个方向, 如我们上面图展示的; 上面图的横坐标都是时间, 纵坐标表示振幅;

频域: 而频域就是不随时间变动的那个方向;

从我们出生，我们看到的世界都以时间贯穿，股票的走势、人的身高、汽车的轨迹都会随着时间发生改变。这种以时间作为参照来观察动态世界的方法我们称其为时域分析。而我们也想当然的认为，世间万物都在随着时间不停的改变，并且永远不会静止下来。但如果我告诉你，用另一种方法来观察世界的话，你会发现世界是永恒不变的，你会不会觉得我疯了？我没有疯，这个静止的世界就叫做频域。

傅里叶级数

傅里叶级数 Fourier series 最简单的理解方式就是任何周期函数都可以分解成一堆正弦函数，这里的正弦是 $A\sin(\omega x + \varphi)$. 首先这里的一堆可以是无穷个，然后又因为 $\sin(\alpha + \beta)=\sin\alpha\cos\beta + \sin\beta\cos\alpha$, 所以我们说可以把任何周期函数分解成一堆正弦与余弦函数。

看下面这一个例子:

如图是一个方波:

那么怎么用正弦波去做出这样一个方波呢?

• step1: 一个正弦波

$$y = \frac{4}{\pi}\sin x$$

• step2: 两个正弦波

在上面的基础上再加上一个正弦波 $$y = \frac{4}{3\pi}\sin 3x$$

将他们叠加 $$y = \frac{4}{\pi}\sin x + \frac{4}{3\pi}\sin 3x$$

• step3: 再叠加 $$y = \frac{4}{\pi}\sin x + \frac{4}{3\pi}\sin 3x + \frac{4}{5\pi}\sin 5x$$

• step?:再加再加 $$y = \frac{4}{\pi}\sin x + \frac{4}{3\pi}\sin 3x + \frac{4}{5\pi}\sin 5x + \cdots $$

当我叠加了50项的时候, 他就很接近方波了

而声波就可以用这样的傅里叶级数来描述, 下面我们就在时域和频域中展示傅里叶级数

图示傅里叶级数(以3项为例)

$$y = \frac{4}{\pi}\sin x + \frac{4}{3\pi}\sin 3x + \frac{4}{5\pi}\sin 5x$$

• 时域频域两边看

图中黑色的就是y, 其他颜色的就是他的组成部分, 在时域上看的很复杂的东西(y), 其实在频域上就是正弦波, 而这个些正弦波就是永恒不变的;

• 单看频域

但是很显然的, 只知道特定的频率的正弦波, 不一定能组成唯一的声波; 因为初相位不一样, 也就是初值不一定一样, 所以虽然后面的叠加都是周期性的;

我们可以再用一个叫做相位谱的东西来看看他的相位;

通过找一条基线, 我这里找的是时间为0, 作为基线;

然后找在这个基线后, 找到这个正弦波在该基线后的最大值所在的时间, 将该时间与基线的时间做差, 得到时间差, 再除周 期再乘$2\pi$

$$相位差 = \frac{时间差}{周期} \times 2\pi$$

然后我们说时域是从正面看, 频域是从侧面看, 相位谱就是从下面看

• 时间差

• 相位差

我们可以看到相位差其实是一样的, 因为我们的级数组成的初相位本来就是0;

$$ y = A\sin(\omega x + \varphi), \varphi 就是初相位$$

好了, 傅里叶级数就是这样刻画任何函数的, 其实也就那样, 对叭;

而傅里叶变换就是将一个傅里叶级数的结果(y), 拆成很多个正弦波函数;

然后上面所有的画图代码如下:

#%% 图示傅里叶级数
import matplotlib.pyplot as plt
import numpy as np


x = np.array(np.linspace(0.1, 10, 100))
y = [1] * 30 + [-1] * 30 + [1] * 30 + [-1] * 10
y1 = 4/np.pi*np.sin(x)
y2 = 4/(3*np.pi)*np.sin(3*x)
y3 = y1 + y2
y4 = 4/(5*np.pi)*np.sin(5*x)
y5 = y3 + y4
y6 = 4/(7*np.pi)*np.sin(7*x)
y7 = y5 + y6
plt.rcParams['axes.facecolor']='black'
plt.figure(figsize=(10,6))
plt.plot(x, y1)
plt.plot(x, y2, 'g')
plt.plot(x, y3, 'r')
plt.plot(x, y5, 'y')
plt.plot(x, y7, 'gray')
plt.show()

for i in range(3, 100, 2):
    y1 += 4/(i*np.pi)*np.sin(i*x)

plt.figure(figsize=(10,6))
plt.plot(x, y, 'gray')
plt.show()
#%% 时域频域两边看
plt.rcParams['axes.facecolor']='snow'
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.view_init(elev=30, azim=335)


ax.plot(x, [1]*len(x), y1, 'red')
ax.plot(x, [3]*len(x), y2, 'green')
ax.plot(x, [5]*len(x), y4, 'cyan')
ax.plot(x, [0]*len(x), y5, 'black')
plt.xlabel('时域')
plt.ylabel('频域')
ax.set_zlabel('振幅')
plt.show()
#%% 时间差
plt.rcParams['axes.facecolor']='snow'
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.view_init(elev=45, azim=330)

x1 = [i for i in x if i <= np.pi/2]
x2 = [i for i in x if i <= np.pi/6]
x4 = [i for i in x if i <= np.pi/10]

ax.plot([0]*6, range(0, 6), 0, 'blue')
ax.plot(x1, [1]*len(x1), 0, 'red')
ax.plot(x2, [3]*len(x2), 0, 'green')
ax.plot(x4, [5]*len(x4), 0, 'cyan')
# ax.plot(x, [0]*len(x), y5, 'black')
ax.set_xlabel('时域')
ax.set_xlim([0, 2])
ax.set_ylabel('频域')
ax.set_ylim([0, 5.5])
ax.set_zlabel('振幅')
ax.set_zlim([0, 1])
plt.show()

#%% 相位谱
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.view_init(elev=45, azim=330)

ax.plot([0]*6, range(0, 6), 0, 'blue')
ax.plot(x1, [1]*len(x1), 0, 'red')
ax.plot(x1, [3]*len(x2)*3, 0, 'green')
ax.plot(x1, [5]*len(x4)*5, 0, 'cyan')
# ax.plot(x, [0]*len(x), y5, 'black')
ax.set_xlabel('时域')
ax.set_xlim([0, 2])
ax.set_ylabel('频域')
ax.set_ylim([0, 5.5])
ax.set_zlabel('振幅')
ax.set_zlim([0, 1])
plt.show()

对上面的无条件做傅里叶变换

在python的spicy库中有一个叫fft的函数, 中文叫快速傅里叶变换, 因为不是电子信息专业的, 拿来用就行了, 不纠结他的算法了, 欧拉公式、虚数啥的也不好搞;

然后matplotlib库中有一个叫specgram的函数, 可以把时域和频域同时画出来;

这里只截取了30秒的片段, 因为多了也就那样, 只是做个示例;

然后音频剪切用的是~~Adobe Audition~~, 当然是python啦;

是pydub库中一个叫AudioSegment的类, 可以支持'mp3', 'wav', 'raw', 'ogg' or other ffmpeg的格式;

话不多说, 上代码!!!

• 音频剪辑代码

#%%音乐剪辑
from pydub import AudioSegment

# 1秒=1000毫秒
SECOND = 1000
# 导入音乐
song = AudioSegment.from_wav("无条件-陈奕迅.wav")

# 取33秒到63秒间的片段
song = song[33*SECOND:63*SECOND]

# # 入场部分提高6分贝, 退场部分减少5分贝
# ten_seconds = 10 * SECOND
# last_five_seconds = -5 * SECOND
# beginning = song[:ten_seconds] + 6
# ending = song[last_five_seconds:] - 5

# # 形成新片段
# new_song = beginning + song[ten_seconds:last_five_seconds] + ending
new_song = song
# 导出音乐
new_song.export('无条件-陈奕迅(30秒).wav', format='wav')

• fft傅里叶变换代码

#%% 无条件 傅立叶变换
from scipy import fft #fft是傅立叶变换
from scipy.io import wavfile
from matplotlib.pyplot import specgram
import matplotlib.pyplot as plt

path = r'C:\Users\潘登\Documents\python全系列\人工智能\无条件-陈奕迅(30秒).wav'
(sample_rate,x) = wavfile.read(path)
print(sample_rate,x.shape)   #sample_rate每秒采样多少样本   x是(1395072, 2)双元组表示这首歌是双通道 用x/sample_rate可以得到这首歌的时长（s)



def plotSpec(file_name):
    plt.subplot(1,2,1)
    sample_rate,x = wavfile.read(file_name)
    x = x[:,0] # 将双通道变成单通道
    specgram(x,Fs = sample_rate,xextent = (0,30))
    plt.xlabel('time')
    plt.ylabel('frequency')
    plt.grid(True,linestyle='-',color = '0.25')
    plt.title('时域--无条件-陈奕迅(30秒)')
    plt.subplot(1,2,2)
    plt.xlabel('frequency')
    plt.xlim(0, 4000)
    plt.ylabel('amplitude')
    plt.title('FFT of 无条件-陈奕迅(30秒)')
    plt.plot(fft(x,sample_rate))
    plt.show()


plt.figure(figsize = (18,9),dpi = 80, facecolor = 'w', edgecolor = 'k')
plotSpec(path)

结果如下:

音乐分类项目

前面铺垫了这么多了, 终于到了实战了;

• 目标：实现音乐分类, 将输入的音乐分为这几类:

$$\begin{bmatrix} classical & jazz & pop & blues & country\ metal & rock& hiphop & disco & reggae\ \end{bmatrix}$$

• 数据：每种歌曲各100首

先看歌曲的时域

很显然, 我们看不出上面所以然, 就一堆乱七八糟的声波;

from scipy import fft #fft是傅立叶变换
from scipy.io import wavfile
from matplotlib.pyplot import specgram
import matplotlib.pyplot as plt

(sample_rate,x) = wavfile.read(r'C:\Users\潘登\Documents\python全系列\人工智能\genres\blues\converted\blues.00000.au.wav')
print(sample_rate,x.shape)   #sample_rate每秒采样多少样本   x是(661794,)单元组表示这首歌是单通道 用x/sample_rate可以得到这首歌的时长（s)

# plt.figure(figsize = (10,4),dpi = 80)
# plt.xlabel('time')
# plt.ylabel('frequency')
# plt.grid(True,linestyle='-',color = '0.25')
# specgram(x,Fs = sample_rate,xextent = (0,30))
# plt.show()

def plotSpec(g,n):
    file_name = r'C:\Users\潘登\Documents\python全系列\人工智能\genres' + '\\' + g + '\converted' + '\\' + g + '.' + n + '.au.wav'
    sample_rate,x = wavfile.read(file_name)
    specgram(x,Fs = sample_rate,xextent = (0,30))
    plt.xlabel('time')
    plt.ylabel('frequency')
    plt.grid(True,linestyle='-',color = '0.25')
    plt.title(g+'_'+n[-1])

plt.figure(figsize = (18,9),dpi = 80, facecolor = 'w', edgecolor = 'k')
plt.subplot(6,3,1);plotSpec('classical','00001')
plt.subplot(6,3,2);plotSpec('classical','00002')
plt.subplot(6,3,3);plotSpec('classical','00003')
plt.subplot(6,3,4);plotSpec('jazz','00001')
plt.subplot(6,3,5);plotSpec('jazz','00002')
plt.subplot(6,3,6);plotSpec('jazz','00003')
plt.subplot(6,3,7);plotSpec('pop','00001')
plt.subplot(6,3,8);plotSpec('pop','00002')
plt.subplot(6,3,9);plotSpec('pop','00003')
plt.subplot(6,3,10);plotSpec('rock','00001')
plt.subplot(6,3,11);plotSpec('rock','00002')
plt.subplot(6,3,12);plotSpec('rock','00003')
plt.subplot(6,3,13);plotSpec('country','00001')
plt.subplot(6,3,14);plotSpec('country','00002')
plt.subplot(6,3,15);plotSpec('country','00003')
plt.subplot(6,3,16);plotSpec('metal','00001')
plt.subplot(6,3,17);plotSpec('metal','00002')
plt.subplot(6,3,18);plotSpec('metal','00003')

plt.tight_layout(pad = 0.4,w_pad = 0,h_pad = 1)
plt.show()

再看歌曲的频域

这个频域是有区分度的, 但是对我们来说还是太难了;

#%%图示傅里叶变换
def plotFFT(g,n):
    file_name = r'C:\Users\潘登\Documents\python全系列\人工智能\genres' + '\\' + g + '\converted' + '\\' + g + '.' + n + '.au.wav'
    sample_rate,x = wavfile.read(file_name)
    plt.plot(fft(x,sample_rate))
    plt.xlabel('frequency')
    plt.xlim(0,3000)   #太高频率的 我们听不见
    plt.ylabel('amplitude')
    plt.title(g+'_'+n[-1])

plt.figure(num = None,figsize = (10,8),dpi = 80,facecolor = 'w',edgecolor = 'k')
plt.subplot(6,2,1);plotSpec('classical','00001')
plt.subplot(6,2,2);plotFFT('classical','00001')
plt.subplot(6,2,3);plotSpec('jazz','00001')
plt.subplot(6,2,4);plotFFT('jazz','00001')
plt.subplot(6,2,5);plotSpec('country','00001')
plt.subplot(6,2,6);plotFFT('country','00001')
plt.subplot(6,2,7);plotSpec('pop','00001')
plt.subplot(6,2,8);plotFFT('pop','00001')
plt.subplot(6,2,9);plotSpec('rock','00001')
plt.subplot(6,2,10);plotFFT('rock','00001')
plt.subplot(6,2,11);plotSpec('metal','00001')
plt.subplot(6,2,12);plotFFT('metal','00001')
plt.show()

做傅里叶变换, 提取特征

将所有音频都做傅里叶变换, 将结果保存起来, 后面使用;

#提取特征
import numpy as np
def creat_fft(g,n):
    file_name = r'C:\Users\潘登\Documents\python全系列\人工智能\genres' + '\\' + g + '\converted' + '\\' + g + '.' + str(n).zfill(5) + '.au.wav'
    sample_rate,x = wavfile.read(file_name)
    fft_features = abs(fft(x)[:1000])  #不需要特别多的数据，再往后都是噪音
    sad = r'C:\Users\潘登\Documents\python全系列\人工智能\trainset' + '\\' + g + '.' + str(n).zfill(5) + '.fft'
    np.save(sad,fft_features)

genre_list = ['classical','jazz','country','pop','rock','metal']
for g in genre_list:
    for n in range(100):
        creat_fft(g,n)

将特征转化为x, 标签转化为y

• 对刚才的傅里叶变换后的特征设置为x, 做一些简单的筛选, 然后将所属的类别转化为数字编码, 设置为y;

• 划分训练集和测试集;

• 然后训练模型, 保存方便使用;

#%%读取傅里叶变换之后的数据集，将其转换成机器学习所需要的x和y
import numpy as np
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
import pickle
genre_list = ['classical','jazz','country','pop','rock','metal']
x = []
y = []
for g in genre_list:
    for n in range(100):
            file_name = r'C:\Users\潘登\Documents\python全系列\人工智能\trainset' + '\\' + g +  '.' + str(n).zfill(5) + '.fft.npy'
            fft_features = np.load(file_name)
            x.append(fft_features)
            y.append(genre_list.index(g))

x = np.array(x)
y = np.array(y)

x_train,x_test,y_train,y_test = train_test_split(x,y,test_size = 0.4,random_state=5)

#训练模型并且保存模型
model = LogisticRegression(multi_class='multinomial',solver = 'sag',max_iter=1000)#这里multinomial 最好与 sag一起使用，用默认的求解方法会报错
model.fit(x_train,y_train)

output = open('model.pkl','wb')  #保存模型到model.pkl  wb是以二进制写入
pickle.dump(model,output)
output.close()

读取模型, 进行预测

import pickle
from pprint import pprint 
from sklearn.metrics import confusion_matrix

pkl_file = open('model.pkl','rb')
model_loaded = pickle.load(pkl_file)
pprint(model_loaded)
pkl_file.close()

temp = model_loaded.predict(x_test)
print(confusion_matrix(y_test,temp,labels = range(len(genre_list))))
print(np.trace(confusion_matrix(y_test,temp,labels = range(len(genre_list))))/180)

但是这样的预测结果不是我们想要的, 因为一个数据集里面的, 要不就是很经典的, 要不就是调性一样的, 预测起来肯定结果会比较好, 于是选择自己找点歌曲来康康;

找点歌来试试

贝多芬的月光曲

还是老办法, 先把他剪辑, 套上面的代码就成;

代码如下:

#%%从网上下载音乐来查看model
print('Starting read wavfile...')

music_name = '月光曲(30秒).wav'
sample_rate,x = wavfile.read(music_name)

print(x.shape)
x = np.reshape(x,(1,-1))[0]  #把双通道的音频转化成单通道

test_fft_features = abs(fft(x)[:1000])

temp = model_loaded.predict([test_fft_features])
print(genre_list[int(temp)])

看到结果, 欧！！！！

再来一首, 无条件要有始有终

#%%从网上下载音乐来查看model
print('Starting read wavfile...')
music_name = '无条件-陈奕迅(30秒).wav'

sample_rate,x = wavfile.read(music_name)

print(x.shape)
x = np.reshape(x,(1,-1))[0]  #把双通道的音频转化成单通道

test_fft_features = abs(fft(x)[:1000])

temp = model_loaded.predict([test_fft_features])
print(genre_list[int(temp)])

看到结果, 欧！！！！

实战音乐多分类项目就在欢声笑语中结束了, 继续下一章吧！pd的Machine Learning