Pythonによる音声認識 MFCC編

前回記事では，フォルマント分析を用いて母音の認識を行いました．株式会社サイシードのインターンシップでMFCCに関して学び，分類精度の検証を行ったので，本記事は前回の続きとして音声認識によく使われているMFCCに関してまとめていきたいと思います．

##目次

• 背景
• MFCCとは
  • MFCCの導出過程
• MFCC導出プログラムの実装
  • librosaに関して
  • MFCCの導出
  • MFCCを用いた音素のクラス分類
• 連続音声認識に向けたプログラムの実装
  • 動的差分に関して
  • 一次差分，二次差分の導出
  • 二次差分の特徴まで含めた音素のクラス分類
• 考察
• まとめ
• 参考

背景

前回記事のフォルマント分析ではスペクトルの包絡を求め，その中で強調されている周波数を低い方からF1，F2として特徴に用いて母音の分析に用いていました．母音に関しては日本では5音素しかなく，また，基本的に声道の形に応じて音素が分類できたためフォルマント分析が有効でした．しかし，文章中の音素や子音に関しては動的変化や発声に使う器官の多さ・違いからフォルマント分析による特徴だけでは表現しきれません．そこで，今回は声道のスペクトルの包絡をより細かく表現することができるMFCCやその動的差分を用いた音素分類をしていきます．

MFCCとは

フォルマントは音声が声道を通る時に強調される複数の共鳴周波数を指し，母音は声道や舌の形により変化することから母音の分析に用いられていました．MFCCでは人の聴覚特性に合わせて音素を変換し，スペクトルの包絡をフォルマント分析より多くの特徴量を用いて表現することでフォルマント分析では捉えきれなかった音素の特徴を捉えることができます．
また，MFCC(メル周波数ケプストラム係数)はこの後細かく説明していくメルフィルタバンクの対数パワースペクトルのケプストラムのことを指します．

MFCCの導出過程

MFCCの導出手順は次のようにして行います．

1. 前処理
   1. 高域強調(プリエンファシス)
   2. 窓関数
2. パワースペクトルをメルフィルタバンクを用いて圧縮する．
3. 圧縮したパワースペクトルを対数パワースペクトルに変換．
4. 対数パワースペクトルを逆離散コサイン変換を適用してケプストラムを導出．
   　ケプストラムの低次成分(スペクトルの包絡を表す係数)を声道のスペクトルとして音声認識に利用

それぞれの詳細に関して説明していきます．

1. 前処理

まず，前処理に関しては音声のパワーが高域になるほど減衰するため，それを補償するために高域強調の処理を行います．加えて，不連続なデータに対してその波形の両端を減衰させるように窓関数をかけます．
音韻の違いにより音(強調される周波数)の大きさが変わることからパワースペクトルを用います(振幅スペクトルも音圧の表現方法の一つなのでそちらで解析している方もいると思います)．

2. パワースペクトルをメルフィルタバンクを用いて圧縮する．

人間の聴覚は周波数によって聞こえ方が変わります．高い周波数であればあるほど音の高さを聞き分けづらくなっていきます．実際の周波数と聴覚上の周波数の関係を実験的に求めたものをメル尺度と呼びます．
メル尺度で周波数軸を考えた時に，均等な幅になるようにメルフィルタバンクの個数に応じたフィルタを作成します．それを用いて，高速フーリエ変換して求めた音声波形のパワースペクトルをメルフィルタバンクの個数のグループに分けます．

3. 対数パワースペクトルに変換

メルフィルタバンクの出力(パワースペクトル)を対数パワースペクトルにするのには二つの理由があります．
一つ目に，人間が感じる音の大きさは音圧の対数に比例していて，音が大きくなるにつれて，音の大きさの違いを知覚しづらくなります．この人間の聴覚特性に対応させるために対数スペクトルにする必要があります．
二つ目に，ある時刻の音声は声門波(音源波)に声道のスペクトル(スペクトルの包絡)を畳み込んだものとなっているので，対数パワースペクトルを取ることによって声門波と声道のスペクトルの線形和に分離させることができます．
以下に，線形和に分離されることを式を追って説明します．

