離散時間システムモデル

離散時間システムモデルは、デジタルフィルターの表現方法です。MATLAB^® 技術計算環境では、以下の節で説明するいくつかの離散時間システムモデルがサポートされています。

伝達関数

"伝達関数" は、デジタルフィルターの基本的な Z 領域表現で、フィルターを 2 つの多項式の比として表現します。この表現は、このツールボックスの主な離散時間モデルです。デジタルフィルターの差分方程式の Z 変換に対する伝達関数モデルは、以下のとおりです。

$Y (z) = \frac{b (1) + b (2) z^{- 1} + \dots + b (n + 1) z^{- n}}{a (1) + a (2) z^{- 1} + \dots + a (m + 1) z^{- m}} X (z) .$

ここで、定数 b(i) および a(i) はフィルター係数であり、フィルターの次数は n と m の最大値になります。MATLAB 環境では、これらの係数を 2 つのベクトル (通常は行ベクトル) に格納します。一方の行ベクトルには分子が、もう一方には分母が格納されます。伝達関数形式の詳細については、フィルターと伝達関数を参照してください。

零点-極-ゲイン

伝達関数の因数分解型または "零点-極-ゲイン" 型は以下のとおりです。

$H (z) = \frac{q (z)}{p (z)} = k \frac{(z - q (1)) (z - q (2)) ... (z - q (n))}{(z - p (1)) (z - p (2)) ... (z - p (n))} .$

通常、多項式の係数は行ベクトルに、多項式の根は列ベクトルに格納されます。つまり、零点-極-ゲイン形式では、伝達関数の分子と分母の零点および極の位置は列ベクトルで表されます。因数分解した伝達関数のゲイン k は、MATLAB のスカラーとなります。

関数 poly および roots により、多項式表現と零点-極-ゲイン表現間の変換が行われます。たとえば、次の簡単な IIR フィルターを考えてみましょう。

b = [2 3 4];
a = [1 3 3 1];

このフィルターの零点と極は、次のようになります。

q = roots(b)
p = roots(a)
% Gain factor
k = b(1)/a(1)

元の多項式へと戻すと、次のようになります。

bb = k*poly(q)
aa = poly(p)

ここで、b と a は次の伝達関数を表しています。

$H (z) = \frac{2 + 3 z^{- 1} + 4 z^{- 2}}{1 + 3 z^{- 1} + 3 z^{- 2} + z^{- 3}} = \frac{z (2 z^{2} + 3 z + 4)}{z^{3} + 3 z^{2} + 3 z + 1} .$

b = [2 3 4] に対し、関数 roots では、z が 0 となる零点が捉えられません。実際、入力となる伝達関数で零点と極の数が等しくない場合、この関数は z が 0 となる零点と極は捉えられません。ほとんどの場合これは問題とはなりませんが、この問題を回避するには、関数 roots を使う前に、0 を追加してベクトルを同じ長さにします。たとえば b = [b 0] のようにします。

状態空間

デジタルフィルター、または連立差分方程式のシステムは、常に、1 階差分方程式の集合として表現することができます。行列形式または "状態空間" 形式では、以下のように方程式を書くことができます。

$\begin{matrix} x (n + 1) = A x (n) + B u (n) \\ y (n) = C x (n) + D u (n), \end{matrix}$

ここで、u は入力、x は状態ベクトル、y は出力です。単一チャネルシステムの場合、A は m 行 m 列の行列 (m はフィルターの次数)、B は列ベクトル、C は行ベクトル、D はスカラーです。状態空間表記法は、入力 u と出力 y がベクトルとなり、B、C、および D が行列となるマルチチャネルシステムで特に便利です。

状態空間表現は、MATLAB 環境に簡単に拡張できます。A、B、C および D は方形配列で、MATLAB 関数では個体変数として扱われます。

状態空間方程式を Z 変換して組み合わせると、状態空間形式と伝達関数形式は等価であることがわかります。

$Y (z) = H (z) U (z), where H (z) = C {(z I - A)}^{- 1} B + D$

線形システムの状態空間表現を熟知していなくても心配する必要はありません。一部のフィルター設計アルゴリズムでは状態空間形式を内部的に使用していますが、これらのアルゴリズムを使用するにあたって状態空間の概念に関する知識は必要ありません。ただし、自分のアプリケーションで状態空間ベースの信号処理を広範に使用する場合は、Control System Toolbox™ 製品に付属の状態空間ツールの総合的なライブラリを参考にしてください。

部分分数展開 (留数型)

伝達関数には、対応する "部分分数展開"、すなわち "留数" 形式の表現が存在します。これは次の式によって与えられます。

$\frac{b (z)}{a (z)} = \frac{r (1)}{1 - p (1) z^{- 1}} + ... + \frac{r (n)}{1 - p (n) z^{- 1}} + k (1) + k (2) z^{- 1} + ... + k (m - n + 1) z^{- (m - n)}$

ただし、H(z) には重複極がないものとします。ここで、n は有理伝達関数 b(z)/a(z) の分母多項式の次数です。r が多重度 s_r の極の場合、H(z) は次の型の項を含むことになります。

$\frac{r (j)}{1 - p (j) z^{- 1}} + \frac{r (j + 1)}{{(1 - p (j) z^{- 1})}^{2}} ... + \frac{r (j + s_{r} - 1)}{{(1 - p (j) z^{- 1})}^{s_{r}}}$

Signal Processing Toolbox™ の関数 residuez では、伝達関数形式と部分分数展開形式間の相互変換が行われます。関数 residuez の末尾の "z" は、"z" 領域、すなわち離散領域を表しています。関数 residuez は、極を列ベクトル p、極に対応する留数を列ベクトル r、また、元の伝達関数の不規則な部分を行ベクトル k で返します。関数 residuez は、2 つの極の大きさの差がいずれか一方の極の大きさの 0.1 パーセント未満であれば、その 2 つは同一の極であるとします。

