octave/IIRフィルタの伝達関数設計

FrontPage

テキストのグラフを再現してみよう

2012/03/17からのアクセス回数 30872

テキストのグラフを再現してみよう †

教科書や雑誌にでているグラフを自分で再現してみることで、理解が深まり応用の幅も広まります。

ここでは、octaveを使って雑誌インターフェース2006/09号の８章の図１３を再現して、IIRフィルタの伝達関数の設計方法を再現していきます。

↑

アナログ・フィルタからデジタル・フィルタへの変換 †

８章IIRフィルタの伝達関数の設計方法にある図１３*1を以下に引用します。

アナログ・フィルタからデジタル・フィルタへの変換には、以下のs-z変換を使って伝達関数をH(s)からH(z)に変換します。

$$ s = \frac{2}{T}\frac{1 - z^{-1}}{1 + z^{-1}} $$

しかし、この変換のみだと図１３の(b)のようにカットオフ周波数が想定した30Hzではなく、25Hzあたりずれてしまいます。

これは、アナログ周波数$ \omega_A $ とデジタル周波数$\omega_D $が歪んでいるために生じます。*2

そこで、このひずみを考慮したカットオフ周波数を計算し（プリワーピングと呼ぶ）アナログ・フィルタを作成します。

プリワーピングの式は、以下の様になります。

$$ \omega_D = \frac{2}{T} tan(\omega_A T / 2) $$

↑

周波数応答関数の修正 †

octaveで提供されている周波数応答表示関数freqs, freqzでは振幅がdB単位で周波数もナイキスト周波数を１に正規化した表示でプロットされるので、これを振幅と周波数でプロットするようにMyFreqs, MyFreqz関数を作成します。

function H = MyFreqs(b, a, Ft)
    Fn = Ft/2;
	w = linspace(0, 2*pi*Fn, 512);
	f = w/(2*pi);
	% アナログフィルターの周波数特性を取得
	H = freqs(b, a, w);
	% プロット
	figure; 
	subplot(2, 1, 1);
	plot(f, abs(H));
	grid on;
	title('Filter response')
	xlabel('Frequency(Hz)')
	ylabel('Magnitude')
	
	subplot(2, 1, 2);
	plot(f, angle(H)*180);
	grid on
	xlabel('Frequency(Hz)')
	ylabel('Phase(degree)')
endfunction

function [H, W] = MyFreqz(b, a, Ft)
	[H, W] = freqz(b, a, 512, "half", Ft);
	Fn = Ft/2;
	w = linspace(0, 2*pi*Fn, 512);
	f = w/(2*pi);
	figure; 
	subplot(2, 1, 1);
	plot(f, abs(H));
	grid on;
	title('Filter response')
	xlabel('Frequency(Hz)')
	ylabel('Magnitude')
	
	subplot(2, 1, 2);
	plot(f, angle(H)*180);
	grid on
	xlabel('Frequency(Hz)')
	ylabel('Phase(degree)')
endfunction

↑

1次ローパス・フィルタの仕様 †

図１３で使っている1次ローパス・フィルタは、以下の様な仕様になっています。

サンプリング周波数 $F_T$は、100 Hz
カットオフ周波数$ F_C $は、30 Hz

でアナログ・フィルタの伝達関数は以下の様になります。

$$ G(s) = \frac{1}{ 1 + \frac{s}{2 \pi F_C}} $$

これをoctaveで以下の様に計算します。

Ft = 100;   % サンプリング周波数(Hz)
T = 1/Ft;   % サンプリング周期(sec)
Fn = Ft/2;  % ナイキスト周波数(Hz)
Fc = 30;    % 遮断周波数(Hz)
Wa = 2*pi*Fc;
a = [1/Wa 1]; b = 1;

% アナログフィルターの周波数特性を表示
MyFreqs(b, a, Ft);

仕様通り30Hzで振幅が0.7に減衰しています。

次にこのアナログ・フィルタをs-z変換して周波数特性を表示します。