$y(n)$をある時刻の音声信号，$v(n)$をある時刻の声門波，$h(n)$を声道のインパルス応答とすると
$$y(n) = \sum_{m=-\infty}^{\infty}v(m) \cdot h(n-m) = v(n)*h(n)$$

音声信号をフーリエ変換すると
$$Y(k) = V(k) \cdot H(k)$$

となります．この時音声信号のパワースペクトル$S(k)$次のようになります．
$$S(k) = |Y(k)|^2 = |V(k) \cdot H(k)|^2 $$
これを対数変換すると
$$\log S(k) = 2\log |V(k)| + 2\log |H(k)| $$
となり，線形和に分離することができます．

4. ケプストラムの導出

得られた対数パワースペクトルを時間信号のように捉え，逆離散フーリエ変換をすることでケプストラムを得ることができます．
では，先ほど用いた対数パワー$logS(k)$を用いてケプストラムの導出過程を式で追っていきます．
逆離散フーリエ変換の式はこちらに載っています．
ケプストラム$c(n)$ は対数パワーの逆離散フーリエ変換を指しているので

\begin{align}
c(n) &= \frac{1}{N}\sum_{k=0}^{N-1}\log S_{k}e^{(i\frac{2\pi kn}{N})}\\
     &= \frac{1}{N}\sum_{k=0}^{N-1}\log S_{k} \cos (\frac{2\pi kn}{N}) + i \log S_{k} \sin (\frac{2\pi kn}{N})\\
\end{align}

また，このパワースペクトルはナイキスト周波数を境に線対称になっています．
つまり，

$$\frac{1}{N}\sum_{k=0}^{N-1}i \log S_{k}\sin (\frac{2\pi kn}{N}) = 0$$

となります．
よって，ケプストラム$c(n)$は次のように逆離散コサイン変換で表すことができます．

\begin{align}
c(n) &= \frac{1}{N}\sum_{k=0}^{N-1}\log S_{k} \cos (\frac{2\pi kn}{N})\\
     &= \frac{2}{N}\sum_{k=0}^{N-1}\log V_{k} \cos (\frac{2\pi kn}{N}) 
        + \frac{2}{N}\sum_{k=0}^{N-1}\log H_{k} \cos (\frac{2\pi kn}{N}) \\
\end{align}

$v(n)$は声門波なので複雑な変化(高周波成分)が多く含まれています．
一方で，$h(n)$は声道のインパルス応答なので滑らかな変化(低周波成分)が多く含まれています．
よって，このケプストラム$c(n)$の低次成分が声道のスペクトル$H_{k}$を表しています．

MFCC導出プログラムの実装

librosaに関して

librosaは音楽や音声を分析するためのPythonのパッケージになっています．今回のプログラムではMFCC，対数パワーの出力のために利用しています．詳細に関してはドキュメントがあるのでこちらを参照してください．

###MFCCの導出
使用するパッケージは次のものを用います．

#使用するパッケージ
import cis
import librosa
import sklearn
import numpy as np
from collections import defaultdict
import scipy.signal

librosaは pip install --upgrade sklearn librosa で自分の環境にインストールしてください．
また，本プログラムで用いている cis はPythonで学ぶ実践画像・音声処理入門に記載されているパッケージになります．こちらに置いてあるものを利用するか，音声ファイルの読み込みに使っているだけなので他のパッケージで代用してください．

MFCC導出までのプログラムは次のようになっています．

mfcc.py

#音声区間の全部の平均特徴をベクトルとして利用
mfcc_data = []
boin_list = ["a","i","u","e","o"]
nobashi_boin = ["a:","i:","u:","e:","o:"]
remove_list = ["silB","silE","sp"]

#高域強調
def preEmphasis(wave, p=0.97):
    # 係数 (1.0, -p) のFIRフィルタを作成
    return scipy.signal.lfilter([1.0, -p], 1, wave)
#mfccの計算
def mfcc(wave):
    mfccs = librosa.feature.mfcc(wave, sr = fs, n_fft = 512)
    mfccs = np.average(mfccs, axis = 1)
    #一次元配列の形にする
    mfccs = mfccs.flatten()
    mfccs = mfccs.tolist()
    #mfccの第1次元と14次元以降の特徴はいらないから消す
    mfccs.pop(0)
    mfccs = mfccs[:12]
    mfccs.insert(0,label)
    return mfccs
  
