FrontPage

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

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

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

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

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

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

fig13.png

アナログ・フィルタからデジタル・フィルタへの変換には、以下の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

fig11.png

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

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

$$ \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

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

図13で使っている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);

fig1.png

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

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

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

fig2.png

図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);

fig3.png

求める図13の(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(バンドリジェクト・フィルタ)

を作っていきます。

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

fig19.png

ここで問題になるのは、バタワースフィルタの次数をどうするかです。 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;
% 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の周波数応答は以下の様になります。

fig4.png

バタワースの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の周波数応答は以下の様になります。

fig5.png

バタワースの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の周波数応答は以下の様になります。

fig6.png

バタワースの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の周波数応答は以下の様になります。

fig7.png

フィルタの縦結合

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

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);    % プロット

fig8.png

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

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を使った変換
% 次数を1とした場合
F = linspace(0, pi, 512);
[B, A] = invfreq(H, F, 1, 1);
MyFreqz(B, A, Ft);

次数を1とした場合

fig9.png

次数を20とした場合

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

fig10.png

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

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

コメント

選択肢 投票
おもしろかった 19  
そうでもない 1  
わかりずらい 0  

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


(Input image string)


*1 Interface2006/09 p127より引用
*2 Interface2006/09 p127より引用
*3 Interface2006/09 p130より引用

添付ファイル: filefig9.png 611件 [詳細] filefig10.png 625件 [詳細] filefig8.png 607件 [詳細] filefig7.png 600件 [詳細] filefig6.png 622件 [詳細] filefig5.png 649件 [詳細] filefig4.png 692件 [詳細] filefig19.png 668件 [詳細] filefig3.png 581件 [詳細] filefig2.png 673件 [詳細] filefig1.png 698件 [詳細] filefig11.png 740件 [詳細] filefig13.png 783件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2016-02-05 (金) 14:58:57 (411d)
SmartDoc