MATLABでIR(Impulse Response) Speaker Simulator VST Plugin 2

#MATLABの機能を使ってImpulse Response取得～VSTプラグイン作成
Impulse Response(IR)から応答特性のシステム同定を行って、同等の特性のエフェクトを作ってみる話の第2弾。
前の記事

##取得したImpulse Responseの調整

取得したIR特性はゲインもレイテンシもまちまちなので、全てのデータを揃える。プラス、IRデータは44100Hzx4secもあってフィルタにするには巨大すぎるので、無信号部分の不要なデータを削除して可能な限り短く（フィルタのTap数を短く）する。
調整したIRデータをFIRフィルタの係数として使用することで、スピーカーシミュレーターとなる。

###IR調整用コード
以下がデータ調整用のMATLABコード

IRデータ調整

%% IRデータ読み込み、レイテンシ調整、FIRフィルタ化

IRFileName = 'measured_ir_data_AmpsSpeaker.mat';
% IRFileName = 'measured_ir_data_MicroCAB.mat';
load(IRFileName)
IRAmpThreshold = 0.005;     % この振幅を超えた時点のIRを使用する
IRBuffer = 3;               % IRAmpThresholdを超えて、IRBufferサンプル前のデータから採用する
sz = size(measurementData); % データ数　N x 14
NormalizedAmp = 1;  % ノーマライズする振幅
Fs = measurementData.SampleRate(1);
%% Cut data
data = cell(1,sz(1));

figure, hold on
for n = 1:sz(1)
    data{n}.Name = measurementData.Properties.RowNames{n};
    data{n}.AmpTemp = measurementData.ImpulseResponse(n).Amplitude;
    % TimeTemp = measurementData.ImpulseResponse(n).Time;
    data{n}.AmpAbs = abs(data{n}.AmpTemp);
    [data{n}.pks, data{n}.locs] = findpeaks(data{n}.AmpAbs, 'MinPeakHeight',max(data{n}.AmpAbs*0.9),...
        'NPeaks', 1);
    
    % 振幅ノーマライズ
    if data{n}.AmpTemp(data{n}.locs) < 0
        adjustCoef = -1*NormalizedAmp/data{n}.pks;
    else
        adjustCoef = NormalizedAmp/data{n}.pks;
    end
    data{n}.AmpTemp = adjustCoef.* data{n}.AmpTemp;
    data{n}.AmpAbs = abs(data{n}.AmpTemp);
    % エンベロープ Upperサイドのみ
    data{n}.AmpAbsEnv = envelope(data{n}.AmpTemp, 10, 'Peak');
    
    % Peak検出
    % IRAmpThresholdの少なくとも100倍の振幅
    % 数は1つだけ
    % Peakの前と後とで別に処理する
    data{n}.indAmpBeforePeak = find(data{n}.AmpAbs(1:data{n}.locs) >= IRAmpThreshold);
    data{n}.indAmpAfterPeak = find(data{n}.AmpAbs(data{n}.locs:end) >= IRAmpThreshold)+data{n}.locs;
    
    % Threshold部分だけ切り出し、振幅ノーマライズ
    % 位相が反転（IRのPeakが下から始まる）している場合は正転させる
    data{n}.Amp = data{n}.AmpTemp(data{n}.indAmpBeforePeak(1):data{n}.indAmpAfterPeak(end));
    % Peak点からレイテンシを求める
    data{n}.Latency = data{n}.locs-data{n}.indAmpBeforePeak(1);
    % Time = TimeTemp(indAmpBeforePeak(1):indAmpAfterPeak(1));
    % 最初と最後のデータは0補正
    % data{n}.Amp = [0; data{n}.Amp(1:end); 0];
    data{n}.NTap = length(data{n}.Amp);
    data{n}.Time = [0:data{n}.NTap-1]/Fs;
    
    % Plot
    subplot(3,1,1)
    plot(data{n}.Time, data{n}.Amp), hold on, grid on
    title('Impulse Response')
    
    subplot(3,1,2)
    [h, f] = freqz(data{n}.Amp,1, 1024*2,  Fs);
    semilogx(f, 20*log10(abs(h))), hold on, grid on
    title('Magnitude Response')
    xlim([40 Fs/2]), ylim([-40 25])
    
    subplot(3,1,3)
    semilogx(f, angle(h)*180/pi), hold on, grid on
    title('Phase Response')
    xlim([40 Fs/2])
end

subplot(3,1,1)
legend(measurementData.Properties.RowNames{1:end})
set(gcf, 'Position', [10 190 870 790]), shg

