前回はFIRローパスフィルタの係数を求め、求めた係数を使って「単位」インパルス入力信号をフィルタして時間応答を求め、さらにそれをFFTして周波数特性のプロットをいたしました。まどろっこしいデス。フイルタ係数を求めるときに「ついでに」周波数特性のプロットも見て、それからOKしたいです。知らんけど。今回は微妙な改善?
※「手習ひデジタル信号処理」投稿順 Indexはこちら
参照させていただいているScilabの日本語版HELPページ
以下のページに「FIRおよびIIRフィルタの振幅」を求めるための関数 frmag についての説明があります。
この関数を使えば、前回の結果から周波数振幅特性が得られるんでないかい。
微妙な改善後の関数
Scilabのffilt関数に「1枚かぶせて」サンプリング周波数Fsと実際の(正規化でない)カットオフ周波数で指定できるようにしたフィルタ係数算出関数を前回作成してみました。今回、振幅周波数特性も同時に得られるように改善してみたものが以下に。
// Calculate and Plot FIR Lowpass Filter // Fs: sampling frequency [Hz] // nFlt: number of filter order // Fc: cutoff frequency [Hz] // Plt: 1 for plot, 0 for no plot function coeff=genFIRLP2(Fs, nFlt, Fc, Plt) Fcnorm = Fc / Fs; coeff = ffilt("lp", nFlt, Fcnorm); if Plt == 1 then [hzm,fr]=frmag(coeff,256) clf(); plot2d(fr*Fs, 20*log10(hzm)) xtitle('Magnitude plot','frequency [Hz]','magnitude [dB]'); end endfunction
上記の関数を実行してみた結果
前回同様のサンプリング周波数10kHz、FIRフィルタの次数100、カットオフ周波数500Hzで振幅周波数プロットあり(1)のときの関数呼び出し例が以下に。
Fs=10000; coeff=genFIRLP2(Fs, 100, 500, 1);
そしてそのとき出力されるプロットはこんな感じ。
前回のは縦軸が素のままだったですが、今回はdBに改まっております。いい感じでないかい。。。ホントか?