octave/IIRフィルタの伝達関数設計
をテンプレートにして作成
[
トップ
] [
新規
|
一覧
|
単語検索
|
最終更新
|
ヘルプ
]
開始行:
[[FrontPage]]
#contents
2012/03/17からのアクセス回数 &counter;
* テキストのグラフを再現してみよう [#ec5c600a]
教科書や雑誌にでているグラフを自分で再現してみることで、...
ここでは、octaveを使って雑誌インターフェース2006/09号の８...
** アナログ・フィルタからデジタル・フィルタへの変換 [#k91...
８章IIRフィルタの伝達関数の設計方法にある図１３((Interfac...
&ref(fig13.png);
アナログ・フィルタからデジタル・フィルタへの変換には、以...
$$
s = \frac{2}{T}\frac{1 - z^{-1}}{1 + z^{-1}}
$$
しかし、この変換のみだと図１３の(b)のようにカットオフ周波...
これは、アナログ周波数\( \omega_A \) と デジタル周波数\(\...
&ref(fig11.png);
そこで、このひずみを考慮したカットオフ周波数を計算し（プ...
プリワーピングの式は、以下の様になります。
$$
\omega_D = \frac{2}{T} tan(\omega_A T / 2)
$$
** 周波数応答関数の修正 [#j59ffee6]
octaveで提供されている周波数応答表示関数freqs, freqzでは...
#pre{{
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次ローパス・フィルタの仕様 [#wc3cea2a]
図１３で使っている1次ローパス・フィルタは、以下の様な仕様...
- サンプリング周波数 \(F_T\)は、100 Hz
- カットオフ周波数\( F_C \)は、30 Hz
でアナログ・フィルタの伝達関数は以下の様になります。
$$
G(s) = \frac{1}{ 1 + \frac{s}{2 \pi F_C}}
$$
これをoctaveで以下の様に計算します。
#pre{{
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);
}}
&ref(fig1.png);
仕様通り30Hzで振幅が0.7に減衰しています。
次にこのアナログ・フィルタをs-z変換して周波数特性を表示し...
#pre{{
% そのまま双一次s-z変換をすると
[Zb, Za] = bilinear(b, a, 1/Ft);
% 変換結果の周波数特性を表示
MyFreqz(Zb, Za, Ft);
}}
&ref(fig2.png);
図１３の(b)のようにカットオフ周波数が25Hzにずれてしまいま...
そこで、プリワーピングを使って \(\omega_D\)を新しいカット...
$$
\omega_D = \frac{2}{T} tan(\omega_A T /2)
$$
#pre{{
% プリワーピイングを行う
Wd = 2/T * tan(Wa*T/2);
a = [1/Wd 1]; b = 1;
[Zb, Za] = bilinear(b, a, 1/Ft);
% プリワーピイングを行った結果の周波数特性を表示
MyFreqz(Zb, Za, Ft);
}}
&ref(fig3.png);
求める図１３の(c)の周波数特性を持つ伝達関数ができました。
** 本格的なフィルタの作成 [#l50fb268]
アナログ・フィルタからデジタル・フィルタへの変換方法が分...
以下の手順で進めます。
- カットオフ周波数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の次数 [#a7451550]
それではバタワースフィルタを使って、
- LPF（ローパス・フィルタ）
- HPF（ハイパス・フィルタ）
- BPF（バンドパス・フィルタ）
- BRF（バンドリジェクト・フィルタ）
を作っていきます。
求めるフィルタは、図１９のようにします。((Interface2006/0...
&ref(fig19.png);
ここで問題になるのは、バタワースフィルタの次数をどうする...
octaveには、buttordという関数が用意されており、これを使っ...
#pre{{
% 基準LPFの次数
Wp = 20/50;
Ws = 30/50;
Rp = 2; % 2dB
Rs = 20; % 20dB
[n, Wc] = buttord(Wp, Ws, Rp, Rs)
}}
結果は、
#pre{{
n = 5
Wc = 0.41636
}}
のようにでます。
これで、次数を5とすることにします。
*** バタワースのLPFの作成 [#b442c271]
プリワーピングで20Hzに対応するデジタルのカットオフ周波数...
#pre{{
% バタワースの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の周波数応答は以下の様になります。
&ref(fig4.png);
*** バタワースのHPFの作成 [#lc0625a4]
バタワースのHPFは、sftransのstopをtrueにして周波数変換し...
#pre{{
% ハイパスデジタルフィルタへの変換
% 基準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の周波数応答は以下の様になります。
&ref(fig5.png);
*** バタワースのBPFの作成 [#se62df13]
バタワースのBPFは、Wpfp, Wphpとstop=falseにして周波数変換...
#pre{{
% バンドパスフィルタへの変換
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の周波数応答は以下の様になります。
&ref(fig6.png);
*** バタワースのBRFの作成 [#xd0cd078]
バタワースのBRFは、BPFの設定に対してstop=trueにして周波数...
#pre{{
% バンドエリミネート・フィルタへの変換
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の周波数応答は以下の様になります。
&ref(fig7.png);
*** フィルタの縦結合 [#t84a2e9c]
求めたデジタル・フィルタを結合して複雑なフィルタを作成す...
20Hzと30Hzを遮断するフィルタを作成してみます。
#pre{{
% フィルターの縦結合
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); % プロット
}}
&ref(fig8.png);
** 任意の周波数特性を持つフィルタ [#q5d835f6]
octaveには、周波数特性を与えて、その伝達関数を作成するた...
これを使って最初の1次LPFを再現してみます。
#pre{{
% 目的の周波数特性を作成
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とした場合
&ref(fig9.png);
次数を20とした場合
#pre{{
% 次数を20とした場合
[B, A] = invfreq(H, F, 20, 20);
MyFreqz(B, A, Ft);
}}
&ref(fig10.png);
如何でしょう、最初のアナログフィルタのカーブがよく表現で...
このようにoctaveを使うとIIRのフィルタがとても簡単に作成で...
是非、自分でoctaveを動かしてみてください。
** コメント [#gc9aeafb]
#vote(おもしろかった[37],そうでもない[1],わかりずらい[3])
皆様のご意見、ご希望をお待ちしております。
- バンドパスフィルタを試そうとしたらこのようなエラーが出...
#comment_kcaptcha
終了行:
[[FrontPage]]
#contents
2012/03/17からのアクセス回数 &counter;
* テキストのグラフを再現してみよう [#ec5c600a]
教科書や雑誌にでているグラフを自分で再現してみることで、...
ここでは、octaveを使って雑誌インターフェース2006/09号の８...
** アナログ・フィルタからデジタル・フィルタへの変換 [#k91...
８章IIRフィルタの伝達関数の設計方法にある図１３((Interfac...
&ref(fig13.png);
アナログ・フィルタからデジタル・フィルタへの変換には、以...
$$
s = \frac{2}{T}\frac{1 - z^{-1}}{1 + z^{-1}}
$$
しかし、この変換のみだと図１３の(b)のようにカットオフ周波...
これは、アナログ周波数\( \omega_A \) と デジタル周波数\(\...
&ref(fig11.png);
そこで、このひずみを考慮したカットオフ周波数を計算し（プ...
プリワーピングの式は、以下の様になります。
$$
\omega_D = \frac{2}{T} tan(\omega_A T / 2)
$$
** 周波数応答関数の修正 [#j59ffee6]
octaveで提供されている周波数応答表示関数freqs, freqzでは...
#pre{{
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次ローパス・フィルタの仕様 [#wc3cea2a]
図１３で使っている1次ローパス・フィルタは、以下の様な仕様...
- サンプリング周波数 \(F_T\)は、100 Hz
- カットオフ周波数\( F_C \)は、30 Hz
でアナログ・フィルタの伝達関数は以下の様になります。
$$
G(s) = \frac{1}{ 1 + \frac{s}{2 \pi F_C}}
$$
これをoctaveで以下の様に計算します。
#pre{{
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);
}}
&ref(fig1.png);
仕様通り30Hzで振幅が0.7に減衰しています。
次にこのアナログ・フィルタをs-z変換して周波数特性を表示し...
#pre{{
% そのまま双一次s-z変換をすると
[Zb, Za] = bilinear(b, a, 1/Ft);
% 変換結果の周波数特性を表示
MyFreqz(Zb, Za, Ft);
}}
&ref(fig2.png);
図１３の(b)のようにカットオフ周波数が25Hzにずれてしまいま...
そこで、プリワーピングを使って \(\omega_D\)を新しいカット...
$$
\omega_D = \frac{2}{T} tan(\omega_A T /2)
$$
#pre{{
% プリワーピイングを行う
Wd = 2/T * tan(Wa*T/2);
a = [1/Wd 1]; b = 1;
[Zb, Za] = bilinear(b, a, 1/Ft);
% プリワーピイングを行った結果の周波数特性を表示
MyFreqz(Zb, Za, Ft);
}}
&ref(fig3.png);
求める図１３の(c)の周波数特性を持つ伝達関数ができました。
** 本格的なフィルタの作成 [#l50fb268]
アナログ・フィルタからデジタル・フィルタへの変換方法が分...
以下の手順で進めます。
- カットオフ周波数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の次数 [#a7451550]
それではバタワースフィルタを使って、
- LPF（ローパス・フィルタ）
- HPF（ハイパス・フィルタ）
- BPF（バンドパス・フィルタ）
- BRF（バンドリジェクト・フィルタ）
を作っていきます。
求めるフィルタは、図１９のようにします。((Interface2006/0...
&ref(fig19.png);
ここで問題になるのは、バタワースフィルタの次数をどうする...
octaveには、buttordという関数が用意されており、これを使っ...
#pre{{
% 基準LPFの次数
Wp = 20/50;
Ws = 30/50;
Rp = 2; % 2dB
Rs = 20; % 20dB
[n, Wc] = buttord(Wp, Ws, Rp, Rs)
}}
結果は、
#pre{{
n = 5
Wc = 0.41636
}}
のようにでます。
これで、次数を5とすることにします。
*** バタワースのLPFの作成 [#b442c271]
プリワーピングで20Hzに対応するデジタルのカットオフ周波数...
#pre{{
% バタワースの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の周波数応答は以下の様になります。
&ref(fig4.png);
*** バタワースのHPFの作成 [#lc0625a4]
バタワースのHPFは、sftransのstopをtrueにして周波数変換し...
#pre{{
% ハイパスデジタルフィルタへの変換
% 基準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の周波数応答は以下の様になります。
&ref(fig5.png);
*** バタワースのBPFの作成 [#se62df13]
バタワースのBPFは、Wpfp, Wphpとstop=falseにして周波数変換...
#pre{{
% バンドパスフィルタへの変換
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の周波数応答は以下の様になります。
&ref(fig6.png);
*** バタワースのBRFの作成 [#xd0cd078]
バタワースのBRFは、BPFの設定に対してstop=trueにして周波数...
#pre{{
% バンドエリミネート・フィルタへの変換
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の周波数応答は以下の様になります。
&ref(fig7.png);
*** フィルタの縦結合 [#t84a2e9c]
求めたデジタル・フィルタを結合して複雑なフィルタを作成す...
20Hzと30Hzを遮断するフィルタを作成してみます。
#pre{{
% フィルターの縦結合
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); % プロット
}}
&ref(fig8.png);
** 任意の周波数特性を持つフィルタ [#q5d835f6]
octaveには、周波数特性を与えて、その伝達関数を作成するた...
これを使って最初の1次LPFを再現してみます。
#pre{{
% 目的の周波数特性を作成
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とした場合
&ref(fig9.png);
次数を20とした場合
#pre{{
% 次数を20とした場合
[B, A] = invfreq(H, F, 20, 20);
MyFreqz(B, A, Ft);
}}
&ref(fig10.png);
如何でしょう、最初のアナログフィルタのカーブがよく表現で...
このようにoctaveを使うとIIRのフィルタがとても簡単に作成で...
是非、自分でoctaveを動かしてみてください。
** コメント [#gc9aeafb]
#vote(おもしろかった[37],そうでもない[1],わかりずらい[3])
皆様のご意見、ご希望をお待ちしております。
- バンドパスフィルタを試そうとしたらこのようなエラーが出...
#comment_kcaptcha
ページ名:
SmartDoc
![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