#データの読み込み, 音素毎にmfccを計算(使用データは500ファイル分)
for i in range(1,500,1):
    data_list = []
    open_file = "wav/sound-"+str(i).zfill(3)+".lab"
    filename = "wav/sound-"+str(i).zfill(3)#サンプリング周波数は16kHz
    v, fs = cis.wavread(filename+".wav")
    with open(open_file,"r") as f:
        data = f.readline().split()
        while data:
            data_list.append(data)
            data = f.readline().split()
        for j in range(len(data_list)):
            label =  data_list[j][2]
            if label in boin_list:
                start = int(fs * float(data_list[j][0]))
                end = int(fs * float(data_list[j][1]))
                voice_data = v[start:end]
                #短すぎるとうまく分析できないので飛ばす．
                if end - start <= 512:
                    continue
                # ハミング窓をかける
                hammingWindow = np.hamming(len(voice_data))
                voice_data = voice_data * hammingWindow
                p = 0.97
                voice_data = preEmphasis(voice_data, p)
                mfcc_data.append(mfcc(voice_data))

先ほどの説明の通り，高次成分に関しては声門波を表す特徴となっているため通常音素の認識には12次元程度用いられています．1次元目はデータの直交成分を表しているので除きます．
また，本プログラムで開いているlabファイルはjuliusのsegmentation-kitによって音素の開始位置・終了位置・音素の種類が入っているファイルになっています．前回記事に一例を載せているので気になる方はご確認ください．

###MFCCを用いた音素のクラス分類
音素の分類はSVM(サポートベクターマシン)を用いて比較していきたいと思います．

svm.py

#必要なパッケージ
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.svm import SVC
from sklearn.metrics import accuracy_score

#データセットの読み込み
df = pd.DataFrame(mfcc_data)

x = df.iloc[:,1:]#mfccで得た特徴点
y = df.iloc[:,0]#母音のラベル
#ラベルは一旦数字に変更
label = set(y)
label_list = list(label)
label_list.sort()
for i in range(len(label_list)):
    y[y == label_list[i]] =i
y = np.array(y, dtype = "int")

#教師データとテストデータの境目を決める
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size = 0.4, random_state = 1)

#データの標準化
sc = StandardScaler()
sc.fit(x_train)
x_train_std = sc.transform(x_train)
x_test_std = sc.transform(x_test)

#SVMのインスタンスを生成
model_linear = SVC(kernel='linear', random_state = 1)
model_poly = SVC(kernel = "poly", random_state = 1)
model_rbf = SVC(kernel = "rbf", random_state =1)
model_linear.fit(x_train_std, y_train)
model_poly.fit(x_train_std, y_train)
model_rbf.fit(x_train_std, y_train)

pred_linear_train = model_linear.predict(x_train_std)
pred_poly_train = model_poly.predict(x_train_std)
pred_rbf_train = model_rbf.predict(x_train_std)
accuracy_linear_train =accuracy_score(y_train, pred_linear_train)
accuracy_poly_train =accuracy_score(y_train, pred_poly_train)
accuracy_rbf_train =accuracy_score(y_train, pred_rbf_train)
print("train_result")
print("Linear : "+str(accuracy_linear_train))
print("Poly : "+str(accuracy_poly_train))
print("RBF : "+ str(accuracy_rbf_train))

pred_linear_test = model_linear.predict(x_test_std)
pred_poly_test = model_poly.predict(x_test_std)
pred_rbf_test = model_rbf.predict(x_test_std)
accuracy_linear_test = accuracy_score(y_test, pred_linear_test)
accuracy_poly_test = accuracy_score(y_test, pred_poly_test)
accuracy_rbf_test = accuracy_score(y_test, pred_rbf_test)
print("-"*40)
print("test_result")
print("Linear : "+str(accuracy_linear_test))
print("Poly : "+str(accuracy_poly_test))
print("RBF : "+ str(accuracy_rbf_test))

