MATLABでGUI付きアプリ開発 ～ App Designer を使おう（3）～タイマー割り込みでスペクトル表示（非同期FIFOを使う）
MATLABでGUI付きアプリ開発 ～ App Designer を使おう（1）～フィルター特性を変えてグラフ表示
MATLABでGUI付きアプリ開発 ～ App Designer を使おう（2）～ストリーミングオーディオをフィルター処理
MATLABでGUI付きアプリ開発 ～ App Designer を使おう（4）～ フィルター周波数・バンド幅をマウスで操作
まえがき
（１）、（２）に引き続き、今回はフィルター処理前後のパワースペクトル表示を追加してみます。
毎フレーム描画をすると動作が重くなるので、タイマー割り込みを使って0.5秒毎に表示をします。レイテンシーが大きくならないよう、フィルター処理とパワースペクトル計算を非同期FIFOを使って分けて処理をします。
要素技術
パワースペクトル推定
今回はパワースペクトル推定に、Signal Processing Toolbox の pwelch 関数を使います。
パワースペクトル推定は色々とややこしくて私も解説できないので詳しい説明は割愛させていただきますが、pwelch を使うと「ややこしい」ことを自動でやってくれます。
（pwelch自体もパラメータ指定パターンがありすぎてややこしいですが・・）
具体的には、オーバーラップフレーム分割～窓掛け～FFT～パワー算出～窓重み補正（PSDの場合は＋サイドローブ補正）～負の周波数成分折り返し～時間平均、辺りを自動でやってくれます。
デフォルトはパワースペクトル密度(PSD)で、パワースペクトルが欲しい場合は 'power' オプションを指定します。
オーバーラップ等を行わない、単純な例で見てみましょう。
入力信号は、100Hzの正弦波＋ノイズとしています。
rng default
Fs = 1000;
N = 512;
t = (0:N-1)/Fs;
freq = 0:Fs/N:Fs/2;
s = cos(2*pi*100*t);
x = s + randn(size(t));
figure
[psdx, f] = pwelch(x,rectwin(N),0,[],Fs,'power');
pdB1 = 10*log10(psdx);
plot(f,pdB1)
grid on
title('Power Spectrum Using pwelch')
xlabel('Frequency (Hz)')
ylabel('Power (dB)')
これと同じと思われることをFFTでやってみます。
figure
xdft = fft(x);
xdft = xdft(1:N/2+1);
psdx = abs(xdft/N).^2;
psdx(2:end-1) = 2*psdx(2:end-1);
pdB2 = 10*log10(psdx);
plot(freq,pdB2)
grid on
title('Power Spectrum Using FFT')
xlabel('Frequency (Hz)')
ylabel('Power (dB)')
差分を見てみましょう。
figure
plot(pdB1' - pdB2)
良さそうですね。
フーリエ変換の定義式から、f=0、fs/2以外では正負の周波数成分に分かれるので、（pwelchを使わずFFTでやる場合）パワースペクトル推定ではそれらのパワーを足し合わせる必要があります。入力が実数であれば正負の大きさは同じ（複素共役）なので、2倍すればOKです。
離散フーリエ変換結果（係数）はフレーム長を偶数Nとすると、0～N/2まで（N/2+1個） がベースバンド信号で、それより上は折り返し成分であり、負の周波数成分と等価です。
詳しくはこちらをご覧ください。
折り返し雑音とは？～なぜ車のホイールは逆回転するのか～デジタルはアナログの近似ではない
実際に係数を見てみましょう。
rng(5)
N = 16;
x = randn(N,1);
F = fft(x);
figure
y = abs(F);
stem(0:N-1,y)
grid
title('abs')
hold on
for k=1:N/2-1
plot([k, N-k],[y(k+1), y(k+1)],'r:')
end
0 を中心にシフトした方が分かりやすいかもしれません。
rng(5)
N = 16;
x = randn(N,1);
F = fftshift(fft(x));
figure
y = abs(F);
stem(-N/2:N/2-1,abs(y))
grid
title('abs')
hold on
for k=1:N/2-1
plot([k-N/2, N-k-N/2],[y(k+1), y(k+1)],'r:')
end
figure
stem(-N/2:N/2-1,real(F))
grid
title('実数部')
figure
stem(-N/2:N/2-1,imag(F))
grid
title('虚数部')
また離散信号解析を行う場合、周波数分解能（FFT bin）は Fs/FFTn となり周辺のパワーが足され、フレーム長を変えると定常信号（ノイズ等）のレベルは変わってしまいます。
そのため、定常信号の場合は周波数分解能で正規化する PSD が用いられることも多いと思います。
FFTで求める場合は、(Fs*FFTn)で割れば良いことになります。
figure
[psdx, f] = pwelch(x,rectwin(N),0,[],Fs);
pdB3 = 10*log10(psdx);
plot(f,pdB3)
grid on
title('PSD Using pwelch')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')
figure
xdft = fft(x);
xdft = xdft(1:N/2+1);
psdx = abs(xdft).^2 / (Fs*N);
psdx(2:end-1) = 2*psdx(2:end-1);
pdB4 = 10*log10(psdx);
plot(freq,pdB4)
grid on
title('PSD Using FFT')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')
figure
plot(pdB3' - pdB4)
こちらも良さそうです。
今回はノイズと信号を見たいわけではなく音楽信号を想定していますし、フィルタリング前後の相対比較ができれば良いのでどちらでも良いとは思います。
このようにMATLABのハイレイヤー関数を使う場合は、一度それを使わずに自分で組んで確認してみることをお勧めします。
正解がそこにあるので自分に勘違い等ないかすぐ確認することができ、実践的で効率的な勉強ができますし、組み込み系に移行する場合なども即対応することができます。
非同期バッファ
DSP System Toolbox には dsp.AsyncBuffer という非同期FIFOオブジェクトがあります。
FIFO とは "First In First Out" の略で、データを入れた順番に読み出すことができるようになっています。
スペクトル解析するにはある程度のサンプル数がないと周波数分解能が確保できません。
こちらは Fig. 2（N=512） を N=128 で計算し直した場合のグラフです。
周波数分解能の違いがお分かりになるかと思います。
一方長すぎると、今度は時間分解能が落ちるため過渡信号の解析ができず、処理も重くなりますし、レイテンシーも大きくなりGUIの反応が鈍くなってしまいます。
ここでは、フィルター処理・音声出力用は256サンプルずつ処理を行い、スペクトル解析用は4096サンプル溜まってから行っています。（pwelchの中で、50%オーバーラップの最大8ブロックに分割されてFFTされます）
この様な非同期処理の場合、普通はデータがどこまで溜まっているか・残っているか等をケアしなければいけませんが、AsyncBuffer を使うとデータの入出力だけすればよくなります。
実際の使い方を見ていただいた方が早いと思うので、以下のコードをご参照ください。（256, 4096は別途指定しています）
buf = dsp.AsyncBuffer;
while app.PLAY
u = hAudioSource(); % 256サンプルずつ読み込み
% フィルター・音声出力処理
write(buf,u); % 256サンプルをFIFOに追加
if buf.NumUnreadSamples >= pFrameSize % 4096サンプル溜まったら実行
in = read(buf,pFrameSize); % FIFOから、最初に入れた4096サンプルを読み出し
% スペクトル解析
end
endFIFOに書き込むサンプル数は、毎回違っても大丈夫です。
今回はpwelch自体でオーバーラップ処理ができるのでそちらしか使っていませんが、オーバーラップ読み出しもできます。
アプリ
では前回作ったアプリにさらに、今回の機能を加えていきましょう。
まずはキャンバスサイズを大きくします。縦を680にしてみました。
UIAxes は小さくします。この辺はお好きなように。
位置やサイズはマウスでも変えられますし、右の「位置」-> Position でも変えられます。左下を原点とする、[x, y, W, H] で指定します。
GUIコンポーネントの追加
・座標軸（UIAxes2）
Title String：Spectrum
XLavel String：Frequency (Hz)
YLabel String：Power (dB)
とします
今回はこれのみです。お好きなサイズに変更してください。
プロパティ
以下を追加します。
SPCTimer; % timer object
timerSet = false;
initFlag = true; % 1st draw of the Spectrum
pFrameSize = 4096; % welch frame size
buf1 = dsp.AsyncBuffer; % FIFO for BYPASS signal
buf2 = dsp.AsyncBuffer; % FIFO for filtered signal
pwelch_F = zeros(0,1), pwelch_P1 = ones(0,2), pwelch_P2 = ones(0,2); % for pwelch()
ax1, ax2; 内部関数
前回追加した drawFrequencyResponse() の下に、タイマー割り込みで実行される関数 SPCTimerFcn() を追加します。2回目以降は y 軸データの中身（YData）のみを更新することによって、描画の高速化を図っています。
% timer function
function SPCTimerFcn(app,~,~)
if app.initFlag
app.ax1 = semilogx(app.UIAxes2,app.pwelch_F,10*log10(app.pwelch_P1),':',LineWidth=1);
hold(app.UIAxes2,'on')
app.ax2 = semilogx(app.UIAxes2,app.pwelch_F,10*log10(app.pwelch_P2));
hold(app.UIAxes2,'off')
grid(app.UIAxes2,'on')
xlim(app.UIAxes2,[20 app.fs/2])
ylim(app.UIAxes2,[-120 0])
legend(app.UIAxes2,'L-BYPASS','R-BYPASS','L-Filtered','R-Filtered',Location='northeast')
else
% update data only
app.ax1(1).YData = 10*log10(app.pwelch_P1(:,1));
app.ax1(2).YData = 10*log10(app.pwelch_P1(:,2));
app.ax2(1).YData = 10*log10(app.pwelch_P2(:,1));
app.ax2(2).YData = 10*log10(app.pwelch_P2(:,2));
app.initFlag = false;
end
endstartupFcn
drawFrequencyResponse(app); の下に、スペクトル表示用変数の初期化と、タイマーオブジェクトの設定を追加します。
% intialize for 1st drawing of the Spectrum
[app.pwelch_P1,app.pwelch_F] = pwelch(zeros(app.pFrameSize,2));
[app.pwelch_P2,~] = pwelch(zeros(app.pFrameSize,2));
if ~app.timerSet
% Create timer object
app.SPCTimer = timer(...
'ExecutionMode', 'fixedRate', ... % Run timer repeatedly
'Period', 0.5, ... % Period in second
'BusyMode', 'queue',... % Queue timer callbacks when busy
'TimerFcn', @app.SPCTimerFcn); % Specify callback function
end
app.timerSet = true;UIFigureCloseRequest
delete(app) の前に以下を追加します。タイマーが走っていれば止めて、delete します。
if strcmp(app.SPCTimer.Running, 'on') % if timer is running then stop
stop(app.SPCTimer);
end
delete(app.SPCTimer); % delete timer object
t = timerfindall;
if ~isempty(t)
delete(t)
end最後の4行は通常必要ありませんが、デバッグ中で途中で止まってしまったときなど、タイマーが残って増えてしまうことがあります。
そのため、一応残ったタイマーオブジェクトがあれば消しています。
コールバック
・PLAYButton
追加はそれほど多くありませんが、数カ所に分かれているので全体を載せておきます。タイマーの ON/OFF と、256サンプルずつFIFOに書き込み（buf1：BYPASS信号用、buf2：フィルター処理済み信号用）、4096サンプル溜まったらパワースペクトル計算を行っています。
再生中にアプリを閉じてしまったときの処理を最後に入れてあります。
function PLAYButtonValueChanged(app, event)
value = app.PLAYButton.Value;
if ~app.setFile % file isn't set
app.PLAYButton.Value = false;
return;
end
app.PLAY = value;
if app.PLAY
app.PLAYButton.Text = 'STOP';
% If timer is not running, start it
if strcmp(app.SPCTimer.Running, 'off')
start(app.SPCTimer);
end
else
app.PLAYButton.Text = 'PLAY';
% Stop the timer
stop(app.SPCTimer);
end
if app.PLAY
hAudioSource = dsp.AudioFileReader(app.inputFilename,'SamplesPerFrame',app.FrameSize);
app.fs = hAudioSource.SampleRate;
hAudioOut = audioDeviceWriter('DeviceName', app.Output_Device, 'SampleRate', app.fs,...
'SupportVariableSizeInput', true, 'BufferSize', app.BufferSize);
while app.PLAY
u = hAudioSource();
if size(u,2) >= 2
u_st = u(:,1:2);
else % mono
u_st = [u(:,1) u(:,1)];
end
wetOut = app.octFilt(u_st); % filtering
write(app.buf1,u_st);
write(app.buf2,wetOut);
if app.buf1.NumUnreadSamples >= app.pFrameSize
in = read(app.buf1,app.pFrameSize);
[app.pwelch_P1,app.pwelch_F] = pwelch(in,[],[],[],app.fs,'power');
in = read(app.buf2,app.pFrameSize);
[app.pwelch_P2,~] = pwelch(in,[],[],[],app.fs,'power');
end
if isDone(hAudioSource) % Loop Play
reset(hAudioSource);
end
if app.BYPASS
out = u_st;
else
out = wetOut;
end
hAudioOut(out); % sound out
drawnow limitrate % for GUI interrupt
if ~isvalid(app); break; end % check unexpected force termination by user
end
reset(hAudioSource);
end
end実行結果
スクリプト
まとめ
これで以前の記事で詳しい説明をしなかったアプリの内容もほぼ説明できたかと思います。
一度慣れてしまうととても簡単にGUIアプリが作れますので、ぜひお試しを！
この記事が気に入ったらサポートをしてみませんか？