MATLABでReverb Pluginを作る！

はじめに

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

💯HOW TO MAKE REVERB💯

空間のコントロールに欠かせない，リバーブサウンドの作り方を解説します😎#DTMerと繋がりたい #DTM #FLStudio #MATLAB

⚠ソースコードはリプ欄に！⚠ pic.twitter.com/iXzAoljK9O

— 9W3R7Y（くわーてぃ） (@qwerty_16180339) September 7, 2021

MATLABでVST？

VSTというと，JUCEを用いたり，C++で1から実装したりすることが思い浮かぶと思います．しかし，実はMATLABでもVSTを作ることができちゃいます！

MATLABでVSTを作るには，Audio Toolboxを用います．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がリバーブの密度を増やしています．

*[2]を改変

パラメタとUIの用意

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

ソースコード

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

% 続く

変数の初期化

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

ソースコード

    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の実装

Comb Filterは入力信号に$m$ サンプルだけ遅らせて出力するものです．この時，出力は $g$ 倍され，フィードバックされます．図にすると以下のようになります．

*[2]より引用

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

• 入力引数

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

• 出力引数

  • Y：出力信号（フレームサイズ×チャンネル数 の行列）
  • Z：過去の信号を蓄積したバッファ（適当な長さ×チャンネル数 の行列）
  • i：バッファのインデックス

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

ソースコード

        %% 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
%続く

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

コマンド

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);

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$ の値が保存されています．

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

ソースコード

        %% 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 の実装

Allpass Filterは周波数特性がフラット（"Allpass"という名称の由来です）です．そのため．直列に繋いでもピークが現れません．仕組みはComb Filterをすこしいじったものとなっています．図にすると以下のようになります．

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

ソースコード

        %% 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 の実装

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

ソースコード

        %% 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
%続く

メインプロセスの実装

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

ソースコード

        %% 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

実装は以上です．

ソースコード全体

ソースコード全体を置いておきます．

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

VSTの生成

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

>> generateAudioPlugin("Schroeders.m")

終わりに

お読みいただきありがとうございました．VSTを自分で作れると楽しいので，是非作ってみてください．この資料とかMATLABのドキュメントが参考になります．

参考文献

[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