正規化周波数表示はカッコいいけど実際の周波数でないと実感がわかないなどとブーたれて、前回は伝達関数からゲインプロットを実周波数で表示する関数を作ってみました。今回は、位相プロットで実周波数表示を行う関数を作ってみたいと思います。素人が良く理解しないまま試行錯誤しているのだけれども大丈夫か？自分。
※「手習ひデジタル信号処理」投稿順 Indexはこちら
※処理のターゲットは教科書の例題の１次のIIRフィルタです。伝達関数および、ブロック線図などはお手数ですが前回をご覧ください。
Scilab上での伝達関数の定義
前回同様、syslin関数で 伝達関数 Hz を定義しているところが以下に。
a1=0.996;
b0=1-a1;
Hz = syslin('d', b0/(1.0-a1*%z));
上記の伝達関数 Hz をそのままボード線図を描く関数 bodeに与えればゲインプロットも位相プロットも一度に得ることができます。
bode(Hz);
実際のプロットが以下に。周波数軸はサンプリング周波数を「1Hz」としたとき（ナイキスト周波数0.5Hz)の「正規化周波数」でプロットされます。こんな感じ。
これでもいいしょ、と思わないこともないけれど、年寄は具体的な周波数が入ってほしい（サンプリング周波数を乗ずればよいのだけれど。）
位相プロット用の関数
前回の振幅プロット用と同様に、第１引数にサンプリング周波数、第２引数に伝達関数をとるプロット用の関数が以下に。今回、第３引数にオプション追加。オプションが０だと前回同様周波数軸はリニアに、１だと対数となります。
// Phase plot
// Fs: sampling frequency [Hz]
// Hz: Transfer Function
// Opt: 0=Linear 1=Log
function plotPH(Fs, Hz, opt)
[frq, rep]=repfreq(Hz);
clf();
gopt='nn'
if opt then
gopt='ln'
end
plot2d(gopt, frq * Fs, (180/%pi)*atan(imag(rep), real(rep)));
xtitle('Phase plot','frequency [Hz]','Phase [degree]');
endfunction
さて上記の関数を使って、前述の伝達関数Hzの位相プロットを描いてみます。まずは周波数軸対数の場合。
plotPH(10000, Hz, 1); // Phase Log
こんな感じ。
つづいて、周波数軸をリニアにとってみました。
plotPH(10000, Hz, 0); // Phase Lin
今度はこんな感じ。
振幅プロット用の関数にもオプション
前回の振幅プロット用関数は周波数軸リニアだけだったので、上記の位相用と同じく第３引数にオプション追加。オプションが０だと周波数軸はリニアに、１だと対数。
// Magnitude plot
// Fs: sampling frequency [Hz]
// Hz: Transfer Function
// Opt: 0=Linear 1=Log
function plotMAG(Fs, Hz, opt)
[xm, fr]=frmag(Hz, 4096);
clf();
gopt='nn'
if opt then
gopt='ln'
end
plot2d(gopt, fr * Fs, 20*log10(xm));
xtitle('Magnitude plot','frequency [Hz]','magnitude [dB]');
endfunction
まずは周波数軸対数の場合のゲインプロット。
plotMAG(10000, Hz, 1); // Mag Log
こんな感じ。
つついて、周波数軸リニアの場合のゲインプロット。
plotMAG(10000, Hz, 0); // Mag Lin
これで前回と同じパターンを得ます。
グラフを描くと、あれもこれもと追加機能が欲しくなりますが、とりあえず当初の目論見の実周波数でプロットすることは出来たみたい。大丈夫か？