MATLABでReverb Pluginを作る！

# はじめに

こちらのツイートが思ったより伸びたので，簡単にまとめておきます．

<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">💯HOW TO MAKE REVERB💯<br><br>空間のコントロールに欠かせない，リバーブサウンドの作り方を解説します😎<a href="https://twitter.com/hashtag/DTMer%E3%81%A8%E7%B9%8B%E3%81%8C%E3%82%8A%E3%81%9F%E3%81%84?src=hash&amp;ref_src=twsrc%5Etfw">#DTMerと繋がりたい</a> <a href="https://twitter.com/hashtag/DTM?src=hash&amp;ref_src=twsrc%5Etfw">#DTM</a> <a href="https://twitter.com/hashtag/FLStudio?src=hash&amp;ref_src=twsrc%5Etfw">#FLStudio</a> <a href="https://twitter.com/hashtag/MATLAB?src=hash&amp;ref_src=twsrc%5Etfw">#MATLAB</a> <br><br>⚠ソースコードはリプ欄に！⚠ <a href="https://t.co/iXzAoljK9O">pic.twitter.com/iXzAoljK9O</a></p>&mdash; 9W3R7Y（くわーてぃ） (@qwerty_16180339) <a href="https://twitter.com/qwerty_16180339/status/1435214089089798149?ref_src=twsrc%5Etfw">September 7, 2021</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>

