#03 Vocoder編―MATLAB & Simulinkでオーディオエフェクト

#はじめに
今回はVocoderを作る．いきなりハードルが上がるが，なるべくわかりやすく説明する．

#背景知識
まず，Vocoderについて説明する．

##Vocoderとは何か
Vocoderでは二つの入力を合成する．例えば，声とシンセサイザーの音を合成する．以下の動画が分かりやすい．
https://www.youtube.com/watch?v=Inb8sUgpdVk

##Vocoderの仕組み
[Logic Pro Xのマニュアル](https://support.apple.com/kb/PH27761?viewlocale=ja_JP&locale=ja_JP)にVocoderの仕組みが細かく書いてあったので，それを基に説明する．

入力は分析オーディオソースと合成オーディオソースの二つである．分析オーディオソースには主に人間の声が入力され，合成オーディオソースには主にシンセサイザーの音が入力される．まず，それぞれの入力はフィルタバンクに入力される．フィルタバンクは複数のバンドパスフィルタから構成されており，入力を帯域ごとの音声に分ける．例えば，入力を三つの帯域に分けるフィルタバンクでは，低音域，中音域，高音域の三つに分けられる．次に，分析オーディオソースから，帯域毎のエンヴェロープが抽出される．エンヴェロープとは，音量の推移である．その後，合成オーディオソースのフィルタバンクの出力に対して，分析オーディオソースから抽出された帯域ごとのエンヴェロープが適応される．エンヴェロープ適応後の各帯域の音声をすべて足すと，分析オーディオソースと合成オーディオソースが合成された音声が出力される．

#実装
早速実装に移る．

##フィルタバンクの作成
はじめにフィルターバンクを作成する．

まずは以下のようにブロックを組む．フィルタバンクの出力はそのままスピーカーから出力できないので，代わりにSpectrum Analyzerで可視化する．Colored Noiseではダブルクリックで開く設定画面からPink NoiseをWhite Noiseに変更し，Outoput sample timeは1/44.1/10^3に設定する (サンプリング周波数の逆数である)．MATLAB Function ブロックはFilter Bankに名前を変更しておく．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/397506/1492460d-9e1b-1840-c1c2-2a6d037dc2da.png)


次にFilter Bank ブロックのコードを書く．MATLABのAudio Toolboxには`octaveFilter()`という関数がるので，これを利用する（Filter bankについて教えてくださった[Takacie
](https://twitter.com/Takacie815)さん，ありがとうございました!). オクターブフィルターでは1オクターブ事に音域を分割するが，引数で`1/3 octave`や`1/6 octave`などと指定することで，分数オクターブごとに分割することもできる．分割数が多ければ多いほど，声が聞きやすくなるが，その分重くなる．まずは`1 octave`に設定する．尚，フィルターは過去の音声も参照して音声を処理するので，永続のオブジェクトとする必要がある．以下がフィルターバンクのコードである．

```matlab
function y = fcn(u)
    %フィルタバンクオブジェクトの作成
    persistent ofb
    if isempty(ofb)
        ofb = octaveFilterBank('1 octave');
    end

    %フィルタバンクの適用
    y = ofb(u);
```
これを実行し，Spectrum Analyzerを確認すると，少し荒いが帯域ごとに音を分割できていることが分かる．高音域の方が綺麗に分かれており，一番右はオレンジ，次は青，その次が黄色，というように音が分割できている．

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/397506/96f447ec-e4ae-52b1-710c-ced6b4481f69.png)

##エンヴェロープの実装
次にエンヴェロープを実装する．今回は簡単な方法でエンヴェロープを実装する．ある時刻での入力の値を $x_t$，ある時刻でのエンヴェロープの値を $y_t$ ，$0$ 以上 $1$ 以下の定数を $k$ とすると以下のような処理になる．

```math
\left\{
\begin{array}{ll}
y_t \leftarrow |x_t| \quad (y_{t-1} \leqq x_t) \\
y_t \leftarrow k\times y_t \quad (otherwise)
\end{array}
\right.
```

要するに，入力が現在のエンヴェロープの値よりも大きいときは入力の値の絶対値をそのままエンヴェロープの値とし，小ささいときはゆっくりと小さくする処理である．要素数が帯域分割数となるような永続な配列`y`を作成し，各帯域ごとに入力の全サンプルにわたってこの処理を行えば，エンヴェロープが抽出できる．

