MATLAB/Simulink で気軽に始める音響信号処理 ～ App Designer で GUIアプリも簡単開発！
MATLAB EXPO楽しかったですね～。（￣ー￣）
以前、ポスター発表して来場者人気投票で2位を頂いたこともあるのですが、今はちょっと他のことする余裕がないので、いずれまた・・。(-_☆)
今回、特にライトニングトークに触発され、MATLABでの音響信号処理を布教すべく、”初心者向け” MATLABを使ったオーディオ信号処理について書いてみました。信号処理自体ではなく、「信号処理アルゴリズム開発環境」的な説明です。結果的に、App Designer の説明がメインになりましたが・・。
以前の記事、
MATLABとSimulinkと私 (思考を止めない開発ツール)
夏休みはVST！ MATLABで楽々（？）VSTプラグイン開発
辺りも参考になるかと思います。
必要なツールボックス
MATLAB本体だけでも基本的なことはできますが、Audio Toolboxと、個人的にはSimulink（機能ブロック置いて繋ぐだけ！）もあった方が良いと思います。思い付いたことをすぐ試すにはSimulinkの方が便利です。
Audio Toolboxの動作には、Signal Processing Toolbox、DSP System Toolboxも必要です。
一般的な機械学習、ディープラーニング用にはそれぞれ、Statistics and Machine Learning Toolbox、Deep Learning Toolbox があります。
その他、単独アプリ化（実行にMATLABライセンス不要）にはMATLAB Compilerが、C出力やJUCEプロジェクト（VST用）出力にはMATLAB Coderが必要です。
残念ながら MATLAB Compilerと MATLAB Coder はHOME版にはありません。
製品詳細についてはMathWorksのサイトをご参照ください。
MATLAB
Simulink
Audio Toolbox
Signal Processing Toolbox
DSP System Toolbox
Statistics and Machine Learning Toolbox
Deep Learning Toolbox
Mac版SimulinkでiOSサポートパッケージ（ノンサポート）を入れれば、機能は限定的ですが iOSアプリ化もできます。
Simulink Support Package for Apple iOS Devices
ラズパイ向けサポートパッケージもあります。
Simulink Support Package for Raspberry Pi Hardware
メディアファイル読み書き
MATLABで、「ステレオオーディオファイル読み込み -> モノラル（ここでは2chでL=R）に変換 -> 再生とファイル保存」は、
[y,fs] = audioread(readFileName);
m = sum(y,2)/2;
sound(m,fs)
audiowrite([m m], writeFileName, fs); % L/R 両chに書き込みこれだけです！対応フォーマットはOSがサポートしている形式、Windowsであれば、mp3、m4a、aac、wav、mp4、flac、ogg、webm等を読み込めます。ビデオファイルからもそのままオーディオ部分を読み込め、m4a、wav、mp4、flac、ogg等で書き込みもできます。
Simulinkだとこんな感じ。
"Add" ブロックは、"加算範囲"：指定した次元、"次元"：2にします。
最初にモデルを作る時は、新規 -> Simulinkモデル で "Basic Audio Player" を選ぶと楽です。
これだと出力が1chになりますが、「2chでL=R」にしたい場合は "Matrix
Concatenate" ブロックで次元を増やします。ついでにL/Rも分離してからやってみましょう。
「ディレクトリ内のflacファイルを読み込みwavに変換して保存」、するにはMATLABでは以下のようにします。
dirlist = dir('*.flac'); % ファイル名一覧取得
for n=1:length(dirlist)
readFile = dirlist(n).name; % 各ファイル名取得
info = audioinfo(readFile); % 各オーディオファイルの属性取得
bps = info.BitsPerSample;
disp(readFile)
[y,fs] = audioread(readfile);
[~,writeFile,~] = fileparts(dirlist(n).name); % 拡張子削除
audiowrite([writeFile '.wav'], y, fs, 'BitsPerSample', bps);
disp('done')
clear y fs
end他の開発環境ではwavファイルを読むだけで大変だったり、読み込んで規定の処理をして再生するだけなら簡単なのに、なにかオリジナルの処理をして再生しようとすると途端に難しくなったりしますが、MATLABならそんなことはありません。
ファイルからでなく、マイク入力等外部入力も同等に扱えます。
Audio Toolboxがあれば、ASIOやMIDIも扱えます。（VST以外）
上にも書きましたが、他の開発環境では、「ステレオオーディオを入力 -> （オリジナルの信号処理を入れる場所） -> ステレオオーディオ出力」というフレームワークを作る（あるいはそのやり方を理解する）だけで結構苦労します。そこが瞬間的に作れるのも MATLAB の利点です。
いかがでしょう？
これだけでもMATLABを使ってみたくなりませんか？（￣ー￣）
可視化
plot、bar、semilogx、semilogy、loglog、spectrogram（Signal Processing Toolbox）、waterfall、・・等々、各種グラフ表示が一発です。一部の拡大やデータ数値表示、waterfall等の3Dグラフはマウスでグリグリ回転もできます。
この辺はとても説明しきれないので "MATLAB グラフ"とかで検索してみてください。
SimulinkでScope（波形の時間軸表示）を使う場合は、後で説明しますがサンプル／フレームベースを設定する必要がありますのでご注意を。
私はほとんど使ったことがないのですが、「信号アナライザー」というアプリも用意されています。ワークスペース上のデータをドロップするだけで、各種グラフ表示～スクリプトの生成ができます。
ただし、信号アナライザーは1曲丸ごととかやるとかなり時間が掛かるようなので、あらかじめ一部を切り取った方が良さそうです。
今回は割愛しますが、「ライブエディター」もメッチャ便利です！
フレーム処理
PC等でリアルタイム処理を行うには、基本的にはフレーム単位での処理を行います。
サンプル単位の処理だと、次のサンプルが入ってくるまでに1サンプル出力するための全ての処理を毎サンプル必ず行わないと間に合いません。
しかしフレーム処理であれば、例えば各chごと1024サンプルずつ読み込んで処理をすれば、1024/fsの時間、CD音源であれば1024/44100=約23msの間に処理が終わりさえすればよくなります。PCの割り込み等が入ってたまに処理が間に合わないことがあっても、平均して1サンプル分の処理が1サンプル分以内に終われば、リアルタイムに処理を行って音切れなく音が出せるということです。
サンプルベース処理
毎サンプル必ず時間T内に処理を終える必要がある
フレームベース処理
時間 N×T以内にN個のデータ処理が終われば良い
（フレーム分～のレイテンシーが発生）
もちろんそのフレーム分のレイテンシー（遅れ）が発生するので、長くすれば良いというものでもありませんが。
Simulinkでは各ブロックの処理がサンプル単位かフレーム単位選べますので、例えば単にディレイさせてフィードバックとかなら、ディレイブロックをフレーム単位にすることにより実行できます。
どうしてもサンプル単位の処理が必要な場合は、"From multimedia File" ブロックの「オーディオ チャネルあたりのサンプル数」を1に設定します。
この場合、ちょっと重い処理をするとすぐ音切れしたりしますので、インプリには注意が必要ですが、できないことはありません。（結構やってます）
MATLABスクリプトの場合も、DSP System Toolbox の System object群を使ってできる限りフレーム単位の処理とすることで高速化ができます。サンプル単位で行いたい部分はフレーム長のforループ中で処理をします。
MEX化
タップ数の多いFIR等はどうしても処理時間が掛かるので、ネイティブコードにして高速化して組み込むという手もあります。
手動でも頑張ればできますが、"MATLAB Coder" があると、mファイルからとても簡単に作れます。
例えば userFIR.m というスクリプトから、
codegen userFIR -args {coder.typeof(0, [16384 2], [1 0]), ～とかやると、double型で、次元1は最大16384の可変サイズ、次元2はサイズ2固定のパラメータを引数とするネイティブコードが一発で生成され、あとは普通の関数としてMATLABから呼び出せます。この辺の設定は対話形式でも行えます。
HOME版は "MATLAB Coder" 自体がありませんが、Simulinkの高速実行モードもないし、この機能だけでも付けてくれると嬉しいのですが・・。（；￢＿￢）
まあ、手軽な高速化手段はVST化ですかね。こちらはAudio Toolboxが動く環境であれば良いので、HOME版でもできます。
App Designer
簡単にGUIアプリが作れます。
英語ですがここに簡単な紹介ビデオがありますので、目を通されると良いかと思います。
Getting Started with App Designer
MATLAB Compiler があればスタンドアローンアプリ化できますが、HOME版にはMATLAB Compiler製品自体がないので、HOME版では毎回 appdesigner から、もしくはMATLABコマンドウィンドウからmlappファイルを起動（拡張子なしのファイル名で起動）する形になります。
> appdesignerで起動し、目的に近いテンプレート、もしくは空のテンプレートから作ることができます。
ここで短い対話式チュートリアルを見ることもできます。
メニューの 新規 -> アプリ でも起動できます。
作成済みを編集したい場合は、
> appdesigner ファイル名でOKです。
アプリ例
入出力デバイス設定、タイマー設定、グラフ表示、オーディオファイル読み込み・再生、と最低限必要と思われる要素を盛り込んだ簡単な例を置いておきます。（R2021aで作成、要Audio Toolbox）
（"RST Dev" をクリックする度にタイマーが増えていたのを修正 2021/10/18）
> appdesigner LR_SPLで起動します。
基本的にはpropertiesとstartupFcnに初期設定等を書き、あとはGUIパーツを置いてそのコールバックを書くだけです。
コールバックは最初はないので、
各コンポーネントを右クリック -> コールバック
で必要なコールバックルーチンを追加してからコードを書き込みます。
作成後は、
各コンポーネントを右クリック -> コールバック -> ～コールバックに移動
でそのコード部分に飛べます。
コールバック以外の関数を作りたい場合は、
左のコードブラウザー -> 関数 -> ➕
で追加します。
コードビューで勝手に書いてしまうと、コードブラウザーと整合が取れなくなることがあるので注意してください。
プロパティも同様に
コードブラウザー -> プロパティ -> ➕
で追加します。
コード説明
－プロパティ、割り込みルーチン
変数宣言とタイマー割り込みで実行する SPLTimerFcn() の設定を行っています。
properties (Access = private)
SPLTimer % Timer object
hAudioSource
Leq_L
Leq_R
Fs
FFTn
FrameSize
TimeInterval
TimeWeighting
freqWeight
bandWidth
timerSet = false;
end
methods (Access = private)
function SPLTimerFcn(app,~,~)
global pwelch_F pwelch_P Lspl_L Lspl_R xt XTickLabel graphSel
switch graphSel
case 1
% welch
% Graph1
semilogx(app.UIAxes,pwelch_F,10*log10(pwelch_P))
[yupper,~] = envelope(pwelch_P,30,"peak");
yupper(yupper<=0) = NaN;
hold(app.UIAxes,'on')
semilogx(app.UIAxes,pwelch_F,10*log10(yupper))
hold(app.UIAxes,'off')
grid(app.UIAxes,'on')
xlim(app.UIAxes,[20 app.Fs/2])
legend(app.UIAxes,'L','R','L Peak','R Peak',"Location","best")
case 2
% SPL meter
bar(app.UIAxes,mean(Lspl_L))
hold(app.UIAxes,'on')
bar(app.UIAxes,mean(Lspl_R),0.5)
hold(app.UIAxes,'off')
xticks(app.UIAxes,xt)
xticklabels(app.UIAxes,XTickLabel)
xlabel(app.UIAxes,'Center frequency (Hz)')
ylabel(app.UIAxes,'Gain (dB)')
grid(app.UIAxes,'on')
end
end
end－startupFcn()
起動時に実行されます。
入力デバイス設定、splMeterオブジェクト、timerオブジェクト生成を行っています。
function startupFcn(app)
% set input devices
global Input_Device Ex_In_Fs sa_in_Name_s2 fileInput graphSel meter Output_Device inputGain
inputGain = 0; % 0dB
% set input devices
info = audiodevinfo;
dev_in_n = length(info.input);
sa_in_Name_s2 = cell(dev_in_n,1);
sa_in = struct2cell(info.input);
sa_in_Name = sa_in(1,:,:);
sa_in_Name_s = squeeze(sa_in_Name);
for i=1:dev_in_n
idx = strfind(sa_in_Name_s{i},'(');
if isempty(idx)
idx = length(sa_in_Name_s{i});
else
idx = idx(end) - 1;
end
sa_in_Name_s2{i+1} = strtrim(sa_in_Name_s{i}(1:idx));
end
sa_in_Name_s2{1} = 'File';
app.InputDevDropDown.Items = sa_in_Name_s2;
fileInput = 2; % External
Input_Device = sa_in_Name_s2{fileInput};
app.InputDevDropDown.Value = Input_Device;
Ex_In_Fs = 48000;
app.Fs = Ex_In_Fs;
% set output devices
dev_n = length(info.output);
sa = struct2cell(info.output);
sa_Name = sa(1,:,:);
sa_Name_s = squeeze(sa_Name);
for i=1:dev_n
idx = strfind(sa_Name_s{i},'(');
if isempty(idx)
idx = length(sa_Name_s{i});
else
idx = idx(end) - 1;
end
sa_Name_s{i} = strtrim(sa_Name_s{i}(1:idx));
end
app.OutputDevDropDown.Items = sa_Name_s;
Output_Device = sa_Name_s{1};
app.OutputDevDropDown.Value = Output_Device;
app.LoadSongButton.BackgroundColor = [0.75 0.75 0.75];
app.LoadSongButton.FontWeight = 'normal';
if ~app.timerSet
app.FFTn = 8192;
app.FrameSize = floor(Ex_In_Fs/app.FFTn)*app.FFTn*1;
graphSel = 1; % FFT graph
app.bandWidth = '1/3 octave';
app.TimeWeighting = 'Fast';
app.TimeInterval = 0.1;
app.freqWeight = 'Z-weighting';
meter = splMeter('SampleRate',app.Fs,...
'Bandwidth',app.bandWidth,...
'TimeWeighting',app.TimeWeighting,...
'TimeInterval',app.TimeInterval,...
'FrequencyWeighting',app.freqWeight);
% Create timer object
app.SPLTimer = timer(...
'ExecutionMode', 'fixedRate', ... % Run timer repeatedly
'Period', 1, ... % Period in second
'BusyMode', 'queue',... % Queue timer callbacks when busy
'TimerFcn', @app.SPLTimerFcn); % Specify callback function
end
app.timerSet = true;
end－SPLMeterUIFigureCloseRequest()
アプリが閉じられたときの処理を書きます。
タイマーが走っていれば停止してから削除、念のため他に走ってるタイマーも全て削除しています。これは必要に応じて。
function SPLMeterUIFigureCloseRequest(app, event)
% Stop timer, then delete timer and app
if strcmp(app.SPLTimer.Running, 'on')
stop(app.SPLTimer);
end
delete(app.SPLTimer);
delete(timerfindall)
delete(app)
end－InputDevDropDownValueChanged()
InputDevのコールバックです。
'File' が指定されたかどうかで 'Load Song'ボタンのアトリビュートを変更しています。
function InputDevDropDownValueChanged(app, event)
global fileInput Input_Device sa_in_Name_s2
Input_Device = app.InputDevDropDown.Value;
f = strcmp(sa_in_Name_s2,Input_Device);
fileInput = find(f);
if fileInput ~= 1 % from device
app.LoadSongButton.BackgroundColor = [0.75 0.75 0.75];
app.LoadSongButton.FontWeight = 'normal';
else
app.LoadSongButton.BackgroundColor = [1 1 1]; % from file
app.LoadSongButton.FontWeight = 'bold';
end
end－STARTButtonValueChang()
STARTのコールバックです。
割り込みタイマーをスタートさせ、mファイル内の calcSPL() 関数を呼び出しています。
function STARTButtonValueChanged(app, event)
global START fileInput
START = app.STARTButton.Value;
if START
% If timer is not running, start it
if strcmp(app.SPLTimer.Running, 'off')
start(app.SPLTimer);
end
cmdstr = sprintf('calcSPL(%d,%d,%d,%d)',...
fileInput,app.FFTn,app.FrameSize,app.Fs);
run(cmdstr)
else
% Stop the timer
stop(app.SPLTimer);
reset(app.hAudioSource);
end
end－StoreButtonPushed()
Storeのコールバックです。
その時表示しているグラフを下に描き写しています。
オブジェクトのコピーができそうな気もするのですがうまくいかなかったので・・。うまくいかなかったついでに重ね描きにしています。
function StoreButtonPushed(app, event)
global pwelch_F pwelch_P Lspl_L Lspl_R graphSel xt XTickLabel
switch graphSel
case 1
% welch
% Graph1
semilogx(app.UIAxes_2,pwelch_F,10*log10(pwelch_P))
[yupper,~] = envelope(pwelch_P,30,"peak");
yupper(yupper<=0) = NaN;
hold(app.UIAxes_2,'on')
semilogx(app.UIAxes_2,pwelch_F,10*log10(yupper))
% hold(app.UIAxes_2,'off')
grid(app.UIAxes_2,'on')
xlim(app.UIAxes_2,[20 app.Fs/2])
legend(app.UIAxes_2,'L','R','L Peak','R Peak',"Location","best")
case 2
% SPL meter
bar(app.UIAxes_2,mean(Lspl_L))
hold(app.UIAxes_2,'on')
bar(app.UIAxes_2,mean(Lspl_R),0.5)
% hold(app.UIAxes_2,'off')
xticks(app.UIAxes_2,xt)
xticklabels(app.UIAxes_2,XTickLabel)
xlabel(app.UIAxes_2,'Center frequency (Hz)')
ylabel(app.UIAxes_2,'Gain (dB)')
grid(app.UIAxes_2,'on')
end
end－ClearButtonPushed()
Clearのコールバックです。
重ね描きをoffにしてクリアしています。
function ClearButtonPushed(app, event)
hold(app.UIAxes_2,'off')
cla(app.UIAxes_2,'reset')
end－GraphDropDownValueChanged()
Graphのコールバックです。
f特／SPL表示を切り替えます。
function GraphDropDownValueChanged(app, event)
global meter graphSel
value = app.GraphDropDown.Value;
f = strcmp(app.GraphDropDown.Items,value);
graphSel = find(f);
cla(app.UIAxes,'reset')
if graphSel == 2
release(meter);
meter = splMeter('SampleRate',app.Fs,...
'Bandwidth',app.bandWidth,...
'TimeWeighting',app.TimeWeighting,...
'TimeInterval',app.TimeInterval,...
'FrequencyWeighting',app.freqWeight);
end
end－WeightingDropDownValueChanged()
Weightingのコールバックです。
聴覚重み付け、A/C/Z（重み付けなし）を切り替えます。
function WeightingDropDownValueChanged(app, event)
value = app.WeightingDropDown.Value;
app.freqWeight = value;
global meter
release(meter);
meter = splMeter('SampleRate',app.Fs,...
'Bandwidth',app.bandWidth,...
'TimeWeighting',app.TimeWeighting,...
'TimeInterval',app.TimeInterval,...
'FrequencyWeighting',app.freqWeight);
end－OctoveDropDownValueChanged()
Octoveのコールバックです。
splMeterの周波数分解能を切り替えます。
function OctoveDropDownValueChanged(app, event)
value = app.OctoveDropDown.Value;
app.bandWidth = value;
global meter
release(meter);
meter = splMeter('SampleRate',app.Fs,...
'Bandwidth',app.bandWidth,...
'TimeWeighting',app.TimeWeighting,...
'TimeInterval',app.TimeInterval,...
'FrequencyWeighting',app.freqWeight);
end－LoadSongButtonPushed()
LoadSongのコールバックです。
音声ファイルのファイル名をセットして、情報を取得します。
function LoadSongButtonPushed(app, event)
global inputFilename
[FileName, PathName, ~] = uigetfile({'*.mp3;*.m4a;*.aac;*.wav;*.mp4;*.flac;*.ogg;*.webm','Music Files';'*.*','All Files'},'Select Song', 'MultiSelect', 'off');
if FileName ~= 0
inputFilename = [PathName FileName];
m_info = audioinfo(inputFilename);
m_length = floor(m_info.Duration);
app.Fs = m_info.SampleRate;
end
end－RSTDevButtonPushed()
RSTDevのコールバックです。
オーディオI/Fの設定を再読込します。
function RSTDevButtonPushed(app, event)
audiodevreset
startupFcn(app);
end－OutputDevDropDownValueChanged()
OutputDevのコールバックです。
出力デバイスを切り替えます。
function OutputDevDropDownValueChanged(app, event)
global START Output_Device
if ~START % while STOP state only
Output_Device = app.OutputDevDropDown.Value;
else
app.OutputDevDropDown.Value = Output_Device;
end
end－InputGainSliderValueChanged()
InputGainのコールバックです。
function InputGainSliderValueChanged(app, event)
global inputGain
inputGain = app.InputGainSlider.Value;
end－calcSPL()
STARTButtonValueChanged()で呼び出しているmファイルの中身です。
function calcSPL(fileInput,FFTn,FrameSize,Fs)
global START pwelch_P pwelch_F Input_Device inputFilename meter Lspl_L Lspl_R xt XTickLabel graphSel Output_Device inputGain
if (fileInput == 1) % read from file
hAudioSource = dsp.AudioFileReader(inputFilename,'SamplesPerFrame',FrameSize);
else
hAudioSource = audioDeviceReader('Device',Input_Device,'SamplesPerFrame',FrameSize,'SampleRate',Fs); % read from external input
end
hAudioOut = audioDeviceWriter('DeviceName', Output_Device, 'SampleRate', hAudioSource.SampleRate,...
'SupportVariableSizeInput', true, 'BufferSize', 4096);
overlapN = FFTn/4;
while START
u = step(hAudioSource) * db2mag(inputGain);
if size(u,2) >= 2
u_st = u(:,1:2);
else % mono
u_st = [u(:,1) u(:,1)];
end
switch graphSel
case 1
% welch
[pwelch_P,pwelch_F] = pwelch(u_st,ones(FFTn,1),overlapN,FFTn,Fs,'power');
if min( pwelch_P(:,2) ) <= 0
pwelch_P(pwelch_P<=0,2)
end
case 2
% SPL meter
[~,~,Lspl_L,~] = meter(u_st(:,1)); % take Lch
[~,~,Lspl_R,~] = meter(u_st(:,2)); % take Rch
cf = round(meter.getCenterFrequencies,2,'significant'); % get each center frequency of the octave filter
xt = 1:4:size(Lspl_L,2); % reduce the number of ticks
XTickLabel = compose('%d',cf(xt));
end
if (fileInput == 1) % read from file
if isDone(hAudioSource) % Loop Play
reset(hAudioSource);
end
end
hAudioOut(u); % soud out
drawnow limitrate
end
reset(hAudioSource);
end hAudioFile = dsp.AudioFileWriter(outputFilename, ～
と宣言しておいて、
hAudioOut(u);
のところで
hAudioFile(u)
とすれば、同時にファイルに書き込むこともできます。
Graph/Weighting/Octove 各 DropDown メニューは、以下のように設計ビューの方で設定しています。InputDev/Output のように（startupFcn() で設定）
app.InputDevDropDown.Items = （セル形式）;
としてコード上で入れることもできます。
使い方
"InputDev" で入力デバイス選択、あるいは "File"を選んでから"Load Song" でファイルを選び、"START" で解析を実行します。
出力は"OutputDev"で指定します。
再生時は変更できません。
"RST Dev"は、（MATLAB）立ち上げ後にI/F環境を変えた場合に再設定を行います。
"START"の中で、run(cmdstr)として実際の解析用 mファイルを引数付きで呼び出しています。
"Store" で下のグラフに保存します。
この説明のためにざっと作って細かいデバッグはしていないので、バグがあってもご容赦を！
エラーチェックも入れてないので、入力デバイスあるいはファイルをちゃんと設定してからSTARTさせてください。（；￢＿￢）
メインの処理もmlappファイル内に書いても良いのですが、外部関数を呼び出す例として別のmファイルとしています。
全部引数にするのが面倒なので、ほとんどを global 変数にしています。
構造体で持って行けば良いのかもしれませんが・・。
レイテンシーを小さくしたい場合は "FrameSize"、"BufferSize" を小さくする必要がありますが、その値は環境依存で小さくし過ぎると音切れするので、私は普段はこれらもメニューで選べるようにしています。
「ここの書き方はおかしい」「もっとこうした方が良い」等あればぜひ教えてください！
短めのを書くつもりが思いのほか長くなってしまいましたが、いかがでしたでしょうか？(^-^;
アルゴリズム開発からインプリまで、思考を止めることなく一貫して（比較的）簡単にできるのが MATLAB の魅力です。制御や解析系と比べ、オーディオ系で使っているユーザー数は圧倒的に少ない気がするのですが、見て聴いてすぐ確かめられるので、オーディオ系開発にはとても適していると思います。使ったことがない人も、（使える環境であれば）ぜひ使ってみてください！