# MATLABでVST？
VSTというと，[JUCE](https://juce.com/)を用いたり，[C++で1から実装](https://vstcpp.wpblog.jp/)したりすることが思い浮かぶと思います．しかし，実はMATLABでもVSTを作ることができちゃいます！

MATLABでVSTを作るには，[Audio Toolbox](https://jp.mathworks.com/products/audio.html)を用います．Audio Pluginクラスを継承してオリジナルのクラスを作成し，`generateAudioPlugin()`を用いてdllファイルを生成することができます．

UIの設計などに関しては，JUCEやC++での実装と比較して自由度の面で制限があります．したがって，既製品のプラグインのような見た目・操作性を実現することに関しては難しいです．しかし，自分用の簡単なプラグインを作成したり，アルゴリズムのテストをしたりする場合には，非常に強力なツールとなるはずです．

ちなみに，JUCEのプロジェクトを生成することが可能です．なので，MATLABでアルゴリズムを実装し，JUCEでUIを設計するというワークフローも考えられます．

#Schroeder's Reverberator
今回はSchroeder's Reverberator [1]を実装していきます！

## 原理

Schroeder's Reverberator は以下の図のように，**並列のComb Filter**と**直列のAllpass Filter**から構成されています．並列のComb Filterによりリリースの長い残響を生成し，直列のAllpass Filterがリバーブの密度を増やしています．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/397506/c9269c41-8c1d-8f76-b466-f59e51a21205.png)
*[2]を改変

## パラメタとUIの用意
`audioPlugin`クラスを継承してクラスを作成し，パラメタとUIを用意してあげます．ソースコードは以下のようになります．

<details><summary>ソースコード</summary><div>

```MATLAB
classdef Schroeders < audioPlugin
    %#codegen
    
    properties
        % フィルタの個数
        combN = 4
        allN = 2
        
        %フィルタの状態保存用変数
        combZ;
        combidx;
        combm;
        
        allZ;
        allidx;
        allm;
        
        % サンプリング周波数
        Fs = 44100;
        
        
        % 残響時間のパラメタ
        T = 10;
        
        % ステレオ感のパラメタ
        W = 1;
        
        % Dryの音量のパラメタ
        drymix = 1;
        
        % Wetの音量のパラメタ
        wetmix = 1;
    end
    
    % UI
    properties (Constant)
        PluginInterface = audioPluginInterface(...
            audioPluginParameter('T','DisplayName','Time','Mapping',{'pow',2,0,20}),...
            audioPluginParameter('W','DisplayName','Width','Mapping',{'lin',0,1}),...
            audioPluginParameter('drymix','DisplayName','Dry','Mapping',{'lin',0,1}),...
            audioPluginParameter('wetmix','DisplayName','Wet','Mapping',{'lin',0,1})...
            );
    end

% 続く
```

</div></details>

## 変数の初期化
フィルタの初期状態を設定します．ここで，遅延数 $m$の値はComb Filterで30~45[ms]程度，Allpass Filterで5[ms]と1.7[ms]に設定します．ソースコードは以下のようになります．

<details><summary>ソースコード</summary><div>

```MATLAB
    methods
        %% Init
        function p = Schroeders()
            % サンプリング周波数の取得
            p.Fs = getSampleRate(p);
            
            % combフィルタの初期化
            p.combZ = zeros(p.combN,1000,2);
            p.combidx = ones(p.combN,2);
            p.combm = floor((rand(p.combN,2)/2+1)*30*44100/1000);
            
            % allpassフィルタの初期化
            p.allZ = zeros(p.allN,1000,2);
            p.allidx = ones(p.allN,2);
            p.allm = floor([[5;1.7],[5;1.7]]*44100/1000);
        end
%続く
```

</div></details>

## Comb Filterの実装
Comb Filterは入力信号に$m$ サンプルだけ遅らせて出力するものです．この時，出力は $g$ 倍され，フィードバックされます．図にすると以下のようになります．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/397506/4d085ee3-fc34-b363-35fd-3d5f3b06c32f.png)

*[2]より引用

これを実装します．$Z$ にはリングバッファを用いています．

- 入力引数
    - ~：自分自身を表すオブジェクト（Comb Filterでは使わない）
    - X：入力信号（フレームサイズ×チャンネル数 の行列）
    - Z：過去の信号を蓄積したバッファ（適当な長さ×チャンネル数 の行列）
    - i：バッファのインデックス
    - m：遅延数
    - g：フィードバックゲイン

- 出力引数
    - Y：出力信号（フレームサイズ×チャンネル数 の行列）
    - Z：過去の信号を蓄積したバッファ（適当な長さ×チャンネル数 の行列）
    - i：バッファのインデックス
  
ソースコードは以下のようになります．

<details><summary>ソースコード</summary><div>

```MATLAB
        %% comb filter
        function [Y,Z,i] = comb(~, X, Z, i, m, g)

            % 変数の初期化
            nsample = length(X); % サンプル数
            Y = zeros(size(X));  % 出力用の行列
            N = length(Z); % バッファの大きさ

            % インデックスの変換
            ringInd = @(i) mod(i-1+N,N)+1;

            % サンプル毎の処理
            for j = 1:nsample

                % 入力
                x = X(j,:);

                % mサンプル前の値を取得
                z  = Z(ringInd(i-m),:);

                % 出力
                y = z;
                Y(j,:) = y; 

                % 入力のバッファへの格納とフィードバック
                Z(ringInd(i),:) = x + g*z;

                % バッファのインデックスを進める
                i = ringInd(i+1);
            end
        end
%続く
```

</div></details>

作成したComb Filter関数は，遅延ブロックの状態とインデックスの状態を出力し，次回のフレームに引き継ぐことができるようになっています．具体的には，以下のように使用します．

<details><summary>コマンド</summary><div>

```MATLAB
p % プラグイン自身
Z % ゼロで初期化
i % 1で初期化
m % 適当な値に設定
g % 適当な値に設定

% フレーム 1 ------
X %入力信号
[Y,Z,i] = comb(p,X,Z,i,m,g);

% フレーム 2 ------
X %入力信号
[Y,Z,i] = comb(p,X,Z,i,m,g);

% フレーム 3 ------
X %入力信号
[Y,Z,i] = comb(p,X,Z,i,m,g);
```

</div></details>

## Parallel Comb Filter の実装
Schroeder's ReverberatorではComb Filterを並列でつなげます．この時，遅延数 $m$ をそれぞれバラバラにします．また，フィードバックゲイン $g$ は，遅延数 $m$，サンプリング周期$T$，残響時間 $T_r$を用いて以下の式で計算します．この式は「残響音が初期の音量から60dB小さくなるまでの時間」という残響時間の定義から導かれています．
$$
g = 10^{-3mT/T_r}
$$

この時左右のチャンネルそれぞれ別の $m$ の値を用いてステレオ感を出してみました．（追記：[2]を読み進めたら，より効率的な手法も記述されていました）

尚，`p.combN`はCombフィルタの数です．また，`p.combZ`には，各Comb フィルタの遅延ブロックの状態が，`p.combidx`には各Comb フィルタの現在のインデックスが，`p.combm`には各Comb フィルタの遅延数 $m$ の値が保存されています．

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

<details><summary>ソースコード</summary><div>

```MATLAB
        %% Parallel Comb Filter
        function out = paraComb(p,in)
            X = in;
            out = zeros(size(in));
            
            for i = 1:p.combN
                % set parameters
                Z = reshape(p.combZ(i,:,:),[1000,2]);
                idx = p.combidx(i,:);
                m   = p.combm(i,:);
                g   = 10.^(-3*m/p.Fs/p.T);
                
                % filter
                [YL,p.combZ(i,:,1),p.combidx(i,1)] = ...
                    comb(p,X(:,1),Z(:,1),idx(1),m(1),g(1));
                [YR,p.combZ(i,:,2),p.combidx(i,2)] = ...
                    comb(p,X(:,2),Z(:,2),idx(2),m(2),g(2));
                
                % sum
                out = out+[YL,YR];
            end       
        end
%続く
```
</div></details>

## Allpass Filter の実装
Allpass Filterは周波数特性がフラット（"Allpass"という名称の由来です）です．そのため．直列に繋いでもピークが現れません．仕組みはComb Filterをすこしいじったものとなっています．図にすると以下のようになります．
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/397506/70e7cdec-f8de-50ae-f59d-a4d93c68e57d.png)

