![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
2012/03/17からのアクセス回数 30284
教科書や雑誌にでているグラフを自分で再現してみることで、理解が深まり応用の幅も広まります。
ここでは、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次ローパス・フィルタは、以下の様な仕様になっています。
でアナログ・フィルタの伝達関数は以下の様になります。
$$ 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)の周波数特性を持つ伝達関数ができました。
アナログ・フィルタからデジタル・フィルタへの変換方法が分かったので、本格的なデジタル・フィルタを作成してみます。
以下の手順で進めます。
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} \) |
それではバタワースフィルタを使って、
を作っていきます。
求めるフィルタは、図１９のようにします。*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とすることにします。
プリワーピングで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は、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は、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は、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を動かしてみてください。
皆様のご意見、ご希望をお待ちしております。
