wav ファイルから特徴量を求める

**コンペ参加に向けて個人的に調べて用いたものです。**
librosa は知識がなくても各特徴量が簡単に求められます！
知識も経験もないですが，サンプルコードがわりに参照していただければと思います。

# 環境
* Python 3.6.7
* librosa 0.6.2

# 特徴量
ここでは，以下の特徴量の計算手順を記します: 

* stft
* melspectrogram
* melspectrogram delta
* melspectrogram delta2
* mfcc
* mfcc delta
* mfcc delta2

## wav -> stft
短時間フーリエ変換を実行する

```python
import librosa
import os, sys
import pathlib
from tqdm import tqdm
import pickle

# wav ファイルのディレクトリ
dir_raw = pathlib.Path(os.environ['PROJECT_HOME']) / 'raw'
# pickle ファイルの保存先
dir_data = pathlib.Path(os.environ['PROJECT_HOME']) / 'data'

n_fft_sec = 0.025  # ウィンドウサイズ
hop_length_sec = 0.010  # ずらし幅

# train データディレクトリと test データディレクトリそれぞれに実行
for d in ['train', 'test']:
    for _fp in tqdm((dir_raw).glob(f'{d}*/*')):
        # wav ファイルを読み込む (ネイティブサンプリングレートを使用)
        y, sr = librosa.load(_fp, sr=None)

        # stft を実行
        D = librosa.stft(
            y,
            n_fft=librosa.time_to_samples(n_fft_sec,  sr),
            hop_length=librosa.time_to_samples(hop_length_sec,  sr)
        )

        # ディレクトリがなければ生成
        if not (dir_data / d).exists():
            (dir_data / d).mkdir()

        # 結果を pickle で保存
        with (dir_data / d / f'{_fp.name.split(".")[0]}.pickle').open('wb') as f:
            pickle.dump(D, f)
```

## wav -> melspectrogram

```python
import librosa
import librosa.feature
import os, sys
import numpy as np
import pathlib
from tqdm import tqdm
import pickle

# wav ファイルのディレクトリ
dir_raw = pathlib.Path(os.environ['PROJECT_HOME']) / 'raw'
# pickle ファイルの保存先
dir_data = pathlib.Path(os.environ['PROJECT_HOME']) / 'data'

n_fft_sec = 0.025  # ウィンドウサイズ
hop_length_sec = 0.010  # ずらし幅

# train データディレクトリと test データディレクトリそれぞれに実行
for d in ['train', 'test']:
    for _fp in tqdm((dir_raw).glob(f'{d}*/*')):
        # wav ファイルを読み込む (ネイティブサンプリングレートを使用)
        y, sr = librosa.load(_fp, sr=None)

        # メルスペクトログラムを計算
        S = librosa.feature.melspectrogram(
            y,
            n_fft=librosa.time_to_samples(n_fft_sec,  sr),
            hop_length=librosa.time_to_samples(hop_length_sec,  sr),
            sr=sr,
            power=2.0
        )

        # 保存先ディレクトリ名
        d_name = f'mel_spect-{d}'

        # ディレクトリがなければ生成
        if not (dir_data / d_name).exists():
            (dir_data / d_name).mkdir()

        # 結果を pickle で保存
        with (dir_data / d_name / f'{_fp.name.split(".")[0]}.pickle').open('wb') as f:
            pickle.dump(S, f)
```

### グラフ化
それぞれ以下のスクリプトで可視化できる。 (要jupyter)


```python
import librosa.display
import IPython.display
import matplotlib.pyplot as plt

# 上を参考に melspectrogram を求める

print(S.shape)
IPython.display.display(IPython.display.Audio(data=y, rate=sr))
plt.figure(figsize=(10, 4))
librosa.display.specshow(
    librosa.power_to_db(S, ref=np.max),
    y_axis='mel', fmax=8000,
    x_axis='time',
    sr=sr
)
plt.colorbar(format='%+2.0f dB')
plt.title('Mel spectrogram')
plt.tight_layout()
plt.show()
```

## melspectrogram -> delta, delta2