引数や使用法はComb Filterと全く同じなので説明は省きます．ソースコードは以下のようになります．

<details><summary>ソースコード</summary><div>

```MATLAB
        %% allpass filter
        function [Y,Z,i] = allpass(~, X, Z, i, m, g)

            % 変数の初期化
            nsample = length(X); % サンプル数
            Y = zeros(size(X));  % 出力用の行列
            N = length(Z); % バッファの大きさ

            % インデックスの変換
            ringInd = @(i) mod(i-1+N,N)+1;

            % サンプル毎の処理
            for j = 1:nsample

                % 入力
                x = X(j,:);

                % mサンプル前の値を取得
                z  = Z(ringInd(i-m),:);

                % 出力
                y = (1-g^2)*z - g*x;
                Y(j,:) = y; 

                % 入力のバッファへの格納とフィードバック
                Z(ringInd(i),:) = x + g*z;

                % バッファのインデックスを進める
                i = ringInd(i+1);
            end
        end
%続く
```

</div></details>

## Series Allpass Filter の実装
Schroeder's ReverberatorではAllpass Filterを直列でつなぎます．$g$ の値は一律で0.7程度にします．それ以外はParallel Comb Filterの時と同じです．ソースコードは以下のようになります．

<details><summary>ソースコード</summary><div>

```MATLAB
        %% Series Allpass Filter
        function out = seriAllpass(p,in)
            X = in;
            
            % series
            
            for i = 1:p.allN
                % set parameters
                Z = reshape(p.allZ(i,:,:),[1000,2]);
                idx = p.allidx(i,:);
                m   = p.allm(i);
                g   = 0.7;
                
                % filter
                [YL,p.allZ(i,:,1),p.allidx(i,1)] = allpass(p,X(:,1),Z(:,1),idx(1),m,g);
                [YR,p.allZ(i,:,2),p.allidx(i,2)] = allpass(p,X(:,2),Z(:,2),idx(2),m,g);
                
                X = [YL YR];
            end
            
            out = X;
        end
%続く
``` 

</div></details>

## メインプロセスの実装
最後に，作成した関数を使って，並列Comb Filterと直列のALlpass Filterをつなぎます．また，ステレオ感の調整やDry/Wetを混ぜる処理などを実装します．ソースコードは以下のようになります．

<details><summary>ソースコード</summary><div>

```MATLAB
        %% Process
        function out = process(p,in)
            reb = paraComb(p,in);
            reb = seriAllpass(p,reb);
            
            % Apply Width
            Lreb = reb(:,1);
            Rreb = reb(:,2);
            
            w = p.W/2+1/2;
            
            L = Lreb*w + Rreb*(1-w);
            R = Rreb*w + Lreb*(1-w);
            
            % Mix
            dry = in;
            wet = [L R];
            
            out = p.drymix*dry + p.wetmix*wet;
        end
    end
end
```

</div></details>

実装は以上です．

# ソースコード全体
<details><summary>ソースコード全体を置いておきます．</summary><div>