以下のようにブロックを組む．入力をWhite NoiseからAudio Device Readerに変更する．マイクが無い場合はいつものようにFrom Multimedia Fileにしてもよい．Filter Bankの隣にMATLAB Functionブロックを追加し，それをEnvelopeと名付ける．Spectrum AnalizerはTime Scopeへと変更する．MATLAB Functionブロックを以下のように書き換える．

```matlab
function envelope = fcn(x,k)
    %入力の配列のサイズから帯域数n_bandsとframe_sizeを求める
    [n_bands,frame_size]=size(x);

    %現在のエンヴェロープの値を保存する永続な配列yの作成
    persistent y
    if isempty(y)
        y = zeros(n_bands);
    end

    %入力の絶対値を取る
    x = abs(x);

    %エンヴェロープの配列の初期化
    envelope = zeros(n_bands,frame_size);
    
    %全体域にわたる繰り返し
    for band = 1:n_bands
        %全サンプルにわたる繰り返し
        for i = 1:frame_size
            %yの値の更新
            if y(band) <= x(band,i)
                y(band) = x(band,i);
            else
                y(band) = y(band)*k;
            end
            envelope(band,i) = y(band);
        end
    end
end
```
上記のコードを書き，保存するとMATLAB Functionブロックに入力`k`が追加されるので，Constantブロックから値を送る．今回はkの値を`0.8`とする．実行し，マイクに向かってしゃべりながらTime Scopeを確認すると，以下のようにエンヴェロープが抽出できていることが分かる．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/397506/d90f346c-9bf4-e46b-f8c5-8d3a07349deb.png)

しかし，現段階ではエンヴェロープが線として認識できないレベルでなめらかではない．これをなめらかにするために，移動平均フィルタをかける．移動平均フィルタにより，高周波数の信号が遮断されるので，なめらかになる．以下のようにブロックを組み，Moving Averageの`Window length`を500に設定する．多ければ多いほどなめらかになるが，その分遅延が生じる．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/397506/91ef9b6a-eba1-1d2a-8afb-2675472ca192.png)

これを実行すると，以下のように，先ほどよりもエンヴェロープがなめらかになっているのが分かる．

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/397506/321a69fd-62c3-1d9c-a8b9-d00bbf002216.png)

##合成
いよいよ二つの音声の合成処理を実装する．二つの音声の合成の際には，一方のフィルタリングの出力を，もう一方のフィルターバンクの出力のエンヴェロープと乗算した後，すべての帯域の音声を加算すればよい．これはSimulinkのブロックで実現できる．以下のようにブロックを組む．先ほどのブロックにAudio Oscillatorを追加し，それをFilter Bankのコピーと接続する．Audio Oscillatorは150HzのSawtoothに設定し，Samples per frameは1024に設定する．その後，Moving Averageブロックの出力とAudio Oscillator側のFilter Bankの出力をProductブロックに入力する．Productブロックの出力は，`加算範囲`を`指定した次元`に設定し，`次元`を2に設定したSum of elementブロックに入力し，スピーカーから出力する．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/397506/6a178761-9f7b-2d38-f65b-04c574e54ab9.png)
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/397506/67955946-1108-9aa2-89e1-836d21482e35.png)

これを実行し，マイクに向かってしゃべると，何かしゃべろうとしているように聞こえる．しかし，何を言っているのかは聞き取れない．

<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">octave filterだと分割数が足りなくて，何て言ってるのかわからない． <a href="https://t.co/spjdam0eLL">pic.twitter.com/spjdam0eLL</a></p>&mdash; 9W3R7Y (@qwerty_16180339) <a href="https://twitter.com/qwerty_16180339/status/1340784437563334656?ref_src=twsrc%5Etfw">December 20, 2020</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

これは，フィルタバンクの分割数が少ないからである．試しに，2つのフィルタバンクのオクターブフィルタバンクを`1/6 octave`に書き換える．コードは以下のようになる．