部分分数展開は、信号処理において伝達関数の逆 Z 変換を求める方法の 1 つとして使用されます。たとえば、次の式を考えてみましょう。

$H (z) = \frac{- 4 + 8 z^{- 1}}{1 + 6 z^{- 1} + 8 z^{- 2}}$

この式の部分分数展開は次のようになります。

b = [-4 8];
a = [1 6 8];
[r,p,k] = residuez(b,a)

これは次の式に相当します。

$H (z) = \frac{- 12}{1 + 4 z^{- 1}} + \frac{8}{1 + 2 z^{- 1}}$

H(z) の逆 Z 変換を求めるには、H(z) の 2 つの加数の逆 Z 変換の和を求めます。因果的インパルス応答は次により与えられます。

$h (n) = - 12 {(- 4)}^{n} + 8 {(- 2)}^{n}, n = 0, 1, 2, \dots$

これを MATLAB 環境で検証するには、次のように入力します。

imp = [1 0 0 0 0];
resptf = filter(b,a,imp)
respres = filter(r(1),[1 -p(1)],imp)+...
  filter(r(2),[1 -p(2)],imp)

2 次セクション (SOS)

いずれの伝達関数 H(z) にも、"2 次セクション" の表現があります。

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

ここで、L はシステムを記述する 2 次セクションの数です。MATLAB 環境では、離散時間システムに対する 2 次セクションを L 行 6 列の配列 sos として表します。配列 sos の各行には 1 つの 2 次セクションが含まれ、2 次セクションを記述する 3 つの分子係数と 3 つの分母係数が行の要素となります。

$sos = (\begin{array}{l} 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_{0 L} & b_{1 L} & b_{2 L} & a_{0 L} & a_{1 L} & a_{2 L} \end{array})$

2 次セクションでフィルターを表現するには多くの方法があります。極と零点を慎重に組み合わせ、カスケードにおけるセクションの順序を決め、セクションを乗法的にスケーリングすることにより、一部の固定小数点フィルターの実装で量子化ノイズのゲインを抑え、オーバーフローを避けることができます。線形システムの変換で説明している関数 zp2sos および ss2sos により、極-零点の組み合わせ、セクションのスケーリング、および順序の設定が行われます。

メモ

Signal Processing Toolbox のすべての 2 次セクション変換は、デジタルフィルターのみに適用されます。

ラティス構造

離散 "N" 次の全極または全零フィルターが、多項式係数 a(n), n = 1, 2, ..., N+1 によって表現される場合、対応する "N" 個のラティス構造係数 k(n), n = 1, 2, ..., N が存在します。パラメーター k(n) は、フィルターの "反射係数" とも呼ばれます。これらの反射係数が与えられれば、以下に示すような離散フィルターを実装することができます。

FIR ラティスフィルターおよび IIR ラティスフィルターの構造図

Block diagrams representing FIR and IIR lattice filters. The FIR filter has both minimum-phase and maximum-phase outputs. The IIR filter has an all-pole output.

多項式係数 a および b で表される一般的な極-零点 IIR フィルターには、分母 a に対するラティス係数 k(n) と、分子 b に対するラダー係数 v(n) が存在します。ラティス/ラダーフィルターは、次のように実装することができます。

ラティス/ラダーフィルターのブロック線図

Block diagram representing a lattice/ladder filter. The filter has an ARMA output.

ツールボックス関数 tf2latc は、多項式形式の FIR フィルターまたは IIR フィルターを入力引数とし、対応する反射係数を返します。次に、多項式型の FIR フィルターの例をあげます。

b = [1.0000   0.6149   0.9899   0.0000   0.0031  -0.0082];

このフィルターのラティス (反射係数) は、次のように表現されます。

k = tf2latc(b)

IIR フィルターでは、反射係数の大きさでフィルターの安定性を容易にチェックすることができます。多項式に対応するすべての反射係数の振幅が 1 未満であれば、多項式の根はすべて単位円の内部にあります。たとえば、上記の分子多項式 b と、次の分母多項式から成る IIR フィルターを考えます。

a = [1 1/2 1/3];

フィルターのラティス表現は以下のようになります。

[k,v] = tf2latc(b,a);

k のすべての反射係数に対し abs(k) < 1 であるため、フィルターは安定しています。

関数 latc2tf により、ラティス (反射) 係数からフィルターの多項式係数が計算されます。反射係数ベクトル k が与えられると、対応する多項式形式は次のようになります。

b = latc2tf(k);

このラティス係数、またはラティス/ラダー係数は、関数 latcfilt を使用したフィルターの実装に使うことができます。

畳み込み行列

信号処理において、2 つのベクトル、または行列の畳み込みを行うことは、一方の入力オペランドを他方でフィルター処理することと同じです。この関係により、デジタルフィルターを "畳み込み行列" として表現することができます。

任意のベクトルを与えると、ツールボックス関数 convmtx により、もう一方のベクトルとの内積が、2 つのベクトルの畳み込みと等しくなる行列が生成されます。生成された行列は、適切な長さの任意のベクトルに適用可能なデジタルフィルターを表しますが、内積を計算するにはオペランドの内部次元が一致している必要があります。

デジタルフィルターの分子係数を表すベクトル b の畳み込み行列は、次のようになります。

b = [1 2 3];
x = randn(3,1);
C = convmtx(b',3);

b と x を畳み込む 2 つの等価な方法は、次のとおりです。

y1 = C*x;
y2 = conv(b,x);