% そのまま双一次s-z変換をすると
[Zb, Za] = bilinear(b, a, 1/Ft);
% 変換結果の周波数特性を表示
MyFreqz(Zb, Za, Ft);

図１３の(b)のようにカットオフ周波数が25Hzにずれてしまいます。

そこで、プリワーピングを使って $\omega_D$を新しいカットオフ周波数として計算します。 $$ \omega_D = \frac{2}{T} tan(\omega_A T /2) $$

% プリワーピイングを行う
Wd = 2/T * tan(Wa*T/2);
a = [1/Wd 1]; b = 1;
[Zb, Za] = bilinear(b, a, 1/Ft);
% プリワーピイングを行った結果の周波数特性を表示
MyFreqz(Zb, Za, Ft);

求める図１３の(c)の周波数特性を持つ伝達関数ができました。

↑

本格的なフィルタの作成 †

アナログ・フィルタからデジタル・フィルタへの変換方法が分かったので、本格的なデジタル・フィルタを作成してみます。

以下の手順で進めます。

カットオフ周波数1の基準ローパス・フィルタを作成
プリワーピイングでデジタルのカットオフ周波数を計算
sftrans関数を使って基準ローパス・フィルタを変換

sftransの使い方が分からなかったのでソースをみたところ以下の様な変換を行っています。

要求されるフィルタ	変換式
LPF	$ p = \frac{s}{\omega_c} $
HPF	$ p = \frac{\omega_c}{s} $
BPF	$ p = \frac{ s^2 + \omega_0^2}{s \omega_b} $
BRF	$ p = \frac{s \omega_b}{ s^2 + \omega_0^2} $

$ \omega_c = 2 \pi f_c $
$ \omega_b = \omega_{c2} - \omega_{c1} $
$ \omega_0 = \sqrt { \omega_{c1} \omega_{c2} } $

↑

基準LPFの次数 †

それではバタワースフィルタを使って、

LPF（ローパス・フィルタ）
HPF（ハイパス・フィルタ）
BPF（バンドパス・フィルタ）
BRF（バンドリジェクト・フィルタ）

を作っていきます。

求めるフィルタは、図１９のようにします。*3

ここで問題になるのは、バタワースフィルタの次数をどうするかです。 octaveには、buttordという関数が用意されており、これを使ってフィルタを実装するのに必要な次元を求めることができます。

% 基準LPFの次数
Wp = 20/50;
Ws = 30/50;
Rp = 2;     % 2dB
Rs = 20;    % 20dB
[n, Wc] = buttord(Wp, Ws, Rp, Rs)

結果は、

n =  5
Wc =  0.41636

のようにでます。これで、次数を5とすることにします。

↑

バタワースのLPFの作成 †

プリワーピングで20Hzに対応するデジタルのカットオフ周波数を計算し、基準LPFに周波数変換を施します。

% バタワースのLPFの例
Fc = 20;    % 遮断周波数(Hz)
Wc = 2*pi*Fc;
% ５次のアナログ基準LPFを作成する
[z, p, g] = butter(n,1,'s');
% 周波数変換でカットオフWcをセット
Wcp = 2/T*tan(Wc*T/2);   % プリワーピイング
[Z, P, G] = sftrans(z,p,g,Wcp,false);
[b, a] = zp2tf(Z, P, G);
[Zb, Za] = bilinear(b, a, 1/Ft);
MyFreqz(Zb, Za, Ft);    % プロット

求めるLPFの周波数応答は以下の様になります。

↑

バタワースのHPFの作成 †

バタワースのHPFは、sftransのstopをtrueにして周波数変換します。

% ハイパスデジタルフィルタへの変換
% 基準LPFからWcとstop=trueを指定してハイパスフィルターに変換
[Z, P, G] = sftrans(z,p,g,Wcp,true);
[b, a] = zp2tf(Z, P, G);
[Zb, Za] = bilinear(b, a, 1/Ft);
MyFreqz(Zb, Za, Ft);    % プロット

HPFの周波数応答は以下の様になります。

↑

バタワースのBPFの作成 †

バタワースのBPFは、Wpfp, Wphpとstop=falseにして周波数変換します。

