MFCCとは

MFCCは聴覚フィルタに基づく音響分析手法で、主に音声認識の分野で使われることが多いです。

~~最近だとニューラルネットワークに学習させる音声特徴量としてよく使われていますね。~~
2019.5.29訂正
Deep Learning for Audio Signal ProcessingによるとDeep Learningにおいては必要な情報が失われるためMFCCは使わずに、最後の計算ステップである離散コサイン変換を省いたメルスペクトラム(log-mel spectrum)が使われるそうです。MFCCは従来手法である隠れマルコフモデル、混合ガウスモデル、サポートベクターマシンで使われることが多いです。

今回はMFCC「メル周波数」や「ケプストラム」についても説明し、具体的なMFCCの実装方法も見ていきたいと思います。

メル尺度

心理学者のStanley Smith Stevensらによって提案された、人間の音高知覚が考慮された尺度です。

1000Hzの純音の高さの感覚を1000メルと決めた上で、1000メルの半分の高さに感じた音を500メル、1000メルの2倍の高さに感じた音を2000メルという容量で定めたものです。

実験により、可聴域の下限に近い音は高め、上限に近い音は低めに聞こえることがわかっています。

メル尺度は基底膜上の座標とほとんど一致します。1000メルは基底膜の1cmに相当し、臨界帯域幅は約137メルに相当します。

周波数$f$をメル尺度$m$に変換する式は

m = m_0\log(\frac{f}{f_0} + 1)

逆変換は

f = f_0(\exp\frac{m}{m_0} - 1)

で与えられます。

$f_0$は自由パラメータの周波数パラメータであり、$m_0$は「1000Hzは1000メル」という制約から導かれる

m_0 = \frac{1000}{\log(\frac{1000\mathrm{Hz}}{f_0} + 1)}

であらわされる従属パラメータです。

よく用いられるパラメータとして

f_0 = 700, m_0 = 2595

があります。

周波数[Hz]とメル周波数[mel]の関係をグラフに表すとこのようになります。

ケプストラム

音声信号のフーリエ変換の絶対値に対数をかけ、さらに逆フーリエ変換（フーリエ変換とする場合もある）したもの。

