MFCC (メル周波数ケプストラム係数)

# MFCC
　MFCCの手順を簡潔にまとめた。

　実際に使用する際はlibrosaなどのライブラリを用いて１行で実装するのがいいと思う。

## MFCCとは
　音声認識で使用される特徴量抽出の方法。

## 流れ
1. wavファイルから音声波形を読み取る
2. 短い時間で切り出す
3. 高周波数を強調する
4. 窓関数を適応する
5. フーリエ変換する
6. メルフィルタバンクをかける
7. 離散コサイン変換

## 1.wavファイルから音声波形を読み取る

![Figure_1.png](https://qiita-image-store.s3.amazonaws.com/0/221435/2c83f5c2-259f-a84e-3ca1-3dcc4a61aff7.png)


```python3
def read_wave(file_name):
    """
    パス先のwaveファイルを読み込む
    """
    wf = wave.open(file_name, "r")
    fs = wf.getframerate()
    x = wf.readframes(wf.getnframes())
    x = np.frombuffer(x, dtype="int16") / 32768.0
    wf.close()
    return x, float(fs)

wav, fs = read_wave(DATA_DIR + "a.wav")
```

## 2.短い時間で切り出す

![Figure_2.png](https://qiita-image-store.s3.amazonaws.com/0/221435/dca9caa9-80e1-5961-1fde-69667a0c70d2.png)

```python3
def cut_wave(wave_data, cuttime=0.04):
    center = len(wave_data) / 2  # 中心のサンプル番号
    cuttime = 0.04         # 切り出す長さ [s]
    wav = wave_data[int(center - cuttime/2*fs) : int(center + cuttime/2*fs)]

    return wav


wave_data, fs = read_wave(DATA_DIR + "a.wav")
wav = cut_wave(wave_data)
```
※グラフを出力するコードは省略しています。

## 3.高周波数を強調する

![Figure_3.png](https://qiita-image-store.s3.amazonaws.com/0/221435/c1ada63c-a0fd-d82a-4508-50e15819ecf2.png)

```python3
def preEmphasis(wave, p=0.97):
    """プリエンファシスフィルタ"""
    # 係数 (1.0, -p) のFIRフィルタを作成
    return scipy.signal.lfilter([1.0, -p], 1, wave)

wave_data, fs = read_wave(DATA_DIR + "a.wav")
wav = cut_wave(wave_data)
wav = preEmphasis(wav)
```


## 4.窓関数を適応する


![Figure_4.png](https://qiita-image-store.s3.amazonaws.com/0/221435/96b73b67-39eb-d937-e821-cb2944e747e4.png)

```python3
def window(wave):
    """
    窓関数
    """
    hanningWindow = np.hanning(len(wave))
    wave = wave * hanningWindow
    return wave

wave_data, fs = read_wave(DATA_DIR + "a.wav")
wav = cut_wave(wave_data)
wav = preEmphasis(wav)
wav = window(wav)
```

## 5.フーリエ変換する

![Figure_5.png](https://qiita-image-store.s3.amazonaws.com/0/221435/636ae07b-cc56-1f76-5238-8d04cadd487c.png)

```python3
def fft(wave):

    """
    波形に対してフーリエ変換を行う関数

    :param: wave_data

    :return Adft: 振幅スペクトル
    :return Pdft: パワースペクトル
    """

    nfft = 2048  # FFTのサンプル数
    # 離散フーリエ変換
    dft = np.fft.fft(wave, n)
    # 振幅スペクトル
    Adft = np.abs(dft)[:int(nfft/2)]
    # パワースペクトル
    Pdft = np.abs(dft)[:int(nfft/2)] ** 2

    return Adft, Pdft

wave_data, fs = read_wave(DATA_DIR + "a.wav")
wav = cut_wave(wave_data)
wav = preEmphasis(wav)
wav = window(wav)
Adft, Pdft = fft(wav)
```

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

![Figure_6.png](https://qiita-image-store.s3.amazonaws.com/0/221435/342b84b7-fe8a-1d56-78e1-fcf8f8ab1203.png)

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

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


def melFilterBank(spec, fs=8820.0, nfft=2048, numChannels=20):
    """メルフィルタバンクを作成"""
    # ナイキスト周波数（Hz）
    fmax = fs / 2
    # ナイキスト周波数（mel）
    melmax = hz2mel(fmax)
    # 周波数インデックスの最大数
    nmax = nfft / 2
    # 周波数解像度（周波数インデックス1あたりのHz幅）
    df = fs / nfft
    # メル尺度における各フィルタの中心周波数を求める
    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((int(numChannels), int(nmax)))
    for c in np.arange(0, numChannels):
        # 三角フィルタの左の直線の傾きから点を求める
        increment= 1.0 / (indexcenter[c] - indexstart[c])
        for i in np.arange(indexstart[c], indexcenter[c]):
            i = int(i)
            c = int(c)
            filterbank[c, i] = (i - indexstart[c]) * increment
        # 三角フィルタの右の直線の傾きから点を求める
        decrement = 1.0 / (indexstop[c] - indexcenter[c])
        for i in np.arange(indexcenter[c], indexstop[c]):
            i = int(i)
            filterbank[c, i] = 1.0 - ((i - indexcenter[c]) * decrement)

    mspec = np.log10(np.dot(spec, filterbank.T))
    return mspec, fcenters


wave_data, fs = read_wave(DATA_DIR + "a.wav")
wav = cut_wave(wave_data)
wav = preEmphasis(wav)
wav = window(wav)
Adft, Pdft = fft(wav)
Adft, fcenters = melFilterBank(Adft)
Pdft, fcenters = melFilterBank(Pdft)
```


## 7.離散コサイン変換

![Figure_7.png](https://qiita-image-store.s3.amazonaws.com/0/221435/dee0c6ee-d371-245c-3150-868574c8de69.png)


```python3
def dct(spec, nceps=12):
    """
    離散コサイン変換
    """
    spec = scipy.fftpack.realtransforms.dct(spec, type=2, norm="ortho", axis=-1)
    return spec[:nceps]


wave_data, fs = read_wave(DATA_DIR + "a.wav")
wav = cut_wave(wave_data)
wav = preEmphasis(wav)
wav = window(wav)
A_spec, P_spec = fft(wav)
A_spec, fcenters = melFilterBank(A_spec)
P_spec, fcenters = melFilterBank(P_spec)
A_spec = dct(A_spec)
P_spec = dct(P_spec)
```