2012/03/17からのアクセス回数 31448
教科書や雑誌にでているグラフを自分で再現してみることで、理解が深まり応用の幅も広まります。
ここでは、octaveを使って雑誌インターフェース2006/09号の8章の図13を再現して、IIRフィルタの伝達関数の設計方法を再現していきます。
8章IIRフィルタの伝達関数の設計方法にある図13*1を以下に引用します。
アナログ・フィルタからデジタル・フィルタへの変換には、以下のs-z変換を使って伝達関数をH(s)からH(z)に変換します。
$$ s = \frac{2}{T}\frac{1 - z^{-1}}{1 + z^{-1}} $$
しかし、この変換のみだと図13の(b)のようにカットオフ周波数が想定した30Hzではなく、25Hzあたりずれてしまいます。
これは、アナログ周波数\( \omega_A \) と デジタル周波数\(\omega_D \)が歪んでいるために生じます。*2
そこで、このひずみを考慮したカットオフ周波数を計算し(プリワーピングと呼ぶ)アナログ・フィルタを作成します。
プリワーピングの式は、以下の様になります。
$$ \omega_D = \frac{2}{T} tan(\omega_A T / 2) $$
octaveで提供されている周波数応答表示関数freqs, freqzでは振幅がdB単位で周波数もナイキスト周波数を1に正規化した表示でプロットされるので、これを振幅と周波数でプロットするように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
図13で使っている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);
図13の(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);
求める図13の(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} \) |
それではバタワースフィルタを使って、
を作っていきます。
求めるフィルタは、図19のようにします。*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; % 5次のアナログ基準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); % プロット
皆様のご意見、ご希望をお待ちしております。