線形時不変システムの差分方程式

線形時不変システムの差分方程式( linear, time-invariant difference equation )は y = filter (B, A, X) 関数で解くことができます。差分方程式が解けるとデジタルフィルターの実験が簡単にできます。filter で解ける差分方程式は次のような形をしています。下の差分方程式の y(k) の係数 a(k) からなるベクトルが filter 関数の A、x(k) の係数 b(k) のベクトルが B、(関数の引数の順序は差分方程式とは逆なので注意)入力 x(k) の時系列を要素とするベクトルが X になります

 N                   M
SUN a(k+1) y(n-k) = SUM b(k+1) x(n-k)     for 1<=n<=length(x)
k=0                 k=0

ただし、N=length(a)-1、M=length(b)-1 です。これと次の差分方程式は等価です。

          N                   M
y(n) = - SUM c(k+1) y(n-k) + SUM d(k+1) x(n-k)     for i<=n<=length(x)
         k=1                 k=0

ここで、 c = a/a(1)、d = b/a(1) です。

z 変換の用語で言うと、y は 離散時間のシグナル x が次の rational system function を通過したときの出力となります。

           M
          SUM d(k+1) z^(-k)
          k=0
H(z) = ------------------------
              N
         1 + SUM c(k+1) z(-k)
             k=1

フィルター

それでは次の簡単な差分方程式について実験してみましょう。これは入力の微分を出力する簡単なデジタルフィルターになります。

y(n) = ( x(n) - x(n-1) ) / Δt

オクターブでの出力は次のようになります。

octave:> pi2 = 2 * pi;
octave:> t = linspace (0, pi2, 100);
octave:> X = sin (t);
octave:> A = [1];
octave:> B = [1, -1] / ( pi2 / 100 );
octave:> Y = filter (B, A, X);
octave:> plot (t, X, t, Y)

コンボリューション

線形離散時間システムのインパルス応答関数 h(t) が分かっている場合は、h(t) と入力 x(t) のコンボリューションを計算すると線形システムに x(t) を入力したときの出力 y(t) を計算できます。つまり、y(t) = conv ( h(t), x(t) ) です。Octave での実行例は下のようになります。

0 から 10 秒までの時間を 100 等分した時系列で 正弦波の入力を、指数関数的に減少するインパルス応答をもつシステム a に入力したときの出力 y を計算しています。下の実行例で y1 = y(1:100) で計算結果 y のうちインデックスが 1 から 100 までの値に限定していますが、これは y = conv(a, x) で計算される y のインデックスが 101 以降は x の入力がなくなった後のシステムの出力のデータだからです。

octave:> t = linspace (0, 10, 100);
octave:> x = sin (5 * t);
octave:> a = exp(-t);
octave:> plot (t, x, t, a)
octave:> y = conv (a, x);
octave:> y1 = y(1:100);
octave:> plot (t, y1)