openSMILEをPythonで使ってみた

## 前書き

- Pythonで使える音響特徴量の解析ライブラリはないかと探してたら，openSMILEがPythonのライブラリに対応したという記事を見かけた
    - [openSMILE3.0がリリースされた](https://qiita.com/Nahuel/items/e8745bfbd710a90a7ca3)（Nafuel氏による投稿）
- あらこれは素晴らしい，ということでNafuel氏の投稿や[公式リファレンス](https://audeering.github.io/opensmile-python/)を見ながら触ってみることにした


## 動作環境（Anacondaの仮想環境として構築）

- Python 3.7.10 
- opensmile 2.2.0
    - これは公式リファレンスに書いてある通り`pip install opensmile`でインストールしたった

※ちなみに，マシンはWindows10だけど特にopenSMILEのbuildとかをしたわけではない

## 使ってみる

### ファイル名を入力としてみた

リファレンスのトップページに「Code example」とあるので，それを参考にした

```python
import opensmile

smile = opensmile.Smile(
    feature_set=opensmile.FeatureSet.ComParE_2016,
    feature_level=opensmile.FeatureLevel.Functionals,
)
feat = smile.process_file('audio.wav')
print(feat.shape)    # [1 rows x 6373 columns]
```

- まず，`smile = opensmile.Smile(...)`でインスタンス化した関数を作っているわけだが，ここで`feature_set`と`feature_level`を指定する必要がある
    - `feature_set`
        - [opensmile.FeatureSet](https://audeering.github.io/opensmile-python/api-smile.html#featureset)を指定することで，どのような特徴量をセットするのかを決めれる
            - `ComParE_2016`（デフォルト値）
                - INTERSPEECH 2016のパラ言語チャレンジのベースラインで使われていた特徴量っぽい
            - `GeMAPS`
            - `GeMAPSv01a`
            - `GeMAPSv01b`
            - `eGeMAPS`
            - `eGeMAPSv01a`
            - `eGeMAPSv01b`
            - `eGeMAPSv02`
            - `emobase`
    - `feature_level`
        - [opensmile.FeatureLevel](https://audeering.github.io/opensmile-python/api-smile.html#featurelevel)を指定することで，特徴量の計算レベルを決めれる
            - `LowLevelDescriptors`: フレーム単位で各種特徴量を計算する，略称はLLD
            - `LowLevelDescriptors_Deltas`: LLDのデルタ特徴量を計算する
                - 現状，`ComParE_2016`のみに対応しているみたい
            - `Functionals`: LLDの統計値を計算する（デフォルト値）
- featにはpandasのDataFrame形式で特徴量が格納されている
    - 公式リファレンスや`print(feat.shape)`の通り，特徴量は6373次元あるみたい
    - 比較的特徴量次元数の少ない`GeMAPSv01b`や`eGeMAPSv02`でもそれぞれ62, 88次元あるみたい
    - 特徴量の命名規則は[ココ](https://audeering.github.io/opensmile/get-started.html#extracting-your-first-features)とか[ココ](https://audeering.github.io/opensmile/about.html)を見ると良さそう
- 機械学習に入れるなら，PCAやLDAみたいな次元圧縮をかけるか，DNNなどに要・不要まで判断させるといいのかも
    - LDDの結果を使うなら次元圧縮よりはDNNに判断させる方がいいかもしれない
- 分析に使うなら，相関係数とか有意確率に対して閾値を設けて，それでいったん重要そうな音響特徴量だけピックアップする方がいいかも


#### LLD使ってみる

```python
import opensmile

smile = opensmile.Smile(
    feature_set=opensmile.FeatureSet.ComParE_2016,
    feature_level=opensmile.FeatureLevel.LowLevelDescriptors,
)
feat = smile.process_file('audio.wav')
print(feat.shape)    # [858 rows x 65 columns]
```

- LLDだと音響特徴量が統計値でまとめられることはないのでrawな値が65次元出てきている
- 音声は8.62秒のものを使ったが，それに対して858フレームの特徴量が抽出されている
- python上だと何をどうしても見にくいと思うので，`feat.to_csv('feat.csv')`とでもしてcsvファイルに書き出してあげるといいかもしれない
    - 特徴量（`feat`）はpandasのDataFrameとして出力されるので，外部ファイル出力もDataFrameのそれと同じ

### numpy配列の信号を入力としてみた

```python
import numpy as np
from scipy.io import wavfile
import opensmile

smile = opensmile.Smile(
    feature_set=opensmile.FeatureSet.ComParE_2016,
    feature_level=opensmile.FeatureLevel.LowLevelDescriptors,
)
rate, signal = wavfile.read('audio.wav')
feat = smile.process_signal(signal / 2**15, rate)
print(feat.shape)    # [858 rows x 65 columns]
```

- 色々と試してみたが，信号を入力とするときは-1から1の範囲に正規化されたnumpyの配列を渡す必要がある
    - openSMILEはその前提で動いているようなので，それ以外（例えばnp.int16型）だとおかしな結果になる
    - `scipy.io.wavfile.read()`はwaveファイルのヘッダ情報を見て変数の型を決めているようで，ここではnp.int16で読み込まれているため2の15乗で正規化している
    - ちなみに，他の読み込み方だと…
        - librosaは正規化された信号を取得できるが，入力にサンプリング周波数を渡してやる必要があるのでちょっと使いにくいかも
        - 標準ライブラリの`wave`を使う場合や，rawファイルを`open()`で読み込む場合だとbytes型で読み込まれるので，`np.frombuffer()`などでnumpyの配列に直す必要がある

#### VAD前と後でなんか違う？

```python
import numpy as np
from scipy.io import wavfile
import webrtcvad
import opensmile

def VAD(signal, rate, mode=0, frame_duration=10):
    vad = webrtcvad.Vad(mode)
    idx = 0
    length = int(rate / 1000 * frame_duration)
    out = []
    while idx<signal.shape[0]:
        if (signal.shape[0]-idx)<length:
            frame = signal[-length: ]
        else:
            frame = signal[idx: idx+length]
        try:
            if vad.is_speech(frame.tobytes(), rate):
                out.append(signal[idx: idx+length])
        except:
            set_trace()
        idx += length
    return np.concatenate(out).astype(np.int16)

smile = opensmile.Smile(
    feature_set=opensmile.FeatureSet.ComParE_2016,
    feature_level=opensmile.FeatureLevel.Functionals,
)
rate, signal = wavfile.read('audio.wav')    # print(signal.shape[0] / rate) -> 8.6210625 (sec.)
vad_signal = VAD(signal)    # print(vad_signal.shape[0] / rate) -> 6.99 (sec.)
feat = smile.process_signal(signal / 2**15, rate)
vad_feat = smile.process_signal(vad_signal / 2 **15, rate)

print(feat["F0final_sma_amean"].values)    # [166.37993]
print(vad_feat["F0final_sma_amean"].values)    # [166.49992]

```

- [webrtcvad](https://github.com/wiseman/py-webrtcvad)で10ミリ秒ごとにVADをかけてみた
    - VADの結果，おおよそ1.63秒ほどが非音声区間として落とされている
- VAD前後で値が異なることから，何かしらの影響はある模様
- 一方で，LLDの結果で単純にF0を平均してみると104.1329となることから，おそらく音声確率あたりを使って内部ではVADらしい挙動が走っているとみなしていいと思う
    - 数値の違いはwebrtcvadとopensmileの音声判定基準の違い，opensmileのアルゴリズムの問題あたりだと思う
        - opensmileの特徴量にはsmaという名前がついているが，これは公式リファレンス曰く「they were smoothed by a moving average filter with window length 3.」とのことである．
        - window length 3は多分，3フレームの窓ということだろう…（そう信じたい，少なくとも3秒はおかしい）
        - 仮に音声確率などに関係なくsmoothingしているとすると，その影響もあるかもしれない


## まとめ

openSMILEをPythonで触ってみた

- ソースからのbuildなど必要なく，pipだけで簡単に導入できる
- openSMILEに実装されている多種多様な音響特徴量を手軽に調べることができる
- 入力も音声ファイルそのもの・numpy配列のどっちでもいける
    - ただしnumpy配列を使う時は，信号強度を-1から1の範囲に正規化しないといけない
- フレーム単位での特徴抽出もできるし，音声ファイル全体で統計的に丸めることもできる
    - 統計的に丸める場合は，VADのようなことを内部でやってくれているみたい
- 今回は1つのファイルを処理しただけだけど，`opensmile.Smile.process_files()`とか複数ファイルをlist型で与えてどんどん処理するとかもいけそう