双一次変換によるIIRフィルタ設計

カットオフ周波数を決めれば、単純な計算だけでIIRフィルタの係数が求まる方法はないかと思っていたら、双一次変換と呼ばれる方法があったので紹介します。

今回の記事は以下の本を主に参考にしています。

双一次変換によるIIRフィルタ設計手順

双一次変換によるIIRフィルタの設計手順は以下です。

  1. アナログ基準ローバスフィルタを設計する。
  2. 周波数変換で所望のカットオフ周波数 \omega_c にする。
  3. アナログフィルタに双一次変換を行い、伝達関数  H(z) を求める。

今回は、手始めに以下のような2次のIIRフィルタ伝達関数を作ります。

\displaystyle{
H(z) = \frac{b_0+b_1z^{-1}+b_2z^{-2}}{a_0+a_1z^{-1}+a_2z^{-2}}
}

アナログ基準ローバスフィルタの設計

カットオフ角周波数が \omega_c=1 であるローバスフィルタ(LPF)を基準LPFと呼び、最初にアナログ基準LPFの伝達関数を設計します。

ちなみに、デジタル信号処理の世界では暗黙のうちに、周波数 f' をサンプリング周波数 f_s で割った正規化周波数  f'/f_s を周波数  f と呼びます。同様に、角周波数についても正規化角周波数  2\pi f'/f_s を角周波数  \omega と呼びます。

アナログフィルタの設計では、 s=j\omega を代入すると以下の理想LPFの周波数特性を近似するような有理関数 H(s) を求めます。