librosaのmfcc関数にかけると，入力した音素が複数に分割されます．今回，分割されたうちの音素の中央のみを使ったMFCCと音素ラベルの全体の平均を使ったMFCCを比較した時に平均を用いた方が精度がよかったのでそちらの結果を本記事では載せます．比較のため前回のフォルマント分析の結果も載せておきます.
クラス分類には8200程度の母音の音素データを使っています．

フォルマント分析による分類結果

train_result
Linear : 0.8109896432681243
Poly : 0.7206559263521288
RBF : 0.8550057537399309
----------------------------------------
test_result
Linear : 0.7825503355704698
Poly : 0.6932885906040268
RBF : 0.8308724832214766

MFCCによる分類結果

train_result
Linear : 0.885286271290786
Poly : 0.9113482454340243
RBF : 0.9201723784116561
----------------------------------------
test_result
Linear : 0.8833487226839027
Poly : 0.8913511849799939
RBF : 0.9039704524469068

フォルマント分析と比べ**約7％**認識精度は向上しました．
また，特徴量が増えたことによりLinearでも中々いい精度で認識できました．
子音に関しても同様に分類していきます．

フォルマント分析の結果

train_result
Linear : 0.2290266367936271
Poly : 0.20114513318396812
RBF : 0.31292008961911877
----------------------------------------
test_result
Linear : 0.22357723577235772
Poly : 0.1991869918699187
RBF : 0.30720092915214864

mfccの結果

train_result
train_result
Linear : 0.5635076681085333
Poly : 0.647463625639009
RBF : 0.679315768777035
----------------------------------------
test_result
Linear : 0.5396638159834857
Poly : 0.5364199351223827
RBF : 0.6004128575641404

データ量はフォルマント分析(約5000の子音データ)，mfccの結果(約8500の子音データ)程度です．
フォルマント分析と比較すると大幅な認識精度の向上になりますが，認識器としては今ひとつという結果になりました．

##連続音声認識に向けたプログラムの実装
実際の連続音声認識では音声フレームごとにMFCCなどの分析の処理を行い，そのフレームの音声に対する認識処理を行っていきます．そのため，音素区間毎で平均したMFCCの特徴を分類に用いるということは実際にはできません．

###動的差分に関して
一般的に，音声認識ではMFCCの特徴に加えて，メルフィルタバンクの出力の和の対数値(以後，対数パワースペクトル)，対数パワースペクトルとMFCCの一次差分(動的差分)，二次差分を加えた39次元を用いることが多いそうです(MLPシリーズ音声認識)．
動的差分を認識に用いることで子音や半母音のような動的に変化する音韻の特徴を捉えることができるようになります．
また，パワースペクトルに関しても母音と比べ子音は急激に変化するので動的差分を取ることで認識精度の向上につながります．

※ 注意
ここまでの実装では音素区間の特徴を平均したものを用いて精度検証していましたが，動的差分では連続的な特徴の変化を追う必要があります．連続データとして扱うために，ウィンドウ幅毎のMFCCの結果を特徴としてそのまま利用します．
また，それに当たりプログラムの修正も入ります．
結果に関しても上記の結果とは別物として以降の結果を見てください．

###一次差分, 二次差分の導出
では，早速プログラムの実装をしていきたいと思います．
まず，今回求める必要のあるMFCC，対数パワースペクトル，それらの動的差分を導出するための関数をまとめて以下に用意しておきます．

function.py

#高域強調
def preEmphasis(wave, p=0.97):
    # 係数 (1.0, -p) のFIRフィルタを作成
    return scipy.signal.lfilter([1.0, -p], 1, wave)

#mfccの計算
def mfcc(wave):
    mfccs = librosa.feature.mfcc(wave, sr = fs, n_fft = 512)
    mfccs = np.average(mfccs, axis = 1)
    #一次元配列の形にする
    mfccs = mfccs.flatten()
    mfccs = mfccs.tolist()
    #mfccの第1次元と14次元以降の特徴はいらないから消す
    mfccs.pop(0)
    mfccs = mfccs[:12]
    return mfccs

#対数パワースペクトルの計算
def cal_logpower(wave):
    S = librosa.feature.melspectrogram(wave, sr = fs,n_fft = 512)
    S = sum(S)
    PS=librosa.power_to_db(S)
    PS = np.average(PS)
    return PS

