MATLAB SimulinkでVTuberっぽいことをしてみる

#はじめに
MATLABのバージョンR2020bで，ディープニューラルネットワーク (DNN) のモデルをSimulink上で動かすことのできる`Predict`ブロックが追加されました．今回はこの`Predict`ブロックを用いてリアルタイムに母音（？）を認識し，イラストの口を動かすことでVTuberっぽいことをしてみます！

完成形は以下のような感じです．

https://twitter.com/qwerty_16180339/status/1368939469609394177?s=20

#必要なもの
今回の記事では，音声処理をするための`Audio Toolbox`と，ディープラーニング (DL) 用の`Deep Learning Toolbox`，およびそれらの依存するToolboxが必要となります．また，音声を録音するためのマイクや，録音した音声をトリミングするための音声編集ソフトが必要となります．

#システムの概要
システムの概要は以下のようになります．

1) マイクから音声を受け取ります．音声は一定の長さを持ったフレームとして受け渡されます．その後，フレーム毎にプリエンファシスを適用し，窓を掛け，LPC係数を計算します．このLPC係数が，DNNに入力される特徴量となります．

2) 特徴量を計算したら，それを`Predict`ブロックに入力して，そのフレームの音声がどの母音の音声なのかを予測します．DNNの出力は確率分布の配列になります．

3) DNNの出力のうち，最も値の大きい要素と対応するイラストを表示します．![Diagram.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/397506/f12efedd-ee2a-b956-bf0e-74816a74b50e.png)

それでは，具体的な説明に移ります．

#ディレクトリの作成
今回は教師データの元となる音声データや，教師データ，学習済みモデル，表示するイラストなど，多くのファイルを取り扱います．そのため，あらかじめディレクトリを用意しておきます．

以下のような構造のディレクトリを用意してください．

```
Vtuber:.
├─audio
├─data
├─image
└─net
```

各ディレクトリは以下のように使い分けます．

* `audio`  ：音声データの保存先
* `data`   ：教師データの保存先
* `image`  ：画像データの保存先
* `net`    ：学習済みネットワークの保存先

## 音声データの用意
教師データを作成するために，音声データを録音します．

それぞれの母音（？）を発音し、それを録音します。この時、「あ」「い」「う」「え」「お」「ん」「無音」の計7種類の音を録ります。（「ん」と「無音」は母音じゃないですが、面倒なのでこれ以降は全部母音と呼ぶことにします）

この時、なるべく自然な声で発音し、音程も自然な感じに揺らすと良いです。なにかの文章を全て母音に置き換えて読み上げると自然な感じになります。音声のサンプルは[GitHub](https://github.com/qwerty16180339/Be-like-VTuber-with-MATLAB)からダウンロードしてください（かなりの奇声です）

音声は、5秒程度となるようにトリミングし、サンプリング周波数が44.1 kHzのモノラルの音声として保存します。保存先は`audio`フォルダで、各母音ごとに以下のようなファイル名をつけます。

* あ：`a.wav`
* い：`i.wav`
* う：`u.wav`
* え：`e.wav`
* お：`o.wav`
* ん：`n.wav`
* 無音：`-.wav`

# 特徴量の計算と教師データの作成
音声データを用意したら，特徴量を計算するためのスクリプトを書きます．スクリプトの概要は以下のようになります．

1) `audio`ディレクトリから音声を読み込み，細かいフレームに分けます．フレーム長を512サンプルとし，フレームシフト長を128とします．

2) フレーム毎にプリエンファシスを適用し，窓をかけ，LPC係数を計算します．プリエンファシスとは音声の高音域を強調する処理です．「窓をかける」とは，信号に対してHamming窓などの窓関数を乗算する処理のことです．LPCとは，人間の声道を音響管とみたて，そのパラメータを推定することで，音声を符号化するものです（僕もよくわかってません）．LPC係数とは，推定される音響間のパラメータです．

3) LPC係数の2番目の要素から最後の要素までの値と，元音声の母音のラベルをセル配列に保存します．LPC係数の1番目の要素は常に1なので，これは無視します．

4) すべての処理を終えたら，セル配列を`data`ディレクトリに保存します．

ソースコードは以下のようになります．