###IRの調整結果
元のデータは176400点（44100Hzx4sec）あったが、これをおおよそ1700点ぐらいのデータに削減することができた。FIRフィルタのTap数1700点ぐらいであれば実現可能なレベルになっている。

調整したIRの特性が以下

##Amp+Speakerのリスニングテスト
得られたIR特性を、実際にエフェクトとしてギター音に適用する。
適用方法は簡単で、ストリーミング用のFIRフィルタ関数に入れるだけ。

StreamingFilter

hFilter = dsp.FIRFilter('NumeratorSource', 'Input port');
out = step(hFilter, in, Coeff)

MATLABでVST Pluginを読み込めるので、それでフリーのギターアンプシミュレータVST Plugin 「LeXtac」を使って、その後にスピーカーシミュレータを接続してテストする。
VSTプラグイン読み込みは以下で読み込める。parameterTuner関数は、パラメータ設定用のGUIを自動的に作成して起動する機能。

ExternalVSTImport

Hamp = loadAudioPlugin(which('LeXtac.dll'));
parameterTuner(Hamp)
out = process(Hamp, input);

こんなGUI（画面左）を作成して、複数のスピーカーを切り替えながら実験してみた。
（右側のGUIはparameterTuner）

まず、Thruにすると、ギター音を歪ませただけなのでシャーシャーした音。MicroCABは元々シミュレータなので、どこかウソ臭い。VHT2902の特性が一番気に入った。ボトムがきちんと出ていてかつ引き締まっている。今回2発入りキャビネットだったが、4発入りのキャビネットだともっと低音はしっかり出るんだろう。

すごく簡単にスピーカーシミュレーターが作成できてしまった。
VR/3D AudioやIR Reverbなんかも同じようにして作成できそう。
今度トンネルに機材持って行って、特性取りしてReverb作ってみたい・・・

全プログラムはこちら。

clear all, close all force
%%
Fs = 44100;
% Fs = 32000;
Nch = 1;
Spf = 128;      % 音切れが発生したり、レイテンシが大きい場合はこれを調整。
% Spf = 256;
audioSource1 = 'hum_solo1.wav';
audioSource2 = 'hum_backing2.wav';
load('dataAmpSpeaker.mat');  % 1 2 4 5 6のデータを使用する
load('dataMicroCab.mat');    % 3 4 5 6 7のデータを使用する
data{1} = dataAmpSpeaker{1};
data{2} = dataAmpSpeaker{2};
data{3} = dataAmpSpeaker{4};
data{4} = dataAmpSpeaker{5};
data{5} = dataAmpSpeaker{6};
data{6} = dataMicroCab{3};
data{7} = dataMicroCab{4};
data{8} = dataMicroCab{5};
data{9} = dataMicroCab{6};


%% Create Stop button
screenSize = get(0,'ScreenSize');
figSize = [380 440];        % [X Y]
fg1 = figure('MenuBar','none','Toolbar','none',...
    'Name', 'Effect Parameter', 'NumberTitle', 'off',...
    'Position',[10 screenSize(4)-figSize(2)-40 figSize(1) figSize(2)]);
stpBtn = true;

% Stop button
uicontrol(fg1, 'Style', 'pushbutton', 'String', 'Stop',...
    'Position', [20 20 100 40],'Callback', 'stpBtn = false;');

% Input Select button
sourceSel = 0;
bg = uibuttongroup('Visible', 'on','Position',[0.08 0.40 .3 .30],...
    'SelectionChangedFcn', 'rdBtn = bg.SelectedObject.String');
uicontrol(bg, 'Style', 'radiobutton', 'String', 'External In',...
    'Position', [20 73 100 40], 'Callback', 'sourceSel=0;');
uicontrol(bg, 'Style', 'radiobutton', 'String', 'Source1',...
    'Position', [20 38 100 40], 'Callback', 'sourceSel=1;');
uicontrol(bg, 'Style', 'radiobutton', 'String', 'Source2',...
    'Position', [20 3 100 40], 'Callback', 'sourceSel=2;');

% Cabinet Select button
cabSel = 0;
bg = uibuttongroup('Visible', 'on','Position',[0.44 0.20 .5 .74],...
    'SelectionChangedFcn', 'rdBtn = bg.SelectedObject.String');
Space = 30;, initialPos = 280;
uicontrol(bg, 'Style', 'radiobutton', 'String', 'Thru',...
    'Position', [5 initialPos 180 40], 'Callback', 'cabSel=0;');
