手習ひデジタル信号処理(62) Scilab、伝達関数H(z)から位相プロット、Fs対応版

Joseph Halfmoon

正規化周波数表示はカッコいいけど実際の周波数でないと実感がわかないなどとブーたれて、前回は伝達関数からゲインプロットを実周波数で表示する関数を作ってみました。今回は、位相プロットで実周波数表示を行う関数を作ってみたいと思います。素人が良く理解しないまま試行錯誤しているのだけれども大丈夫か?自分。

※「手習ひデジタル信号処理」投稿順 Indexはこちら

※処理のターゲットは教科書の例題の1次の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)の「正規化周波数」でプロットされます。こんな感じ。

NormalizedBode

これでもいいしょ、と思わないこともないけれど、年寄は具体的な周波数が入ってほしい(サンプリング周波数を乗ずればよいのだけれど。)

位相プロット用の関数

前回の振幅プロット用と同様に、第1引数にサンプリング周波数、第2引数に伝達関数をとるプロット用の関数が以下に。今回、第3引数にオプション追加。オプションが0だと前回同様周波数軸はリニアに、1だと対数となります。

// 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

こんな感じ。phasePlotLog

つづいて、周波数軸をリニアにとってみました。

plotPH(10000, Hz, 0); // Phase Lin

今度はこんな感じ。phasePlotLin

振幅プロット用の関数にもオプション

前回の振幅プロット用関数は周波数軸リニアだけだったので、上記の位相用と同じく第3引数にオプション追加。オプションが0だと周波数軸はリニアに、1だと対数。

// 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

こんな感じ。MagPlotLog

つついて、周波数軸リニアの場合のゲインプロット。

plotMAG(10000, Hz, 0); // Mag Lin

これで前回と同じパターンを得ます。

MagPlotLin

グラフを描くと、あれもこれもと追加機能が欲しくなりますが、とりあえず当初の目論見の実周波数でプロットすることは出来たみたい。大丈夫か?

手習ひデジタル信号処理(61) Scilab、伝達関数H(z)からゲイン線図、Fs対応版 へ戻る

手習ひデジタル信号処理(63) Scilab、差分方程式を直接計算、1次IIRフィルタ へ進む