```matlab
function y = fcn(u)
    %フィルタバンクオブジェクトの作成
    persistent ofb
    if isempty(ofb)
        ofb = octaveFilterBank('1/6 octave');
    end

    %フィルタバンクの適用
    y = ofb(u);
```
すると，「MATLABでVocoder」と言っていることが分かるようになった．
<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">1/6 octave filterだと聞こえる <a href="https://t.co/RSv1DejOpi">pic.twitter.com/RSv1DejOpi</a></p>&mdash; 9W3R7Y (@qwerty_16180339) <a href="https://twitter.com/qwerty_16180339/status/1340784577514622976?ref_src=twsrc%5Etfw">December 20, 2020</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

##フォルマントシフト
ここで，各帯域のエンヴェロープを二次元的に可視化する．現時点のエンヴェロープは時間軸を持つ三次元の情報なので，時間軸で和を取って二次元の情報にし，Array Protで可視化する．以下のようにブロックを組む．Moving Averageブロックの直後のSum of elementブロックの次元は1に設定する．また，1×バンド数の行列となっているものをバンド数×1の行列にするためにTransposeブロックを置き，転置する．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/397506/46e4ba45-9523-cfc9-b65f-269443db1f47.png)

これを実行すると，声の各帯域のエンヴェロープが以下のように可視化される．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/397506/2d2a6511-26b2-88af-f62d-1ca2d0376ad3.png)
横軸がエンヴェロープのインデックスで，縦軸がエンヴェロープの値である．多くのVocoderでは，この山の形を丸ごと左右に移動させるフォルマントシフトの機能が備わっている．本記事でも，フォルマントシフトを実装する．

まずは整数個の要素だけ左右に移動する処理を考える．フォルマントシフトの量を`s`とすると，以下の図のように，フォルマントシフト後の配列の`i`番目の要素には，`i-s`番目の要素が来ることが分かる．尚，存在しない帯域のエンヴェロープの値は0と仮定する．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/397506/d272b900-2433-a5a4-8e02-f77141bfe352.png)

次に少数個の要素の移動を考える．当り前だが，配列のインデックスは整数しかない．したがって，少数個の移動を整数のインデックスで表現する必要がある．本記事では，`s`の整数部分番目の要素と，その次の要素を(`s`の少数部分－1)：`s`の少数部分の比で混ぜることで，これを実現する．例えば，1.3要素の移動では，1番目の要素と2番目の要素を0.7：0.3の割合で混ぜれることで再現する．

以下のようにブロックを組む．Array Plotブロックにフォルマントシフト前の要素とフォルマントシフト後の要素を入力し，両者を比べる．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/397506/dd005312-c15a-ff5f-68eb-899e843bed98.png)

ブロックを設置したら，フォルマントシフトのブロックを以下のように書き換える．shiftの値にフィルターの種類の分母を掛けることで単位をオクターブにしたり，行列の左右を0で埋めて余白を取ることで，上記の処理を実装している．

```matlab
function y  = fcn(u,shift)
    %フレームサイズと帯域の分割数を取得
    [frame_size,n_bands]=size(u);
    
    %エンヴェロープの左右に余白を取る（エンヴェロープの配列と同じ大きさのゼロ行列を結合）
    v = cat(2,cat(2,zeros(frame_size,n_bands),u),zeros(frame_size,n_bands));
    
    %フォルマントシフト後の配列の始めの要素が，フォルマントシフト前のどの要素だったのかを取得
    %余白を取った分，n_bandsだけ平行移動している
    %オクターブフィルタが1/6 octaveなので，6を掛けて単位をオクターブにそろえる
    %マトラボでは始めの要素は１なので，１を足す
    p = n_bands-shift*6+1;
    
    %整数部分の計算
    p_int = floor(p);
    
    %小数部分の計算
    p_dec = p - p_int;
    
    %二つの配列のブレンド
    y = v(:,(p_int):(p_int+n_bands-1))*(1-p_dec)+v(:,(p_int+1):(p_int+n_bands))*p_dec;
```
コードを書き換えたら，以下のようにブロックを組み，フォルマントシフトの量をKnobで操作できるようにする．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/397506/ad7bc4d5-41c9-d960-934b-512464518669.png)
これを実行しKnobをいじると，フォルマントシフトによりエンヴェロープが左右に動くことが分かる．また，音の質感も変わることが分かる．
<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">フォルマントシフト！！！（なんでかよくわかんないけど軽量化した！） <a href="https://t.co/b94Bxmnhjm">pic.twitter.com/b94Bxmnhjm</a></p>&mdash; 9W3R7Y (@qwerty_16180339) <a href="https://twitter.com/qwerty_16180339/status/1340811388378234880?ref_src=twsrc%5Etfw">December 21, 2020</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>
以上でVcoderの実装は終わりである．

