暗影中的声音（2）：手搓实验的一些记录

1
2
3
4
5
import csv
import numpy as np
from sklearn.cluster import KMeans
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal

1
2
3
N = 1000  # 样本数
D = 13     # 特征数
K = 5     # 高斯数

1
2
3
4
5
6
7
8
9
10
data = []
with open("aa.csv", "r") as f:
    r = csv.reader(f, delimiter = ",")
    for i, d in enumerate(r):
        if i < N:
            data.append(d)
        else:
            break
data = np.array(data).astype(np.float)
data.shape

1
(1000, 13)

1
2
3
4
z_n_k = np.zeros((N,K))  # 隐变量的后验概率
mu_k = np.zeros((K,D))   # 高斯均值
cov_k = np.zeros((K,D,D))# 协方差
pi_k = np.zeros((K,1))   # 均值矩阵

1
2
3
4
5
6
7
8
9
def k_means(data,N,K,D,MAX_ITERATIONS): # O(样本数*迭代次数)
    cov = np.zeros((K,D,D))    #  协方差
    pi = np.zeros((K,1))       # 高斯权值
    kmeans = KMeans(n_clusters = K).fit(data)   # 寻找均值，高斯核中心
    for k in range(K):
        cov[k] = np.cov(data[kmeans.labels_ == k].T)
        pi[k] = data[kmeans.labels_ == k].shape[0]/N
    print("均值：{0}，协方差：{1}，权值{2}".format(kmeans.cluster_centers_.shape,cov.shape,pi.shape))
    return kmeans.cluster_centers_,cov,pi

1
mu_k,cov_k,pi_k = k_means(data,N,K,D,100)  

1
2
3
4
5
6
def e_step(): # O(样本数*高斯数)
    num = np.zeros((K,1))
    for n in range(N):
        for k in range(K):
            num[k] = pi_k[k]* (multivariate_normal.pdf(data[n],mu_k[k],cov_k[k],allow_singular=True))
        z_n_k[n] = np.reshape(num/np.sum(num),(K,))

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
def special_op(array_1,array_2,N):
    for n in range(N):
        array_1[:,n] *= array_2[n]
    return array_1

def m_step():
    n_k = np.sum(z_n_k,axis = 0)
    #print("--->",n_k)
    # 更新均值 O(高斯数*样本数)
    for i in range(K):
        temp = np.zeros((1,D))
        for n in range(N):
            temp += z_n_k[n][i]*data[n,:]
        mu_k[i] = (1/n_k[i])*temp

    #print("--mu",mu_k)
    # 更新协方差 O(高斯数*样本数)
    for i in range(K):
        cov_k[i] = np.dot(special_op((data - mu_k[i,:]).T,(z_n_k[:,i]),N), (data-mu_k[i])) / n_k[i]
        cov_k[i] = cov_k[i] + 0.001*np.eye(cov_k[i].shape[0])
        
    #print("--cov",cov_k)
    # 更新权值 O(高斯数)
    for i in range(K):
        pi_k[i] = n_k[i]/N
    #print("--pi",pi_k)

1
2
3
4
5
6
7
8
def calculate_likelihood():
    log_likelihood = 0
    for n in range(N):
        temp = 0
        for k in range(K):
            temp += pi_k[k] * (multivariate_normal.pdf(data[n],mu_k[k],cov_k[k],allow_singular=True))
        log_likelihood += np.log(temp)
    return log_likelihood

1
2
3
4
5
6
7
8
9
10
11
12
13
iterations = 0
old_likelihood = 0
new_likelihood = 9999
likelihood = []

while(abs(old_likelihood - new_likelihood)>1 and iterations < 300):
    iterations += 1
    old_likelihood = new_likelihood
    e_step()
    m_step()
    new_likelihood = calculate_likelihood()
    likelihood.append(new_likelihood)
    print("iterations:{0}|likelihood:{1}|==>{2}".format(str(iterations),str(new_likelihood),abs(old_likelihood - new_likelihood)))

1
2
final_dict = {"mean":mu_k,"cov":cov_k,"pi":pi_k,"likelihood":likelihood}
np.save("ubm_file_name",final_dict)