% バンドパスフィルタへの変換
Fpf = 20; Wpf = 2*pi*Fpf;
Fph = 35; Wph = 2*pi*Fph;
% プリワーピイング
Wpfp = 2/T*tan(Wpf*T/2);
Wphp = 2/T*tan(Wph*T/2);
% 基準LPFからバンドパスフィルタへの変換
[Z, P, G] = sftrans(z,p,g,[Wpfp Wphp],false);
[b, a] = zp2tf(Z, P, G);
[Zb, Za] = bilinear(b, a, 1/Ft);
MyFreqz(Zb, Za, Ft);    % プロット

BPFの周波数応答は以下の様になります。

↑

バタワースのBRFの作成 †

バタワースのBRFは、BPFの設定に対してstop=trueにして周波数変換します。

% バンドエリミネート・フィルタへの変換
Fsf = 20; Wsf = 2*pi*Fsf;
Fsh = 35; Wsh = 2*pi*Fsh;
% プリワーピイング
Wsfp = 2/T*tan(Wsf*T/2);
Wshp = 2/T*tan(Wsh*T/2);
% 基準LPFからバンドパスフィルタへの変換
[Z, P, G] = sftrans(z,p,g,[Wsfp Wshp],true);
[b, a] = zp2tf(Z, P, G);
[Zb, Za] = bilinear(b, a, 1/Ft);
MyFreqz(Zb, Za, Ft);    % プロット

BRFの周波数応答は以下の様になります。

↑

フィルタの縦結合 †

求めたデジタル・フィルタを結合して複雑なフィルタを作成することができます。

20Hzと30Hzを遮断するフィルタを作成してみます。

% フィルターの縦結合
Fs1 = 19.5; Ws1 = 2*pi*Fs1;
Fs2 = 20.5; Ws2 = 2*pi*Fs2;
Fs3 = 29.5; Ws3 = 2*pi*Fs3;
Fs4 = 30.5; Ws4 = 2*pi*Fs4;
% プリワーピイング
Ws1p = 2/T*tan(Ws1*T/2);
Ws2p = 2/T*tan(Ws2*T/2);
Ws3p = 2/T*tan(Ws3*T/2);
Ws4p = 2/T*tan(Ws4*T/2);
% 基準LPFからバンドパスフィルタへの変換
[Z, P, G] = sftrans(z,p,g,[Ws1p Ws2p],true);
sys1 = zp2sys(Z, P, G, T);
[Z, P, G] = sftrans(z,p,g,[Ws3p Ws4p],true);
sys2 = zp2sys(Z, P, G, T);

[b, a] = sys2tf(sysmult(sys1, sys2));
[Zb, Za] = bilinear(b, a, 1/Ft);
MyFreqz(Zb, Za, Ft);    % プロット

↑

任意の周波数特性を持つフィルタ †

octaveには、周波数特性を与えて、その伝達関数を作成するための関数invfreqzがあります。

これを使って最初の1次LPFを再現してみます。

% 目的の周波数特性を作成
Wa = 2*pi*Fd;
a = [1/Wa 1]; b = 1;
w = linspace(0, 2*pi*Fn, 512);
H = freqs(b, a, w);
% invfreqzを使った変換
% 次数を１とした場合
F = linspace(0, pi, 512);
[B, A] = invfreq(H, F, 1, 1);
MyFreqz(B, A, Ft);

次数を1とした場合

次数を20とした場合

% 次数を20とした場合
[B, A] = invfreq(H, F, 20, 20);
MyFreqz(B, A, Ft);

如何でしょう、最初のアナログフィルタのカーブがよく表現できています。

このようにoctaveを使うとIIRのフィルタがとても簡単に作成でき、図化して確認することができます。是非、自分でoctaveを動かしてみてください。

↑

コメント †

皆様のご意見、ご希望をお待ちしております。

バンドパスフィルタを試そうとしたらこのようなエラーが出ました。どのようにすればよいですか？error: zpk: second argument 'p' must be a vector or a cell of vectors -- 名無し? 2018-11-24 (土) 12:16:30