#前後何フレームを見るか(通常2~5)
#フレームに対する重み部分
def make_scale(K):
    scale = []
    div = sum(2*(i**2) for i in range(1,K+1))
    for i in range(-K,K+1):
        scale.append(i/div)
    return np.array(scale)

#差分特徴量の抽出
def make_delta(K, scale, feature):
    #自身の位置からK個前までのデータ参照
    before = [feature[0]]*K
    #自身の位置以降のK個のデータ参照
    after = []
    #差分特徴量の保管リスト
    delta = []
    for i in range(K+1):
        after.append(feature[i])      
    for j in range(len(feature)):
        if j == 0:
            match =  np.array(before + after)
            dif_cal =  np.dot(scale, match)
            delta.append(dif_cal)
            after.append(feature[j+K+1])
            after.pop(0)
        #後ろからK+1までは差分としてみる部分がある
        elif j < (len(feature) - K - 1):
            match = np.array(before + after)
            dif_cal = np.dot(scale, match)
            delta.append(dif_cal)
            before.append(feature[j])
            before.pop(0)                
            after.append(feature[j+K+1])
            after.pop(0)
        #データ量-K以降はafterにデータを追加できないため
        else:
            match = np.array(before + after)
            dif_cal = np.dot(scale, match)
            delta.append(dif_cal)
            before.append(feature[j])
            before.pop(0)
            after.append(feature[len(feature)-1])
            after.pop(0)
    return delta

次に，mfccと対数パワースペクトル，また，その二次差分までの特徴の抽出するためのプログラムを以下に記述します．

get_all_feature.py

#音声区間の全部の特徴をベクトルとして利用
phoneme = []
feature_data = []
delta_list = []
delta_2_list = []
nobashi_boin = ["a:","i:","u:","e:","o:"]
remove_list = ["silB","silE","sp"]


#データの読み込み, 音素毎にmfccを計算(使用データは500ファイル分)
for i in range(1,500,1):
    data_list = []
    open_file = "wav/sound-"+str(i).zfill(3)+".lab"
    filename = "wav/sound-"+str(i).zfill(3)#サンプリング周波数は16kHz
    v, fs = cis.wavread(filename+".wav")
    with open(open_file,"r") as f:
        data = f.readline().split()
        while data:
            data_list.append(data)
            data = f.readline().split()
        for j in range(len(data_list)):
            label =  data_list[j][2]
            if label not in remove_list:
                start = int(fs * float(data_list[j][0]))
                end = int(fs * float(data_list[j][1]))
                #伸ばし母音に関して
                if label in nobashi_boin:
                    label = label[0]
                voice_data = v[start:end]
                # ハミング窓をかける
                hammingWindow = np.hamming(len(voice_data))
                voice_data = voice_data * hammingWindow
                p = 0.97
                voice_data = preEmphasis(voice_data, p)
                mfccs = librosa.feature.mfcc(voice_data, sr = fs, n_fft = 512)
                mfccs_T = mfccs.T
                S = librosa.feature.melspectrogram(voice_data, sr = fs,n_fft = 512)
                S = sum(S)
                PS=librosa.power_to_db(S)
                for i in range(len(PS)):
                    feature = mfccs_T[i][1:13].tolist()
                    feature.append(PS[i])
                    feature_data.append(feature)
                    phoneme.append(label)
        K = 3
        scale = make_scale(K)
        delta = make_delta(K,scale,feature_data[len(delta_list):len(feature_data)])
        delta_list.extend(delta)
        second_delta = make_delta(K,scale,delta)
        delta_2_list.extend(second_delta)

これで二次差分まで作ることができたので，phoneme，feature_data，delta_list，delta_2_listを用いてデータフレームを作成してください．また，分類は先ほどのSVMのプログラムに形を合わせれば使えるかと思います．

###二次差分の特徴まで含めた音素のクラス分類
差分特徴量を加えて，SVMを用いてクラス分類した結果は次のようになります．
どちらの結果に関してもウィンドウ幅毎のMFCCの特徴を利用した認識精度です．
クラス分類には母音の音素データ（32500程度)程度を用いています．

MFCC(12次元)のみの結果