```MATLAB:gen_data.m
vowels = ['a', 'i', 'u', 'e', 'o', 'n', '-'];   %音の種類
p = 128;                                        %LPCの次数
N = 512;                                        %フレーム長
nshift = 128;                                   %フレームシフト長
preemph = [1 -0.97];                            %プリエンファシスの係数
w = hamming(N);                                 %窓関数

X = {};                                         %LPC係数保存用のセル配列
Y = {};                                         %ラベル保存用のセル配列
n = 1;                                          %データの番号

%母音毎のループ
for i = 1:length(vowels)
    %音声の読み込み
    [x, Fs] = audioread(append('audio\',vowels(i),'.wav'));
    
    %音声の長さ
    L = length(x);
    
    %フレーム毎のループ
    for t = 1:nshift:L-N
        %フレームの切り出し
        %→プリエンファシス
        %→窓かけ
        frame = filter(preemph,1,x(t:t+N-1)).*w;
        
        %LPC係数計算
        A = lpc(frame,p);
        
        %特徴量をセル配列に保存
        %LPC係数の2番目から最後を取り出す
        %→転置(よくわからないけど転置したらうまくいった)
        X{n,1} = A(2:end).';
        
        %ラベルの保存
        %categoricalな値であると示す
        Y{n,1} = categorical(cellstr(vowels(i)));
        
        %データ番号の更新;
        n = n + 1;
    end
end

%保存
save("data\data",'X',"Y")
```

このスクリプトを一番上の階層に保存し，実行すると，`data`ディレクトリに`data.mat`が保存されます．

#DNNモデルの作成
教師データを準備したら，次はDNNのモデルを作ります．コマンドラインに以下のように入力し，deepNetworkDesignerを開き，「空のネットワーク」を選択してください．

```
>> deepNetworkDesigner
```

そして，以下のようなモデルを組みます．赤字で書いてあるところだけ設定し，あとはデフォルトで大丈夫です．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/397506/6975070f-27b9-44ff-fb81-0e647d6c0fa9.png)


モデルを組んだら，エクスポートボタンを押してモデルをワークスペースに保存してください．すると，ワークスペースに`layers_1`という変数が追加されます．誤って削除すると面倒なので，これをドラッグアンドドロップでディレクトリの最上位の階層に保存し，`layers.mat`という名前に変更してください．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/397506/6ea1b0de-cf8b-8e94-7b61-a3906f734c59.png)

#ネットワークの学習
いよいよネットワークの学習に移ります．

まずは教師データを訓練データ，検証データ，テストデータに三分割します．訓練データは，ネットワークの訓練に，検証データは，ネットワークの検証に用います．テストデータは今回は使いません（Simulink上でリアルタイム評価をするのでいりません）．その後学習に関する設定を行ってから学習を実行し，学習を終えるとネットワークを保存します．

```MATLAB:train.m
%教師データの読み込み
load("data\data.mat")
%未学習のネットワークの読み込み
load("layers.mat")

%インデックスを分割
[trainInd,valInd,testInd] = dividerand(length(X));
%データの分割
Xtrain = X(trainInd);   Xval = X(valInd);   Xtest = X(testInd);
Ytrain = Y(trainInd);   Yval = Y(valInd);   Ytest = Y(testInd);

%設定と学習
miniBatchSize = 27;
options = trainingOptions('adam', ...
    'ExecutionEnvironment','cpu', ...
    'MaxEpochs',100, ...
    'MiniBatchSize',miniBatchSize, ...
    'ValidationData',{Xval,Yval}, ...
    'GradientThreshold',2, ...
    'Shuffle','every-epoch', ...
    'Verbose',false, ...
    'Plots','training-progress');
net = trainNetwork(Xtrain,Ytrain,layers_1,options);

%ネットワークの保存
save("net\net","net");
```
このスクリプトを最上位の改装に保存し，実行すると，以下のようなウィンドウが現れ，学習の進捗が表示されます．学習を終えると，学習済みのモデルが`net`ディレクトリに保存されます．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/397506/adbc6b14-eed2-bac8-dbdb-f8db36b41053.png)

#イラストの用意
表示するイラストを用意します．画力がない場合には神絵師に依頼しましょう．まあ，私は神絵師なので問題ないですが．

母音ごとに口の形を変えた顔を用意して，以下のようなファイル名で`images`ディレクトリに保存します．

* あ：`a.jpg`
* い：`i.jpg`
* う：`u.jpg`
* え：`e.jpg`
* お：`o.jpg`
* ん：`n.jpg`
* 無音：`-.jpg`

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/397506/754a020a-65da-89e8-8594-83cb6fc47cbd.png)

#Simulinkモデルの作成
いよいよ学習したモデルを使って，リアルタイムに画像を表示します．

まずは以下のようにブロックを組みます．`Audio Device Readre`　→　`MATLAB Function` → `Window Function` → `MATLAB Function` → `Predict` →　`MATLAB Function`の順です．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/397506/3976ce4b-b066-509e-ecc3-a6e44b6f2bf1.png)

