アルゴリズム (IIRフィルタ)

IIR(Infinite Impulse Response、無限インパルス応答)フィルタは次の図のように、フィードバックを行うデジタルフィルタです。

デジタルフィルタはよく、差を表す数式で表記され、入力シグナルと出力シグナルの関係を定義します。

$b_0*x(n)+b_1*x(n-1)+ \cdots +b_N*x(n-N)-a_0*y(n)-a_1*y(n-1)- \cdots -a_M*y(n-M) = 0$

このように変形します。

$y(n) = \frac{1}{a_0} \left ( b_0*x(n)+b_1*x(n-1)+ \cdots +b_N*x(n-N)-a_1*y(n-1)-a_2*y(n-2)- \cdots -a_M*y(n-M) \right ) = \frac{1}{a_0} \left ( \sum_{i=0}^N b_i*x(n-i) - \sum_{j=1}^M a_j*y(n-j) \right )$

$a_0 \ne 0$ の中で、Nはフィードフォワードフィルタの順番、 $b_i$ はフィードフォワードの係数、Mはフィードバックフィルタの順番、 $a_i$ はフィードバックの係数、x(n)は入力シグナル、y(n)は出力シグナルを表します。 $\sum_{j=1}^M a_j*y(n-j)$ の表記がフィードバックを表します。

デジタルフィルタ表現

OriginではIIRフィルタの出力表現が4種類あります。

変換関数
IIRフィルタの変換関数はz-領域での2つの多数項の複合式 $z - 1$ で表わされます。IIRフィルタの変換フィルタを見つけるには、上記数式を以下のように並び替えます。

$\sum_{i=0}^N b_i*x(n-i) = \sum_{j=0}^M a_j*y(n-j)$

z-変換は次のように定義されます。

$X(z) = \sum_{}^{} x(n)z^n$

z-変換を数式の両辺から取り除くと次のようになります。

$\sum_{i=0}^N b_i*z^{-i}*X(z) = \sum_{j=0}^M a_j*z^{-j}*Y(z)$

IIRフィルタの変換関数をz-領域内で表わすと、次のようになります。

$H(z) = \frac{Y(z)}{X(z)} = \frac{\sum_{i=0}^N b_iz^{-i}}{\sum_{j=0}^M a_jz^{-j}}$
ゼロ-極-ゲイン
上記変換関数が示すように、分子はゼロの位置を示して分母は極の位置を示しています。つまり、変換関数はゼロ-極-ゲイン形式に書き直す事ができます。

$H(z) = \frac{q(z)}{p(z)} = k*\frac{(z-q_1)(z-q_2) \cdots (z-q_N)}{(z-p_1)(z-p_2) \cdots (z-p_M)} = k* \frac{\prod_{i=1}^N (z-q_i)}{\prod_{j=1}^M (z-p_j)}$

ここで、変換関数のkはゲインを、 $q_i$ と $p_j$ はゼロと極をそれぞれ表しています。
状態空間
フィルタシステムの状態空間は次のように定義されます。

$x(n+1) = Ax(n)+Bu(n)$

$y(n) = Cx(n)+Du(n)$
ここで、u(n)は入力、x(n)は状態ベクトル、y(n)は出力、Aはm×mの行列、mはフィルタの順番、Bは列ベクトル、Cは行ベクトル、Dはスケーラを示しています。
第二順序セクション(Second Order Section, SOS)
デジタルフィルタ変換関数に対応する第二順序セクションは以下のように書き表されます。

$H(z) = g \prod_{k=1}^L H_k(z) = g \prod_{k=1}^L \frac{b_{0k}+b_{1k}*z^{-1}+b_{2k}*z^{-2}}{a_{0k}+a_{1k}*z^{-1}+a_{2k}*z^{-2}}$

ここで、gはゲイン、 $b_{0k}, b_{1k}, b_{2k}$ は分子係数、 $a_{0k}, a_{1k}, a_{2k}$ は分母係数を示しています。もしフィルタ順序mが偶数なら $L = m/2$ に、mが奇数なら $L = (m+1)/2$ になります。すると、SOS は以下のようなL×6行列で表記できます。

$SOS = \begin{bmatrix} b_{01} & b_{11} & b_{21} & a_{01} & a_{11} & a_{21} \\ b_{02} & b_{12} & b_{22} & a_{02} & a_{12} & a_{22} \\ b_{03} & b_{13} & b_{23} & a_{03} & a_{13} & a_{23} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ b_{0L} & b_{1L} & b_{2L} & a_{0L} & a_{1L} & a_{2L} \\ \end{bmatrix}$

IIRフィルタを設計するには

IIRフィルタを設計する一般的な手順は次の通りです。

まず、フィルタの仕様を指定します。
ローパスアナログフィルタのプロトタイプを指定します。OriginでサポートされているプロトタイプはButterworth、Chebyshev Type I、Chebyshev Type II、Ellipticです。

手法大きさの二乗の反応関数アナログフィルタの転嫁関数

Butterworth $|H_a(j \Omega)|^2 = \frac{1}{1+ \Omega^{2N}}$ $H_a(s) = \frac{q(s)}{p(s)} = \frac{g}{(s-p_1)(s-p_2)\cdots (s-p_N)} = \frac{g}{\prod_{k=1}^N (s-p_k)}$