1
2
3
4
5
6
def plot(z_n_k):
    plt.figure()
    color_array = ["r","g","b","y","black"]
    for i,each in enumerate(data):
        plt.scatter(each[0],each[1],c = color_array[np.argmax(z_n_k[i,:])],edgecolor = color_array[np.argmax(z_n_k[i,:])])
    plt.show()

1
plot(z_n_k)

1
2
3
import numpy as np
import csv
from scipy.stats import multivariate_normal

1
2
3
4
5
6
7
8
9
10
11
N = 1000
D = 13
K = 5

z_n_k = np.zeros((N,K))
mu_k = np.zeros((K,D))
mu_new = np.zeros((K,D))
n_k = np.zeros((K,1))
pi_k = np.zeros((K,1))
cov_k = np.zeros((K,D,D))
data = np.zeros((N,D))

1
2
3
4
ubm = np.load("ubm_file_name.npy").item()
mu_k = ubm["mean"]
cov_k = ubm["cov"]
pi_k = ubm["pi"]

1
2
3
4
5
6
7
with open("data/english1.csv") as f_data:
    reader_data = csv.reader(f_data)
    for i,each_data in enumerate(reader_data):
        if i == N:
            break
        data[i,:] = list(map(float,each_data))
data.shape

1
(1000, 13)

1
2
3
4
5
old_likelihood = 99999
new_likelihood = 0
likelihood_list = []
iterations = 0
final_dict = {}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
while(abs(old_likelihood - new_likelihood) > 1e-4 and iterations < 300):
    iterations += 1
    old_likelihood = new_likelihood
    print("epoch:",iterations)

    # 第一步：E步
    num = np.zeros((K,1))
    for n in range(N):
        for k in range(K):
            num[k] = pi_k[k] * (multivariate_normal.pdf(data[n],mu_k[k],cov_k[k]))
        z_n_k[n] = np.reshape(num/np.sum(num),(K,))
    print("finished E step...")  
    
    n_k = np.sum(z_n_k,axis = 0)
    n_k += 1e-10

    for i in range(K):
        temp = np.zeros((1,D))
        for n in range(N):
            temp += z_n_k[n][i]*data[n,:]
        mu_new[i] = (1/n_k[i])*temp
        
    # 第二步：MAP
    # 1.计算修正因子
    adaptation_coefficient = n_k/(n_k + 16.0)
    # 2.更新均值
    print("beginning adaptation...")
    for k in range(K):
        mu_k[k] = (adaptation_coefficient[k] * mu_new[k]) + ((1 - adaptation_coefficient[k]) * mu_k[k])
    
    # 计算对数似然值
    log_likelihood = 0
    print("calculating likelihood")
    for n in range(N):
        temp = 0
        for k in range(K):
            temp += pi_k[k] * (multivariate_normal.pdf(data[n],mu_k[k],cov_k[k]))
        log_likelihood += np.log(temp)
    new_likelihood = log_likelihood
    likelihood_list.append(log_likelihood)

    print(abs(old_likelihood - new_likelihood))
    print("************************")

1
2
final_dict = {"mean":mu_k,"cov":cov_k,"pi":pi_k,"likelihood":np.array(likelihood_list)}
np.save('map_file_name',final_dict)

1
2
3
4
5
6
7
import numpy as np
import python_speech_features as psf
from scipy.io import wavfile
from IPython.lib.display import Audio
import matplotlib.pyplot as plt
import pandas as pd
import time

1
2
3
4
5
6
NO_SPEAKERS = 4   # 说话人个数
'''读取wav'''
spk  = [wavfile.read('data/english' + str(i) + '.wav') for i in range(1, NO_SPEAKERS + 1)]
for i in range(NO_SPEAKERS):
    print("说话人 "+ str(i+1) + " 注册时长约为" + str(round(len(spk[i][1])/spk[i][0],2)))

1
Audio(spk[0][1], rate=spk[0][0])

1
2
3
4
5
6
samples = spk[0][1]
sample_rate = spk[0][0]
plt.figure(figsize=(12, 3))
plt.plot(np.arange(0, len(samples))/sample_rate, samples)
plt.xlabel("Seconds")
plt.show()