次に，プリエンファシスの処理を書きます．さっきと全く同じです．コードは以下のようになります．

```matlab
function frame = fcn(frame)
    %プリエンファシス
    preemph = [1 -0.97];
    frame = filter(preemph,1,frame);
```

次に，LPC係数を計算する処理を書きます．関数`lpc()`はコード生成に対応していないらしいのですが，外部関数とすることで使用することができるみたいです．

```matlab
function X = fcn(frame)
    %lpcを外部関数とする
    coder.extrinsic("lpc");
    
    %LPC係数の計算
    %外部関数から得た数値にはいろいろと制約があるっぽかった
    %（例えば配列のスライスが使えなかった）ので，あらかじめ作成した変数に
    %代入することでこれを回避してます (何かいい方法があればぜひ教えてください！）
    N = 128;
    A = zeros(N+1,1);
    A = lpc(frame,N);
    X = A(2:end);
```
次に，Predictブロックの設定を変更します．さっきほど作成・保存した学習済みのネットワークを指定し，ミニバッチのサイズを1にします．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/397506/cf9706a1-ed7b-2d5b-f9d1-651e63caf698.png)
最後にイラストを表示する処理を書きます．受け取った配列のうち，最も値の大きな値と対応する画像を表示するだけです．コードは以下のようになります．

```matlab
function fcn(prob)
    %コード生成に対応しない関数を外部関数とする
    coder.extrinsic("append","image");

    %永続関数imgsに画像を格納
    persistent imgs
    if isempty(imgs)
        %音の種類
        vowels = ['a', 'i', 'u', 'e', 'o', 'n', '-'];

        %セル配列imgsを作成
        imgs = cell(length(vowels),1);

        %音の種類ごとに画像を読み込み
        for i = 1:length(vowels)
            %画像を格納する配列を作成
            %注：画像のサイズと合わせる．チャンネル数は基本的に3
            img = uint8(zeros(160,160,3));

            %画像の読み込みと保存
            img = imread(append('images\',vowels(i),".jpg"));
            imgs{i} = img;
        end
    end

    %値の最も大きい要素のインデックスを取得
    [~,index] = max(prob);

    %取得したインデックスと対応する画像を描画
    image(imgs{index})
    drawnow
```
以上で完成です！

#実行してみる
これを実行すると，こんな感じで，そこそこ認識できていることがわかります．
<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">Qiita記事用のサンプル（Twitterの使い道を間違えている） <a href="https://t.co/rHNNYEg7XJ">pic.twitter.com/rHNNYEg7XJ</a></p>&mdash; 9W3R7Y@ツイート多め (@qwerty_16180339) <a href="https://twitter.com/qwerty_16180339/status/1368930589496111107?ref_src=twsrc%5Etfw">March 8, 2021</a>
</blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

#おまけ：エラーをごまかす
上の動画を見ると，1フレームだけ認識を誤ることがあることがわかります．これをごまかすためには，移動平均フィルタが有効です．以下のようにブロックを組み，移動平均フィルタを適応します．移動平均フィルタのWindow Lengthは2~3がおすすめです．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/397506/a7256e7f-b3a8-6a68-6654-dc48097727fd.png)
すると，先ほどよりも安定感が増します．（ちょっと「お」の精度が低いですが，データが自然な声じゃなかったことが原因だと考えられます．なるべく自然なデータを取るようにしないと，しゃべった時の音声に対応できなくなります．）
<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">Moving Averageでちょっと安定する <a href="https://t.co/iALkmhBtqB">pic.twitter.com/iALkmhBtqB</a></p>&mdash; 9W3R7Y@ツイート多め (@qwerty_16180339) <a href="https://twitter.com/qwerty_16180339/status/1368933395904950276?ref_src=twsrc%5Etfw">March 8, 2021</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

#おわりに
今回はSimulinkの新ブロック`Predict`を使って，VTuberみたいに口を動かすやつを作りました．リアルタイムで動かす環境をあっという間に作れてしまうところが本当に便利だと感じました．また何か作ったら記事書こうと思います．

#参考文献
* [やさしく音声分析法を学ぶ： ケプストラム分析とLPC分析](https://www.slideshare.net/ShinnosukeTakamichi/lpc-49065650)
* [MATLABによる機械学習の基礎　～特徴抽出・分類器・交差検定～](https://jp.mathworks.com/videos/matlab-for-machine-learning-for-beginners-feature-extraction-classifier-and-cross-validation-117948.html)

#GitHub
以下のリポジトリのv1.0のリリースを参照してください．

https://github.com/9W3R7Y/Be-like-VTuber-with-MATLAB/releases/tag/1.0