Chebyshev Type I $|H_a(j \Omega)|^2 = \frac{1}{1+ \varepsilon^2 T_N^2(\Omega)}$ $H_a(s) = \frac{q(s)}{p(s)} = \frac{g}{(s-p_1)(s-p_2)\cdots (s-p_N)} = \frac{g}{\prod_{k=1}^N (s-p_k)}$

Chebyshev Type II $|H_a(j \Omega)|^2 = \frac{1}{1+ (\varepsilon^2 T_N^2(\frac{1}{\Omega}))^{-1}} = \frac{\varepsilon^2 T_N^2(\frac{1}{\Omega})}{1+ \varepsilon^2 T_N^2(\frac{1}{\Omega})}$ $H_a(s) = \frac{q(s)}{p(s)} = g \frac{(s-q_1)(s-q_2)\cdots (s-q_N)}{(s-p_1)(s-p_2)\cdots (s-p_N)} = g \prod_{k=1}^N \frac{(s-q_k)}{(s-p_k)}$

Elliptic $|H_a(j \Omega)|^2 = \frac{1}{1+ \varepsilon^2 U_N^2(\Omega)}$ $H_a(s) = \frac{q(s)}{p(s)} = g \frac{(s-q_1)(s-q_2)\cdots (s-q_N)}{(s-p_1)(s-p_2)\cdots (s-p_N)} = g \frac{\prod_{i=1}^N (s-q_i)}{\prod_{j=1}^N (s-p_j)}$

上記表で、 $\Omega$ は周波数、Nはフィルタの順番、 $\varepsilon$ は通過帯域の周波数反応内での最大振幅、 $T_N$ はChebyshev多項式、 $U_N$ はJacobian楕円関数、gはスケーラゲイン、sはラプラス変換の面、 $q_k$ または $q_i$ はゼロ、そして $p_k$ または $p_j$ は極を表します。
アナログフィルタのための周波数変換
ローパスフィルタをハイパス、バンドパス、バンドストップフィルタに変換するには希望のカットオフ周波数を使用します。Originでは、状態空間法が周波数変換計算において使用されます。ローパスフィルタの元の変換関数は $H(s')$ として仮定され、返還後の変換関数は $H(s)$ になります。
- ローパスからローパスへは、アナログローパスフィルタをカットオフ周波数1ラジアン/秒でローパスフィルタに流し込み、そのうえで特定のカットオフ周波数を使用して変換します。
  
  $H(s)= H(s')|_{s'=s/\omega_0}$
- ローパスからハイパス
  
  $H(s) = H(s')|_{s'=\omega_0/s}$
- ローパスからバンドパス
  
  $H(s) = H(s')|_{s'= \frac{\omega_0}{B_{\omega}} \frac{(s/\omega_0)^2+1}{s/\omega_0}}$
- ローパスからバンドストップ
  
  $H(s) = H(s')|_{s'= \frac{B_{\omega}}{\omega_0} \frac{s/\omega_0}{(s/\omega_0)^2+1}}$
ここで、 $\omega_0 = sqrt(\omega_1*\omega_2)$ は中央の周波数、 $B_\omega = \omega_2 - \omega_1$ はバンド幅、 $\omega_1$ と $\omega_2$ は下側と上側のハンド幅の角になります。
アナログフィルタをデジタルフィルタに変換します。
アナログフィルタをデジタルフィルタに変換するのに、Originは次式で定義される双線形変換で定義します。

$s = \frac{1-z^{-1}}{1+z^{-1}}$

手法	大きさの二乗の反応関数	アナログフィルタの転嫁関数
Butterworth	$\|H_a(j \Omega)\|^2 = \frac{1}{1+ \Omega^{2N}}$	$H_a(s) = \frac{q(s)}{p(s)} = \frac{g}{(s-p_1)(s-p_2)\cdots (s-p_N)} = \frac{g}{\prod_{k=1}^N (s-p_k)}$
Chebyshev Type I	$\|H_a(j \Omega)\|^2 = \frac{1}{1+ \varepsilon^2 T_N^2(\Omega)}$	$H_a(s) = \frac{q(s)}{p(s)} = \frac{g}{(s-p_1)(s-p_2)\cdots (s-p_N)} = \frac{g}{\prod_{k=1}^N (s-p_k)}$
Chebyshev Type II	$\|H_a(j \Omega)\|^2 = \frac{1}{1+ (\varepsilon^2 T_N^2(\frac{1}{\Omega}))^{-1}} = \frac{\varepsilon^2 T_N^2(\frac{1}{\Omega})}{1+ \varepsilon^2 T_N^2(\frac{1}{\Omega})}$	$H_a(s) = \frac{q(s)}{p(s)} = g \frac{(s-q_1)(s-q_2)\cdots (s-q_N)}{(s-p_1)(s-p_2)\cdots (s-p_N)} = g \prod_{k=1}^N \frac{(s-q_k)}{(s-p_k)}$
Elliptic	$\|H_a(j \Omega)\|^2 = \frac{1}{1+ \varepsilon^2 U_N^2(\Omega)}$	$H_a(s) = \frac{q(s)}{p(s)} = g \frac{(s-q_1)(s-q_2)\cdots (s-q_N)}{(s-p_1)(s-p_2)\cdots (s-p_N)} = g \frac{\prod_{i=1}^N (s-q_i)}{\prod_{j=1}^N (s-p_j)}$