```MATLAB
classdef Schroeders < audioPlugin
    %#codegen

    properties
        % フィルタの個数
        combN = 4
        allN = 2

        %フィルタの状態保存用変数
        combZ;
        combidx;
        combm;

        allZ;
        allidx;
        allm;

        % サンプリング周波数
        Fs = 44100;


        % 残響時間のパラメタ
        T = 10;

        % ステレオ感のパラメタ
        W = 1;

        % Dryの音量のパラメタ
        drymix = 1;

        % Wetの音量のパラメタ
        wetmix = 1;
    end

    % UI
    properties (Constant)
        PluginInterface = audioPluginInterface(...
            audioPluginParameter('T','DisplayName','Time','Mapping',{'pow',2,0,20}),...
            audioPluginParameter('W','DisplayName','Width','Mapping',{'lin',0,1}),...
            audioPluginParameter('drymix','DisplayName','Dry','Mapping',{'lin',0,1}),...
            audioPluginParameter('wetmix','DisplayName','Wet','Mapping',{'lin',0,1})...
            );
    end

% 続く
    methods
        %% Init
        function p = Schroeders()
            % サンプリング周波数の取得
            p.Fs = getSampleRate(p);

            % combフィルタの初期化
            p.combZ = zeros(p.combN,1000,2);
            p.combidx = ones(p.combN,2);
            p.combm = floor((rand(p.combN,2)/2+1)*30*44100/1000);

            % allpassフィルタの初期化
            p.allZ = zeros(p.allN,1000,2);
            p.allidx = ones(p.allN,2);
            p.allm = floor([[5;1.7],[5;1.7]]*44100/1000);
        end
%続く
        %% comb filter
        function [Y,Z,i] = comb(~, X, Z, i, m, g)

            % 変数の初期化
            nsample = length(X); % サンプル数
            Y = zeros(size(X));  % 出力用の行列
            N = length(Z); % バッファの大きさ

            % インデックスの変換
            ringInd = @(i) mod(i-1+N,N)+1;

            % サンプル毎の処理
            for j = 1:nsample

                % 入力
                x = X(j,:);

                % mサンプル前の値を取得
                z  = Z(ringInd(i-m),:);

                % 出力
                y = z;
                Y(j,:) = y; 

                % 入力のバッファへの格納とフィードバック
                Z(ringInd(i),:) = x + g*z;

                % バッファのインデックスを進める
                i = ringInd(i+1);
            end
        end
%続く
        %% Parallel Comb Filter
        function out = paraComb(p,in)
            X = in;
            out = zeros(size(in));

            for i = 1:p.combN
                % set parameters
                Z = reshape(p.combZ(i,:,:),[1000,2]);
                idx = p.combidx(i,:);
                m   = p.combm(i,:);
                g   = 10.^(-3*m/p.Fs/p.T);

                % filter
                [YL,p.combZ(i,:,1),p.combidx(i,1)] = ...
                    comb(p,X(:,1),Z(:,1),idx(1),m(1),g(1));
                [YR,p.combZ(i,:,2),p.combidx(i,2)] = ...
                    comb(p,X(:,2),Z(:,2),idx(2),m(2),g(2));

                % sum
                out = out+[YL,YR];
            end       
        end
%続く
        %% allpass filter
        function [Y,Z,i] = allpass(~, X, Z, i, m, g)

            % 変数の初期化
            nsample = length(X); % サンプル数
            Y = zeros(size(X));  % 出力用の行列
            N = length(Z); % バッファの大きさ

            % インデックスの変換
            ringInd = @(i) mod(i-1+N,N)+1;

            % サンプル毎の処理
            for j = 1:nsample

                % 入力
                x = X(j,:);

                % mサンプル前の値を取得
                z  = Z(ringInd(i-m),:);

                % 出力
                y = (1-g^2)*z - g*x;
                Y(j,:) = y; 

                % 入力のバッファへの格納とフィードバック
                Z(ringInd(i),:) = x + g*z;

                % バッファのインデックスを進める
                i = ringInd(i+1);
            end
        end
%続く
        %% Series Allpass Filter
        function out = seriAllpass(p,in)
            X = in;

            % series

            for i = 1:p.allN
                % set parameters
                Z = reshape(p.allZ(i,:,:),[1000,2]);
                idx = p.allidx(i,:);
                m   = p.allm(i);
                g   = 0.7;

                % filter
                [YL,p.allZ(i,:,1),p.allidx(i,1)] = allpass(p,X(:,1),Z(:,1),idx(1),m,g);
                [YR,p.allZ(i,:,2),p.allidx(i,2)] = allpass(p,X(:,2),Z(:,2),idx(2),m,g);

                X = [YL YR];
            end

            out = X;
        end
%続く        %% Process
        function out = process(p,in)
            reb = paraComb(p,in);
            reb = seriAllpass(p,reb);

            % Apply Width
            Lreb = reb(:,1);
            Rreb = reb(:,2);

            w = p.W/2+1/2;

            L = Lreb*w + Rreb*(1-w);
            R = Rreb*w + Lreb*(1-w);

            % Mix
            dry = in;
            wet = [L R];

            out = p.drymix*dry + p.wetmix*wet;
        end
    end
end
```

</div></details>

# VSTの生成
上記のソースコードを`Schroeders.m`という名前で保存し，以下のコマンドを実行するとdllファイルを生成できます．dllファイルをDAWで読み込むと，プラグインとして使えます．

```MATALB
>> generateAudioPlugin("Schroeders.m")
```

# 終わりに
お読みいただきありがとうございました．VSTを自分で作れると楽しいので，是非作ってみてください．[この資料](https://www.mathworks.com/content/dam/mathworks/tag-team/Associates/a/89235_PluginGenerationEbriefMWFormat.pdf)とか[MATLABのドキュメント](https://jp.mathworks.com/help/audio/audio-plugin-creation-and-hosting.html)が参考になります．

# 参考文献
[1] Schroeder, M. R. (1962). “Natural Sounding Artificial Reverberation”. J. Audio Eng. Soc., 10(3).
[2] Kahrs, M., & Brandenburg, K. (Eds.). (1998). Applications of digital signal processing to audio and acoustics. Springer Science & Business Media.

# GitHub
僕の作り掛けのプラグインとかが乗ってるので作曲での活用は推奨しませんが，実装の参考程度に参照していただけたら幸いです．
https://github.com/9W3R7Y/MATLAB_VST