```python
import librosa
import librosa.feature
import os, sys
import numpy as np
import pathlib
from tqdm import tqdm
import pickle

# pickle ファイルの保存先
dir_data = pathlib.Path(os.environ['PROJECT_HOME']) / 'data'

n_fft_sec = 0.025  # ウィンドウサイズ
hop_length_sec = 0.010  # ずらし幅

# train データディレクトリと test データディレクトリそれぞれに実行
for d in ['train', 'test']:
    for _fp in tqdm((dir_data).glob(f'mel_spect-{d}/*')):
        # 上で wav -> melspectrogram の結果を保存したものを読み出す
        with _fp.open('rb') as f:
            S = pickle.load(f)

        # delta* を計算する
        d1 = librosa.feature.delta(S, width=5, order=1)  # delta
        d2 = librosa.feature.delta(S, width=5, order=2)  # delta2

        # 保存先ディレクトリ名
        d1_name = f'del1-{d}'
        d2_name = f'del2-{d}'

        # ディレクトリがなければ生成
        for d_name in [d1_name, d2_name]:
            if not (dir_data / d_name).exists():
                (dir_data / d_name).mkdir()

        # 各結果を pickle で保存する
        with (dir_data / d1_name / _fp.name).open('wb') as f:
            pickle.dump(d1, f)
        with (dir_data / d2_name / _fp.name).open('wb') as f:
            pickle.dump(d2, f)
```


## wav -> mfcc, mfcc_del1, mfcc_del2
librosa では wav から直接それぞれを求めることができる。

```python
import librosa
import librosa.feature
import os, sys
import numpy as np
import pathlib
from tqdm import tqdm
import pickle

# wav ファイルのディレクトリ
dir_raw = pathlib.Path(os.environ['PROJECT_HOME']) / 'raw'
# pickle ファイルの保存先
dir_data = pathlib.Path(os.environ['PROJECT_HOME']) / 'data'

n_fft_sec = 0.025  # ウィンドウサイズ
hop_length_sec = 0.010  # ずらし幅

# train データディレクトリと test データディレクトリそれぞれに実行
for d in ['train', 'test']:
    for _fp in tqdm((dir_raw).glob(f'{d}*/*')):
        # wav ファイルを読み込む (ネイティブサンプリングレートを使用)
        y, sr = librosa.load(_fp, sr=None)

        # mfcc を求める (次元数は13: 13次元が通例(?))
        mfcc = librosa.feature.mfcc(
            y,
            n_fft=librosa.time_to_samples(n_fft_sec,  sr),
            hop_length=librosa.time_to_samples(hop_length_sec,  sr),
            sr=sr,
            power=2.0,
            n_mfcc=13
        )

        # mfcc delta* を計算する
        delta_mfcc  = librosa.feature.delta(mfcc)  # mfcc delta
        delta2_mfcc = librosa.feature.delta(mfcc, order=2)  # mfcc delta2

        # 保存先ディレクトリ名
        d0_name = f'mfcc-{d}'
        d1_name = f'mfcc_del1-{d}'
        d2_name = f'mfcc_del2-{d}'

        # ディレクトリがなければ生成
        for d_name in [d0_name, d1_name, d2_name]:
            if not (dir_data / d_name).exists():
                (dir_data / d_name).mkdir()

        # 各データを pickle で保存
        with (dir_data / d0_name / f'{_fp.name.split(".")[0]}.pickle').open('wb') as f:
            pickle.dump(mfcc, f)
        with (dir_data / d1_name / f'{_fp.name.split(".")[0]}.pickle').open('wb') as f:
            pickle.dump(delta_mfcc, f)
        with (dir_data / d2_name / f'{_fp.name.split(".")[0]}.pickle').open('wb') as f:
            pickle.dump(delta2_mfcc, f)
```

### グラフ化
それぞれ以下のスクリプトで可視化できる。 (要jupyter)

```python
import librosa.display
import IPython.display
import matplotlib.pyplot as plt

# 上を参考に mfcc, mfcc delta, mfcc delta2

print(mfcc.shape, delta_mfcc.shape, delta2_mfcc.shape)

IPython.display.display(IPython.display.Audio(data=y, rate=sr))
plt.figure(figsize=(12, 6))

plt.subplot(3,1,1)
librosa.display.specshow(mfcc)
plt.ylabel('MFCC')
plt.colorbar()

plt.subplot(3,1,2)
librosa.display.specshow(delta_mfcc)
plt.ylabel('MFCC-$\Delta$')
plt.colorbar()

plt.subplot(3,1,3)
librosa.display.specshow(delta2_mfcc, sr=sr, x_axis='time')
plt.ylabel('MFCC-$\Delta^2$')
plt.colorbar()

plt.tight_layout()
plt.show()
```

# まとめ
簡単に各値が求められ，いい感じにグラフに出力できたと思います。
librosa 便利！
あとは頑張ってモデル構築ですね。

また音声ファイルなどを扱う機会があれば改めてまとめたいと思います。

# 参考
* [LibROSA &mdash; librosa 0.6.2 documentation](https://librosa.github.io/librosa/0.6.2/index.html)

# 次に読みたいもの
* [音声分類を色々なモデルや特徴量でやってみた - Qiita](https://qiita.com/kshina76/items/5686923dee2889beba7c)