for n = 1:numel(data)
    uicontrol(bg, 'Style', 'radiobutton', 'String', data{n}.Name,...
        'Position', [5 initialPos-(n*Space) 180 40], 'Callback', ['cabSel=' num2str(n) ';']);
end

% View Select button
viewTSel = 0;
uicontrol(fg1, 'Style', 'pushbutton', 'String', 'TScope',...
    'Position', [140 20 60 40],'Callback', 'if viewTSel == 0; viewTSel = 1; else, viewTSel = 0; end;');
viewFSel = 0;
uicontrol(fg1, 'Style', 'pushbutton', 'String', 'FScope',...
    'Position', [210 20 60 40],'Callback', 'if viewFSel == 0; viewFSel = 1; else, viewFSel = 0; end;');

%% オーディオソースと再生用オブジェクト
haIn = audioDeviceReader('SampleRate', Fs, 'NumChannels', Nch,...
    'Driver', 'ASIO','SamplesPerFrame', Spf, 'Device', 'M-Audio AIR 192 6 ASIO');
haOut = audioDeviceWriter('SampleRate', Fs,...
    'Driver', 'ASIO', 'Device', 'M-Audio AIR 192 6 ASIO');

hfR1 = dsp.AudioFileReader(audioSource1, 'PlayCount', inf,...
    'SamplesPerFrame', Spf);
hfR2 = dsp.AudioFileReader(audioSource2, 'PlayCount', inf,...
    'SamplesPerFrame', Spf);
htS = dsp.TimeScope('SampleRate', Fs, 'TimeSpan', 0.1, 'BufferLength', 50e4,...
    'YLimits', [-1 1], 'ShowGrid', 1);
hsa = dsp.SpectrumAnalyzer('SampleRate', Fs, 'FrequencySpan', 'Start and stop frequencies',...
    'StartFrequency', 0, 'SpectralAverages', 10, 'StopFrequency', Fs/2,...
    'FrequencyResolutionMethod', 'WindowLength', 'PlotAsTwoSidedSpectrum', false,...
    'FrequencyScale', 'Log');

%% 外部VSTの定義
% Guitar Amp
Hamp = loadAudioPlugin(which('LeXtac.dll'));
setParameter(Hamp, 1, 1);     % Input
setParameter(Hamp, 2, 0);     % Ch
setParameter(Hamp, 3, 1);     % Drive
setParameter(Hamp, 4, 0.6);     % Low
setParameter(Hamp, 5, 0.7);     % Mid
setParameter(Hamp, 6, 0.5);     % High
setParameter(Hamp, 7, 0.7);     % Contour
setParameter(Hamp, 8, 1);     % PreEQ
setParameter(Hamp, 9, 1);     % Structure
setParameter(Hamp, 10, 1);     % Plexi
setParameter(Hamp, 11, 1);     % Boost
setParameter(Hamp, 12, 0.88);     % Out
setParameter(Hamp, 14, 1);     % PowerAmp
setParameter(Hamp, 15, 1)       % Mono
parameterTuner(Hamp)

% Hdly
%% MATLAB Effects
% Cabinet response is loaded from IR file
hFilter = dsp.FIRFilter('NumeratorSource', 'Input port');
num{1} = [1;0]; tapLengthMax = length(num{1});
for n = 1:numel(data)
    num{n+1} = data{n}.Amp;
    tapLengthMax = max([length(num{n}), length(num{n+1}), tapLengthMax]);
end
for n = 1:numel(num)
    num{n} = padarray(num{n}, [tapLengthMax-length(num{n}) 0], 'post');
end

%% Stopボタン押されるまでストリーム実行
while stpBtn   % Stream
    % Input Select
    if sourceSel == 0
        sourceSig = step(haIn);
    elseif sourceSel == 1
        sourceSig = step(hfR1);
    else
        sourceSig = step(hfR2);
    end
    temp = process(Hamp, [sourceSig sourceSig]);    % VST Guitar Amp
    % Cabinet Select 1-10. 1=Thru
    temp = step(hFilter, temp, num{cabSel+1}');
    step(haOut, temp);
    if viewTSel == 1
        step(htS, temp);
    end
    if viewFSel == 1
        step(hsa, temp);
    end
    drawnow;
end
stpBtn = true;

%% release
release(haIn)
release(haOut)
release(hfR1)
release(hfR2)

次回は作成したものからVSTプラグインを生成してみようと思う。
MATLABでIR(Impulse Response) Speaker Simulator VST Plugin 3につづく