信号$x(n)$のフーリエ変換を$X(e^{j\omega}$とすると、ケプストラム$c(m)$は

c(m) = \frac{1}{2\pi}\int_{\pi}^{\pi}\log|X(e^{j\omega})|e^{j\omega n}d\omega

高ケフレンシ領域は微細振動（音声でいうと声帯の振動）、低ケフレンシ領域はスペクトル包絡（音声でいうと声道フィルタの特性）を表します。

リフタリングという処理を施すことによってこれらを分離することができます。

変換手順

この記事を参考にしながら実際に変換してみます。

コードはほぼそのままですが、Python3対応させたり、変数名やパラメータをちょこちょこ変えたりしてます。

今回は「あ」の音声波形の定常部分のMFCCを求めてみましょう。（音声ファイルは各自用意してください。）

まずは必要なモジュールのインポートから

import numpy as np
import matplotlib.pyplot as plt
import soundfile as sf

1. 波形を適当な長さに分割し、窓関数をかけてFFTを行う

# 音声の読み込み
master, fs = sf.read('wavfiles/a_1.wav')

t = np.arange(0, len(master) / fs, 1/fs)

# 音声波形の中心部分（定常部）を切り出す
center = len(master)//2 #　中心のサンプル番号
cuttime = 0.04 # 秒
x = master[int(center - cuttime / 2 * fs):int(center + cuttime / 2 * fs)]
time = t[int(center - cuttime/2*fs): int(center + cuttime/2*fs)]

plt.plot(time * 1000, x)
plt.xlabel("time [ms]")
plt.ylabel("amplitude")
plt.show()

上のコードを実行するとこのようなグラフが描けます。

次に窓処理を行って振幅スペクトルを求めます。

# ハミング窓をかける
hamming = np.hamming(len(x))
x = x * hamming

# 振幅スペクトルを求める
N = 2048 # FFTのサンプル数
spec = np.abs(np.fft.fft(x, N))[:N//2]
fscale = np.fft.fftfreq(N, d = 1.0 / fs)[:N//2]

plt.plot(fscale, spec)
plt.xlabel("frequency [Hz]")
plt.ylabel("amplitude spectrum")
plt.show()

2.メルフィルタバンクをかける

フィルタバンクは複数のバンドパスフィルタを並べたものです。

MFCCでは三角形のバンドパスフィルタをオーバラップさせながら並べ、フィルタバンクとしたメルフィルタバンクを用います。

メルフィルタバンクは周波数領域での三角窓がメル尺度上で等間隔になるように並べたものであり、低周波ほど幅が狭く、高周波ほど幅が広くなります。

並べられたバンドパスフィルタの数をチャンネル数といいます。

まずは周波数とメル尺度を変換する関数を実装します。

def hz2mel(f):
    """Hzをmelに変換"""
    return 2595 * np.log(f / 700.0 + 1.0)

def mel2hz(m):
    """melをhzに変換"""
    return 700 * (np.exp(m / 2595) - 1.0)

次にメルフィルタバンクを作る関数を実装します。

def melFilterBank(fs, N, numChannels):
    """メルフィルタバンクを作成"""
    # ナイキスト周波数（Hz）
    fmax = fs / 2
    # ナイキスト周波数（mel）
    melmax = hz2mel(fmax)
    # 周波数インデックスの最大数
    nmax = N // 2
    # 周波数解像度（周波数インデックス1あたりのHz幅）
    df = fs / N
    # メル尺度における各フィルタの中心周波数を求める
    dmel = melmax / (numChannels + 1)
    melcenters = np.arange(1, numChannels + 1) * dmel
    # 各フィルタの中心周波数をHzに変換
    fcenters = mel2hz(melcenters)
    # 各フィルタの中心周波数を周波数インデックスに変換
    indexcenter = np.round(fcenters / df)
    # 各フィルタの開始位置のインデックス
    indexstart = np.hstack(([0], indexcenter[0:numChannels - 1]))
    # 各フィルタの終了位置のインデックス
    indexstop = np.hstack((indexcenter[1:numChannels], [nmax]))
    filterbank = np.zeros((numChannels, nmax))
    print(indexstop)
    for c in range(0, numChannels):
        # 三角フィルタの左の直線の傾きから点を求める
        increment= 1.0 / (indexcenter[c] - indexstart[c])
        for i in range(int(indexstart[c]), int(indexcenter[c])):
            filterbank[c, i] = (i - indexstart[c]) * increment
        # 三角フィルタの右の直線の傾きから点を求める
        decrement = 1.0 / (indexstop[c] - indexcenter[c])
        for i in range(int(indexcenter[c]), int(indexstop[c])):
            filterbank[c, i] = 1.0 - ((i - indexcenter[c]) * decrement)

    return filterbank, fcenters

# メルフィルタバンクを作成
numChannels = 20  # メルフィルタバンクのチャネル数
df = fs / N   # 周波数解像度（周波数インデックス1あたりのHz幅）
filterbank, fcenters = melFilterBank(fs, N, numChannels)

# メルフィルタバンクのプロット
for c in np.arange(0, numChannels):
    plt.plot(np.arange(0, N / 2) * df, filterbank[c])
    
plt.title('Mel filter bank')
plt.xlabel('Frequency[Hz]')
plt.show()

次に振幅スペクトルにメルフィルタバンクをかけます。

振幅スペクトルに対して作成した各フィルタバンクをかけ、フィルタリングされた振幅を足し合わせて対数を取ります。

forループを使って書いてもいいですが、内積の形でも書け、こっちの方がエレガントです。

# 振幅スペクトルにメルフィルタバンクを適用
mspec = np.dot(spec, filterbank.T)

元の振幅スペクトルとフィルタバンクをかけた後のスペクトルを比べてみましょう。

# 元の振幅スペクトルとフィルタバンクをかけて圧縮したスペクトルを表示
plt.figure(figsize=(13, 5))

plt.plot(fscale, 10* np.log10(spec), label='Original Spectrum')
plt.plot(fcenters, 10 * np.log10(mspec), "o-", label='Mel Spectrum')
plt.xlabel("frequency[Hz]")
plt.ylabel('Amplitude[dB]')
plt.legend()
plt.show()

3.離散コサイン変換を行う

メルフィルタバンクをかけた後のスペクトルをケフレンシ空間に戻すために離散コサイン変換を行います。

離散コサイン変換を行う理由として、係数間の相関が減り、特徴量としての性能が向上するなどがあります。

離散コサイン変換の結果から低次の係数を取り出すと、MFCCが得られます。大体12次まで取ることが多いです。

from scipy.fftpack.realtransforms import dct

ceps = dct(10 * log10(mspec), type=2, norm="ortho", axis=-1)
nceps = 12
mfcc = ceps[:nceps]
print(mfcc)

[46.02737864 32.79920243  4.57031165 -2.94696774 -1.54742386 -1.26479934
 -5.16211103 -0.62832775  5.74186803  4.19594078  0.79866862 -1.76984708]

LibrosaでMFCCを求める

上では音声データ全体の中の１フレームのみを用いてMFCCを求めましたが、Librosaを使うと簡単に各フレームごとのMFCCを求めることができます。

import librosa

x, fs = sf.read('wavfiles/a_1.wav')
mfccs = librosa.feature.mfcc(x, sr=fs)
print(mfccs.shape)
print(mfccs[0])

(20, 26)
[-299.48698558 -306.19008252 -281.77019349 -205.61561205 -165.07338147
 -182.40750356 -172.10454081 -175.49027956 -191.51493931 -210.77326972
 -217.29547751 -223.91082068 -232.23909119 -247.65589292 -256.91669445
 -263.66001621 -275.46132992 -299.72892089 -288.90436044 -262.39782466
 -252.39945491 -256.61596037 -260.47453287 -274.17839766 -299.54783037
 -342.87809203]

カラーマップで可視化することも出来ます。

import librosa.display

librosa.display.specshow(mfccs, sr=fs, x_axis='time')
plt.colorbar()

おわりに

今回はMFCCに入門してみました。
GitHubにJupyter Notebookを置いていますので、よかったらこちらもどうぞ。

参考

リアルタイム音声認識（安藤彰夫　著）

音楽と機械学習前処理編 MFCC ~ メル周波数ケプストラム係数

メル周波数ケプストラム係数（MFCC） - 人工知能に関する断創録

メル尺度 (Wikipedia)