train_result
Linear : 0.8297360883797054
Poly : 0.8647708674304418
RBF : 0.8777618657937807
----------------------------------------
test_result
Linear : 0.8230149597238204
Poly : 0.8420406597621788
RBF : 0.8566168009205984

MFCC+対数パワーの二次差分まで加えた結果

train_result
Linear : 0.8631853518821604
Poly : 0.9549918166939444
RBF : 0.9495703764320785
----------------------------------------
test_result
Linear : 0.8588415803605677
Poly : 0.9132336018411967
RBF : 0.9177598772535481

動的差分を用いることで**約6%**の認識精度の向上になります．
連続的なデータを解析に用いているため前後の音からの変化が特徴として効果があることが確かめられました．

続いて，子音に関しても分類していってみましょう．
クラス分類には約27000の子音データを利用しています．

MFCC(12次元)のみの結果

train_result
Linear : 0.418983174835406
Poly : 0.5338332114118508
RBF : 0.5544989027066569
----------------------------------------
test_result
Linear : 0.4189448660510195
Poly : 0.4599067385937643
RBF : 0.49419402029807075

MFCC+対数パワーの二次差分を加えた結果

train_result
Linear : 0.5945501097293343
Poly : 0.8152889539136796
RBF : 0.8201658132162887
----------------------------------------
test_result
Linear : 0.5684374142817957
Poly : 0.6783395812379994
RBF : 0.7101581786595959

子音に関してもMFCC(12次元)の特徴を利用したものと比べて大きく精度の向上(約22%)が見られました．
子音の方が母音と比べて波形の変化が大きくなるため，動的差分が特徴としてより効いているのが確認できました．
一方で，特徴量が増えるにつれてtestとtrainの精度の差が広がっていることからデータ不足が問題として考えられます．

##考察
####母音に関して
母音に関しては，フォルマント分析と比較して，MFCCの利用によるスペクトルの包絡の特徴量を用いることで精度が向上することを確認しました．
文頭の方で，母音のみの場合はフォルマント分析で十分と記述していましたが，それは単一の母音の発音データに限っての場合です．しかし，今回母音データに関しては文章の発話データをjuliusを用いてsegmentしたものを用いていました．こういった母音データの場合は前後の音の影響を受け，波形が単一の母音を発音したものとは異なる部分が出てくる可能性があります．そのため，フォルマント分析では捉えきれていない波形の変化があり，MFCCを用いた方が精度が向上したと考えられます．さらに，前後の音からの変化の影響を受けるため動的差分を入れることによって精度がより向上したと推測されます．

####子音に関して
子音に関しては今回の結果から三つのことが言えます．

1. 子音は母音と比較してより多くの器官を用いて発音するため，フォルマント分析ではスペクトルの包絡により強調される周波数を見ているだけなので認識精度が悪かったことが窺えます．MFCCではスペクトルの包絡をより詳細に表現することができるため精度が向上したと考えられます．
2. 子音の音素波形の中には破裂音のように急激に波形が変化するものがあります．そういった音素に対して動的差分やパワースペクトルを用いることでMFCCと比べてさらに精度が向上したと推測されます．
3. データ不足に関してです．今回全体のデータ数(動的差分時のもの)が27000に対して子音の数は26種類(juliusによる分類)を用いているため，単純に計算して一つの音素に対しては1000強のデータしかありません．それに対して特徴量が39もあるのでデータ数が不足していることが挙げられます．また，今回扱った文章発話のデータは音素が同じようなバランスになるように調節されているものではありません．そのため，少ないと150程度のデータ量のものから多いと4000程度のデータ量のものまであります．データ量のばらつきの調整と全体的なデータの量を増やすことで精度は今回の結果よりさらによくなる可能性があります．

##まとめ
本記事では音響モデルで最もよく用いられているMFCCに関してその理論から実際の分類検証までを記載しています．
MFCCに関してはライブラリ(librosa)が用意されており，簡単に実装していくことができるので，その理論や使用結果の一例を細かく解説していきました．
音声認識の仕組みの理解やその実装の手がかりに本記事が繋がれば幸いです．

今後は連続的な音素に対する音響モデルや言語モデルの方に関しても触れていきたいと思います．

##参考

• Miyazaki's Pukiwiki
• 人工知能に関する断創録
• 音声認識(MLPシリーズ)