1
2
3
4
spk_mfcc = [psf.mfcc(spk[i][1])  for i in range(NO_SPEAKERS)]
print("MFCC的维度：" + str(spk_mfcc[1].shape[1]))
for i in range(NO_SPEAKERS):
    print("说话人 "+ str(i+1) + " 帧数为" + str(len(spk_mfcc[i])))

1
2
3
for i, speaker_mfcc in enumerate(spk_mfcc):
    average = sum(speaker_mfcc) / len(speaker_mfcc)
    spk_mfcc[i] -= average

1
2
3
4
5
6
7
8
9
10
11
12
In [7]:

data = []
lable = []
for i in range(NO_SPEAKERS):
    data.extend(spk_mfcc[i])
    lable.extend([i] * len(spk_mfcc[i]))
data = np.array(data)
lable = np.array(lable)
totals = pd.DataFrame(data)
totals['ID'] = lable
totals.head()

1
2
3
4
5
6
7
8
9
train_size = 0.9 #训练集大小
spk_train_size = []

# lable = []  # 还可以加label

for i in range(NO_SPEAKERS):
    spk_train_size.append(int(train_size * len(spk_mfcc[i])))
print("每个说话人的训练集大小（帧数）：", end=" ")
print(spk_train_size)

1
2
3
4
5
6
7
8
9
10
11
12
13
total_mfcc = []  # 总的背景语料
spk_start = []   # 每个说话人语料的开始
spk_end = []     # 每个说话人语料的结束

for i in range(NO_SPEAKERS):
    spk_start.append(len(total_mfcc))
    total_mfcc.extend(spk_mfcc[i][:spk_train_size[i]])
    spk_end.append(len(total_mfcc))
print("背景语料总帧数：" + str(len(total_mfcc)))
print("start:", end=" ")
print(spk_start)
print("e n d:", end=" ")
print(spk_end)

1
2
3
4
5
6
7
from sklearn.mixture import GaussianMixture
GMM = []
UBM = []
# 每一个人训练一个GMM 和 一个 UBM
for i in range(NO_SPEAKERS):
    GMM.append(GaussianMixture(n_components= 32, covariance_type= 'diag'))
    UBM.append(GaussianMixture(n_components= 32, covariance_type= 'diag'))

1
2
3
4
5
6
7
8
9
10
11
12
for i in range(NO_SPEAKERS):
    print(str(i+1) + "-GMM开始", end=",")
    start = time.time()
    GMM[i].fit(spk_mfcc[i][:spk_train_size[i]])
    end = time.time()
    print("耗时: {}s ".format(round(end-start,2)), end="--> ")
    print(str(i+1) + "-UBM开始", end=",")
    start = time.time()
    UBM[i].fit(total_mfcc[:spk_start[i]] + total_mfcc[spk_end[i]:])
    end = time.time()
    print("耗时: {}s ".format(round(end-start,2)), end="--> ")
    print()

1
2
3
4
5
6
7
8
9
10
11
for i in range(NO_SPEAKERS):
    print("说话人：{}".format(i+1))
    # 似然比检验
    x = GMM[i].score_samples(spk_mfcc[i][spk_train_size[i] + 2 : ]) - UBM[i].score_samples(spk_mfcc[i][spk_train_size[i] + 2 : ])
    total = 0
    correct = 0
    for i in x:
        if i > 0:
            correct +=1
        total += 1
    print("准确率（属于GMM的帧/总帧数） is {}".format(correct/total))

1
2
3
4
5
6
7
8
confusion = np.zeros((NO_SPEAKERS,NO_SPEAKERS),dtype=int)
for i in range(NO_SPEAKERS):
    for j in range(NO_SPEAKERS):
        x = GMM[j].score_samples(spk_mfcc[i][spk_train_size[i] + 2:]) - 
            UBM[j].score_samples(spk_mfcc[i][spk_train_size[i] + 2:])
        for score in x:
            if score > 0:
                confusion[i][j] += 1

1
2
3
4
class_names = ["ID_1","ID_2","ID_3","ID_4"]
from confusion_plt import plot_confusion_matrix
plot_confusion_matrix(confusion, classes=class_names,
                      title='Confusion matrix, without normalization')