\displaystyle{
H_{id}(\omega) = \left\{
\begin{array}{clc}
1  & \omega \lt 1 \\
0  & \omega \geq 1
\end{array}
\right.
}

有理関数というのは以下のように表せる関数です。

\displaystyle{
H(s) = \frac{b_0s^m+b_1s^{m-1}+\cdots+b_m}{s^{n}+a_1s^{n-1}+\cdots + a_n}
}

直接周波数特性を近似するのは難しい(周波数特性は複素数になる)ので、まず理想の2乗振幅特性を有理関数で近似します。アナログフィルタの代表であるバタワースフィルタでは以下の |H_n(\omega)|^2で近似させています。

\displaystyle{
|H_n(\omega)|^2 = \frac{1}{1+\omega^{2n}}
}

図1をみると、確かに理想LPFに近い特性となっています。

図1:バタワースフィルタの周波数特性
図1:バタワースフィルタの周波数特性

2乗振幅特性は有理関数で表せたので、つぎに伝達関数 H(s) を求めていきます。

 |H_n(\omega)|^2平方根を求めて \omega=s/j と置き換えればいいと思いますが、有理関数にならないので工夫が必要です。

 |H_n(\omega)|^2 が以下のように表せることを利用します。

\displaystyle{
\begin{align}
|H_n(\omega)|^2 &= H_n(\omega) H_n^{*}(\omega) \\
&= H_n(\omega) H_n(-\omega)
\end{align}
}

 H_n^{*}(\omega)=H_n(-\omega) は実数をフーリエ変換したときの性質となっています。この性質がないと周波数特性を逆フーリエ変換した結果が複素数になってしまいます。

 |H_n(\omega)|^2 \omega=s/j とおくと、

\displaystyle{
\begin{align}
|H(s)|^2 &= \frac{1}{1+(s/j)^{2n}}  \\
&= \frac{1}{1+(-1)^n s^{2n}}
\end{align}
}

 |H(\omega)|^2 の極を求めてみると、

n が奇数のときは

\displaystyle{
\begin{align}
 1 - s^{2n} &= 0  \\
s &= e^{j2\pi k/2n}
\end{align}
}

n が偶数のときは

\displaystyle{
\begin{align}
 1 + s^{2n} &= 0  \\
s &= e^{j(2k+1)\pi/2n}
\end{align}
}

よって、n=2、n=3のときの |H(s)|^2 の極の位置は図2のようになります。

図2:極の位置
図2:極の位置

 |H(s)|^2=H(s)H(-s) という条件があるので、この条件を満たすように  H(s) の極を選ぶ必要があります。加えて、安定なフィルタを設計するには s 平面の左半平面に極が位置する必要があるので、 H(s) の極は左半平面にある極全部となります。

したがって、n=2のときの伝達関数 H_2(s) は以下のようになります。

\displaystyle{
\begin{align}
 H_2(s) &= \frac{1}{(s-e^{j3\pi/4})(s-e^{j5\pi/4})}  \\
  &= \frac{1}{s^2+\sqrt{2}s + 1} 
\end{align}
}

n=1 のときの伝達関数  H_1(s) も求めてみると以下のようになります。

\displaystyle{
H_1(s) = \frac{1}{s+1}
}

これで、 アナログ基準 LPF の設計ができました。

周波数変換

フィルタのタイプは図3のように4つあります。

図3:フィルタのタイプ
図3:フィルタのタイプ

カットオフ角周波数が  \omega_c である各タイプのフィルタを設計したい場合は表のように s を変換すればよいです。

タイプ 変換
LPF  {\displaystyle s\rightarrow \frac{s}{\omega_c}}
HPF  {\displaystyle s\rightarrow \frac{\omega_c}{s}}
BPF  {\displaystyle s\rightarrow \frac{s^2+\omega_o^2}{s\omega_b}}
BRF  {\displaystyle s\rightarrow \frac{s\omega_b}{s^2+\omega_o^2}}

ここで、 \omega_b=\omega_{c2}-\omega_{c1} です。

H_2 (s) を以上のLPF変換、HPF変換すると以下のようになります。

\displaystyle{
\begin{align}
 H_{LPF}(s) &= \frac{\omega_c^2}{s^2+\sqrt{2}\omega_c s + \omega_c^2}   \\
 H_{HPF}(s) &= \frac{s^2}{s^2+\sqrt{2}\omega_c s + \omega_c^2} 
\end{align}
}

また、H_1 (s) を以上のBPF変換、BRF変換すると以下のようになります。

\displaystyle{
\begin{align}
 H_{BPF}(s) &= \frac{\omega_b s}{s^2+\omega_b s + \omega_o^2}   \\
 H_{BRF}(s) &= \frac{s^2+\omega_o^2}{s^2+\omega_b s + \omega_o^2}  
\end{align}
}

双一次変換

双一次変換は s 変数を以下の式で z 変数にする変換をいいます。

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

このとき、図4のようにアナログ角周波数  \omega_a : [-\infty, \infty] と デジタル角周波数  \omega_d : [-\pi, \pi] が対応するようになっています。

図4:双一次変換による周波数の対応
図4:双一次変換による周波数の対応

双一次変換の式に、 s=j\omega_a z=e^{j\omega_d} を代入すると関係式が以下のようにわかります。

\displaystyle{
\omega_a=2 \tan{\frac{\omega_d}{2}}
}

バタワースフィルタの伝達関数を双一次変換することで、 H(z) を得ることができます。ただし、双一次変換によりカットオフ周波数の位置もずれてしまうため、アナログフィルタ設計の時点で予め考慮する必要があります。

 \omega_a \omega_d の関係式より、デジタルフィルタのカットオフ角周波数を  \omega_c としたいときは、アナログフィルタの設計時ではカットオフ角周波数を以下の \omega'_c にすればよいです。

\displaystyle{
\omega'_c=2 \tan{\frac{\omega_c}{2}}
}

このことはプリワーピングと呼ばれています。

フィルタ係数の決定

双一次変換をバタワースフィルタの伝達関数に行った結果、各フィルタ係数は以下のようになります。 \omega_c はプリワーピングしたカットオフ角周波数です。

表:LPFのフィルタ係数
 b_0  \omega_c^2
 b_1  2\omega_c^2
 b_2  \omega_c^2
 a_0  \omega_c^2+2\sqrt{2}\omega_c+4
 a_1  2\omega_c^2-8
 a_2  \omega_c^2-2\sqrt{2}\omega_c+4
表:HPFのフィルタ係数
 b_0  4
 b_1  -8
 b_2  4
 a_0  \omega_c^2+2\sqrt{2}\omega_c+4
 a_1  2\omega_c^2-8
 a_2  \omega_c^2-2\sqrt{2}\omega_c+4
表:BPFのフィルタ係数
 b_0  2\omega_b
 b_1  0
 b_2  -2\omega_b
 a_0  \hspace{15px} \omega_o^2+2\omega_b+4 \hspace{15px}
 a_1  2\omega_o^2-8
 a_2  \omega_o^2-2\omega_b+4
表:BRFのフィルタ係数
 b_0  \omega_o^2+4
 b_1  2\omega_o^2-8
 b_2  \omega_o^2+4
 a_0  \hspace{15px} \omega_o^2+2\omega_b+4 \hspace{15px}
 a_1  2\omega_o^2-8
 a_2  \omega_o^2-2\omega_b+4

プログラム

Python で記述したIIRフィルタの係数と周波数特性を求めるプログラムは以下です。

import numpy as np
import scipy.signal as sg
import matplotlib.pyplot as plt
import sys

# コマンドラインの引数を取得
args = sys.argv
if len(args)!=4:
    print("\nusage: python main.py  fs  fc  type")
    print("   fs        : sampling frequency")
    print("   fc        : cutoff frequency")
    print("   type      : filter type (lpf, hpf, bpf, brf)")
    raise Exception("Argument error ")

fs = int(args[1])    # サンプリング周波数
fc = float(args[2])  # カットオフ周波数
ftype = args[3]      # フィルタの種類 LPF, HPF, BPF, BRF
fb = 5000            # バンドパスの幅

omega_b  = 2*np.pi*(fb/fs)  # バンドパスの幅
omega_c  = 2*np.pi*(fc/fs)    # カットオフ正規化角周波数
omega_c  = 2*np.tan(0.5*omega_c)  # プリワーピング
omega_c1 = omega_c-0.5*omega_b
omega_c2 = omega_c+0.5*omega_b
omega_o  = np.sqrt(omega_c1*omega_c2) 
a = [0,0,0]
b = [0,0,0] 

# フィルタ係数決定
if(ftype=="lpf"):
    b[0] = omega_c**2
    b[1] = 2*omega_c**2
    b[2] = omega_c**2
    a[0] = omega_c**2+2*np.sqrt(2)*omega_c+4
    a[1] = 2*omega_c**2-8
    a[2] = omega_c**2-2*np.sqrt(2)*omega_c+4
elif(ftype=="hpf"):
    b[0] = 4
    b[1] = -8
    b[2] = 4
    a[0] = omega_c**2+2*np.sqrt(2)*omega_c+4
    a[1] = 2*omega_c**2-8
    a[2] = omega_c**2-2*np.sqrt(2)*omega_c+4
elif(ftype=="bpf"):
    b[0] = 2*omega_b
    b[1] = 0
    b[2] = -2*omega_b
    a[0] = omega_o**2+2*omega_b+4
    a[1] = 2*omega_o**2-8
    a[2] = omega_o**2-2*omega_b+4
else: # brf
    b[0] = omega_o**2+4
    b[1] = 2*omega_o**2-8
    b[2] = omega_o**2+4
    a[0] = omega_o**2+2*omega_b+4
    a[1] = 2*omega_o**2-8
    a[2] = omega_o**2-2*omega_b+4

# 周波数特性の出力
w, h = sg.freqz(b, a)
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(w, np.abs(h), c="red")
ax.set_ylim([0,1.2])
ax.set_xlabel("omega") 
ax.set_ylabel("Amplitude [dB]")
fig.savefig(ftype+".png")

以下のようにコマンドを実行することで、サンプリング周波数 fs=48kHz、カットオフ周波数 fc=12kHz の LPF のフィルタ係数と振幅特性が求まります。BPF や BRF を指定したときはバンド幅は 5kHz で、第二引数の値がバンドの中心となります。

  python main.py 48000 12000 lpf

周波数特性について

サンプリング周波数 fs=48kHz、カットオフ周波数 fc=12kHz を指定してさきほどのプログラムを実行した結果、図5のような特性になりました。

図5:作成したフィルタの振幅特性
図5:作成したフィルタの振幅特性

あまり急峻な特性ではないですが、2次の IIR フィルタなので仕方ないですかね。

白色雑音に対して各フィルタを畳み込んだときのスペクトログラムは図6のようになります。

図6:フィルタを畳み込んだ結果
図6:フィルタを畳み込んだ結果

やはり急峻ではないですが、上手く機能しているのがわかると思います。

おわりに

長くなりましたが、双一次変換によるIIRフィルタ設計についてまとめました。FPGAでフィルタを作成するときなどに利用しようと思います。

参考文献

  1. 陶山 健仁、”ディジタルフィルタ原理と設計法”、科学情報出版株式会社、2018.
  2. 尾知 博、”ディジタル・フィルタ設計入門 各種フィルタの原理とDSPによる実現”、CQ出版社、1990.

お問い合わせフォーム プライバシーポリシー

© 2021 Setoti All rights reserved.