はじめに
こちらのツイートが思ったより伸びたので，簡単にまとめておきます．
💯HOW TO MAKE REVERB💯
— 9W3R7Y（くわーてぃ） (@qwerty_16180339) September 7, 2021
空間のコントロールに欠かせない，リバーブサウンドの作り方を解説します😎#DTMerと繋がりたい #DTM #FLStudio #MATLAB
⚠ソースコードはリプ欄に！⚠ pic.twitter.com/iXzAoljK9O
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