#おまけ
ボコーダーでは和音のシンセを用いることが多い．DAW上で和音のシンセの音を作成し，取り込むと以下のように，いかにもボコーダーらしい音になる．
<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">やっぱりボコーダーは和音っしょ <a href="https://t.co/OPixyamNhW">pic.twitter.com/OPixyamNhW</a></p>&mdash; 9W3R7Y (@qwerty_16180339) <a href="https://twitter.com/qwerty_16180339/status/1340815306273349632?ref_src=twsrc%5Etfw">December 21, 2020</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

バンド数は多ければ多いほどきれいな音になる．試しにオクターブフィルタを`1/48 octave`に変更してみる．バンド数は480バンドになる（
<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">マトラボで強引に1/48 octave フィルタで480バンドのVocoderを作ってみたww リアルタイムでは重すぎて無理なので，オーディオ書き出し！めっちゃ綺麗なんだけどww <a href="https://t.co/KUv7G1ISg4">pic.twitter.com/KUv7G1ISg4</a></p>&mdash; 9W3R7Y (@qwerty_16180339) <a href="https://twitter.com/qwerty_16180339/status/1340821376354021376?ref_src=twsrc%5Etfw">December 21, 2020</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>
とてもきれいな音になる．

今度はフィルタの次数を100にしてみる．フィルターバンクを以下のように書き換える

```matlab
function y = fcn(u)
    %フィルタバンクオブジェクトの作成
    persistent ofb
    if isempty(ofb)
        ofb = octaveFilterBank('1/48 octave',"FilterOrder",100);
    end

    %フィルタバンクの適用
    y = ofb(u);
```
すると，遅延の影響で訳の分からない音ができる．これはこれで面白い．
<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">フィルタの次数を100にしてみた．遅延がえぐいw Padの素材になりそう <a href="https://t.co/EKbS2x6uR1">pic.twitter.com/EKbS2x6uR1</a></p>&mdash; 9W3R7Y (@qwerty_16180339) <a href="https://twitter.com/qwerty_16180339/status/1340828319466418177?ref_src=twsrc%5Etfw">December 21, 2020</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

Delayの時と同様に以下のようにブロックを組めば，ステレオの音声を入力できる．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/397506/bb454360-e3a2-ee0c-d3f5-5a2c3f73e3e1.png)
<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">careerをステレオ対応にした！重すぎる <a href="https://t.co/ge2Eqi0JA8">pic.twitter.com/ge2Eqi0JA8</a></p>&mdash; 9W3R7Y (@qwerty_16180339) <a href="https://twitter.com/qwerty_16180339/status/1340842663545729024?ref_src=twsrc%5Etfw">December 21, 2020</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

#追記：フィルタの次数について
実は，デフォルトのフィルタの次数は2となっているが，いろいろ試したところ4から6程度がちょうどいいことが分かった．したがって，フィルタバンクを以下のように書き換えることを推奨する．

```matlab
function y = fcn(u)
    %フィルタバンクオブジェクトの作成
    persistent ofb
    if isempty(ofb)
        ofb = octaveFilterBank('1/48 octave',"FilterOrder",4);
    end

    %フィルタバンクの適用
    y = ofb(u);
```

#おわりに
今回はボコーダーを実装した．自由課題としてのオーディオエフェクトの実装はこれで最後にする．時間があるときに随時更新する．

#GitHub
以下で完成したオーディオエフェクトを公開している．
https://github.com/qwerty16180339/matlab_audio_effects

#参考文献

* Apple, Logic Pro X: ボコーダーの基礎 https://support.apple.com/kb/PH27761?viewlocale=ja_JP&locale=ja_JP
* MathWorks, octaveFilterBank, https://jp.mathworks.com/help/audio/ref/octavefilterbank-system-object.html?searchHighlight=octaveFilterBank&s_tid=srchtitle

※エンヴェロープに関しては父親が話していたことをそのまま実装した．