この記事ではラプラス変換の記事の続編として、「離散化したラプラス変換」と呼ばれるZゼッド変換について説明します。

はじめに

以前、ラプラス変換の基本について説明する記事を投稿しました。

その後、ラプラス変換の応用としてよく使われる伝達関数について説明する記事も投稿しました。

このように、ラプラス変換が物理学や工学の問題で使えることはお分かりいただけたかと思います。しかし、ラプラス変換が扱うのは連続値の時間 t ∈ [0, ∞) の空間です。現実世界では時間は連続値ですが、計算がコンピューターで行われる現代では、多くのものが離散化されます。時間も t ∈ { 0, 0.01, 0.02, 0.03, ... } のように微小な時間間隔で扱うのが一般的です。いわゆるディジタル化です。

厳密には時間だけでなく、表現したい位置$x$などの値もディジタル化されますが、それは精度の問題であり、計算方法に影響を与えないため、特に説明する必要はありません。特に複雑になるのは時間の離散化です。

離散化のイメージを極端に表現しているものとして、昔の8ビットゲームなどが挙げられます。

（前回の記事に投稿したラプラスの画像をFLUX.1 Kontextによって8ビット風にしたもの。幼少期初めて遊んだポケモンみたいで懐かしい）

このように時間の値を離散的に扱う場合、ラプラス変換も離散化する必要があります。その離散化されたラプラス変換に相当する概念がZ変換と呼ばれる変換です。

尚、フーリエ変換には離散フーリエ変換がありますが、ラプラス変換には「離散ラプラス変換」というものはありません。Z変換は「離散ラプラス変換」と呼ばれることもありますが、実際には単に「ラプラス変換を離散化したもの」というだけでなく、表現方法も大きく異なるため、フーリエ変換と離散フーリエ変換のような単純な関係ではありません。理解するにはより多くの工夫が必要です。

この記事でのZ変換の解説は、既にラプラス変換の記事を読んだことを前提に進めます。

又、Z変換もラプラス変換と同様に伝達関数を定義することができますが、今回は伝達関数の話には触れません。Z変換の伝達関数については今後の記事で書く予定です。

離散化の考え方

時間$t$の関数である何らかの値$x(t)$を考えます。時間が連続値であれば、あらゆる時間に$x$値の情報があり、グラフを描くと滑らかな曲線になるでしょう。

しかし、ディジタルの世界では一定の間隔で値を認識するため、その時間だけ値が存在します。任意の時間の値を求める場合、基本的には直前に認識した値を採用します。

離散化すると、元の連続時間$t$は間隔$\Delta_T$と整数$n$を用いて次のように表現できます。

t \equiv n\Delta_T

すると、時間の関数である値$x$は次のように時系列として書き換えられます。

x(t) \equiv x(n\Delta_T) \equiv x_n

$n$は整数なので、このように下付きで書く方が一般的ですが、括弧$($　$)$で書いても同じ意味です。又$[$　$]$を使う教科書もありますが、ここでは下付きで表記します。

このように、$x$は離散値$n$の関数になります。

\begin{align}
x(\Delta_T) &&&\equiv& x_1 \\
x(2\Delta_T) &&&\equiv& x_2 \\
x(t-\Delta_T) &\equiv& x((n-1)\Delta_T) &\equiv& x_{n-1} \\
x(t+\Delta_T) &\equiv& x((n+1)\Delta_T) &\equiv& x_{n+1}
\end{align}

では、このような関数の積分を求める場合どうなるでしょうか？連続関数であれば積分はグラフの下の面積ですが、離散化するとこのように並んだ長方形になります。

すると、積分は長方形の面積$x_n\Delta_T$の総和になります。

\int_{0}^{\infty}x(t)dt = \sum_{n=0}^{\infty}x_n\Delta_T

次に、微分について見てみましょう。$\dot{x}(t)$を時間$t$に対する値$x(t)$の微分とします。

\dot{x}(t) = \frac{dx(t)}{dt}

すると、$x(t)$の変化は$\dot{x}(t)$と時間の微小変化$dt$の積で表現できます。

dx(t) = \dot{x}(t)dt

すると、時間が経過した後の値$x(t+dt)$は次のようになります。

\begin{align}
x(t+dt) &= x(t) + dx(t) \\
&= x(t) + \dot{x}(t)dt
\end{align}

経過時間が離散化間隔$\Delta_T$であるとすると、

\begin{align}
x(t+\Delta_T) &= x(t) + \dot{x}(t) \Delta_T \\
x(n\Delta_T+\Delta_T) &= x(n\Delta_T) + \dot{x}(n\Delta_T) \Delta_T \\
x_{n+1} &= x_n + \dot{x}_n \Delta_T
\end{align}

そして、次の関係が得られます。

\begin{align}
x_1 &= x_0 + \dot{x}_0\Delta_T \\
x_2 &= x_1 + \dot{x}_1\Delta_T \\
\cdots \\
x_{n} &= x_{n-1} + \dot{x}_{n-1}\Delta_T \\
\end{align}

微分はこのように求められるとわかります。

\dot{x}_{n} = \frac{x_{n+1} - x_{n}}{\Delta_T}

更に$\dot{x}_n$の微分である$\ddot{x}_n$も同様に、

\begin{align}
\ddot{x}_n &= \frac{\dot{x}_{n+1} - \dot{x}_{n}}{\Delta_T} \\
&= \frac{(x_{n+2} - x_{n+1}) - (x_{n+1} - x_{n})}{\Delta_T^2} \\
&= \frac{x_{n+2} - 2x_{n+1} + x_{n}}{\Delta_T^2} \\
\end{align}

同じ方法で更に進めると、

\dddot{x}_n = \frac{x_{n+3} - 3 x_{n+2} + 3x_{n+1} - x_{n}}{\Delta_T^3}

$m$階微分は次のように表現できることがわかります。

\begin{align}
x_n^{(m)} &= \frac{\sum_{k=0}^{m}(-1)^{k}\frac{m!}{k!(m-k)!}x_{n+m-k}}{\Delta_T^m} \\
&= \frac{x_{n+m} - m x_{n+m-1} + (-1)^{2}\frac{m!}{2!(m-2)!}x_{n+m-2} + \cdots +(-1)^{m-1} x_{n+1} +(-1)^m x_{n}}{\Delta_T^m}
\end{align}

以上説明したのはオイラー法を基にした離散化の考え方ですが、実際の離散化方法には他にもルンゲ–クッタ法や台形法などいくつかあります。通常はオイラー法が一般的に使われていますが、他の方法が使われる場合もあるため、方法は一つではないことを意識しておいた方がよいかもしれません。

漸化式と微分方程式の関係

上述のように、連続時間での微分は離散化時間では前の段階の値$x_{n-1}$、$x_{n-2}$を用いた表現になります。すると、微分方程式は前の時間の値が含まれる方程式になります。例えば次のように、

x_n = a x_{n-1} + b x_{n-2} + c x_{n-3} + \cdots

このような方程式は漸化式ぜんかしき又は差分方程式と呼ばれます。

漸化式と微分方程式は、どちらも時間の関数である値の変化を表現する方程式であり、深い関係があります。

ラプラス変換が微分方程式を解くために使われるのに対し、Z変換は漸化式を解くために使うと言えます。

ラプラス変換からZ変換の導出

それでは、Z変換がどこから来たのか見ていきましょう。まずはZ変換の定義を見る前に、ラプラス変換から始めます。

ラプラス変換は次のような積分で定義されています。

\mathcal{L}[x(t)] = X(s) = \int_{0}^{\infty}x(t)e^{-st}dt

もし時間が離散化され、$t \equiv n\Delta_T$、$x(t) \equiv x_n$となった場合、上述のように積分は総和になります。

\begin{align}
X(s) &= \sum_{n=0}^\infty x_n e^{-sn\Delta_T}\Delta_T \\
\end{align}

ここで、

s = \frac{1}{\Delta_T}\log(z)

z = e^{s\Delta_T}

とすると、

X(s) = \sum_{n=0}^{\infty} x_n z^{-n}\Delta_T

それで$X(s)$は新しい変数$z$の関数$X(z)$として書き換えられます。

X(z) = \frac{X(s)}{\Delta_T} = \sum_{n=0}^{\infty} x_n z^{-n}

これがZ変換の定義です。ラプラス変換が$\mathcal{L}[　]$で表記するのと同様に、Z変換は$\mathcal{Z}[　]$と表記されます。

\mathcal{Z}[x_n] = X(z) = \sum_{n=0}^{\infty} x_n z^{-n}

このように、ラプラス変換に離散化を適用すると離散版のラプラス変換が得られますが、そのまま$s$を使うのではなく、$z$という新しい変数での表現に切り替えます。それがZ変換です。

表現方法は少し変わりますが、実質的にZ変換はラプラス変換を離散化したものであることは間違いありません。

片側Z変換と両側Z変換

補足となりますが、ここに書かれている定義は厳密には片側Z変換と呼ばれるものです。他に両側Z変換があります。

\mathcal{Z}[x_n] = X(z) = \sum_{n=-\infty}^{\infty} x_n z^{-n}

違いは0ではなく$-\infty$から始まる点だけです。単にZ変換と言う場合、どちらも指す可能性があります。

実はラプラス変換にも両側ラプラス変換があります。本来のラプラス変換は片側なので、片側Z変換に相当しますが、両側ラプラス変換を離散化すると両側Z変換になります。

通常、考慮したい問題は0から始まれば十分な場合が多く、殆どの場合は負の側は必要なく片側で十分です。そのため、片側Z変換について気にする必要はなく、今後もZ変換と言えば片側Z変換のみを指します。

Z変換の計算の例

例としてこの関数のZ変換を行います。

x_n = 0.5^n

\begin{align}
\mathcal{Z}[0.5^n] &= \sum_{n=0}^{\infty} 0.5^n z^{-n} \\
&= \sum_{n=0}^{\infty} \left(0.5z^{-1}\right)^n \\
&= \frac{1}{1-0.5z^{-1}}
\end{align}

総和が消えて綺麗な形になったのは、収束する等比数列の形になっているからです。

\begin{align}
\sum_{i=0}^{\infty}r^k = \frac{1}{1-r}
\end{align}

（ただし$r<1$は収束条件）

このように、ラプラス変換と同じように、特定の形の関数はZ変換を行うと綺麗な関数になります。

よく使われる変換は変換表に纏められているため、実際に計算するのではなく、変換表を見て比べることで行うことが多いです。

Z変換の表

よく使う関数の変換表をここに掲載しておきます。ただし、$z$空間では$z$そのままで表現するだけでなく、$z^{-1}$を使った表現もよく使われるため、どちらもすぐ使えるように両方載せます。

n空間	z空間
$\delta(n)$	$1$
$1$	$\begin{align} \frac{1}{1-z^{-1}} \end{align}$	$\begin{align} \frac{z}{z-1} \end{align}$
$a^n$	$\begin{align} \frac{1}{1-az^{-1}} \end{align}$	$\begin{align} \frac{z}{z-a} \end{align}$
$na^n$	$\begin{align} \frac{az^{-1}}{(1-az^{-1})^2} \end{align}$	$\begin{align} \frac{az}{(z-a)^2} \end{align}$
$n^2a^n$	$\begin{align} \frac{az^{-1}(1+az^{-1})}{(1-az^{-1})^3} \end{align}$	$\begin{align} \frac{az(z+a)}{(z-a)^3} \end{align}$
$\cos(\omega_0 n)$	$\begin{align} \frac{ 1-z^{-1} \cos(\omega_0) }{ 1-2z^{-1}\cos(\omega_0)+ z^{-2}} \end{align}$	$\begin{align} \frac{ z^2 - z\cos(\omega_0) }{ z^{2}-2z\cos(\omega_0)+ 1} \end{align}$
$\sin(\omega_0)$	$\begin{align} \frac{ z^{-1} \sin(\omega_0) }{ 1-2z^{-1}\cos(\omega_0)+ z^{-2} } \end{align}$	$\begin{align} \frac{ z \sin(\omega_0) }{ z^2-2z\cos(\omega_0)+ 1 } \end{align}$
$a^n\cos(\omega_0)$	$\begin{align} \frac{ 1-a z^{-1} \cos( \omega_0) }{ 1-2az^{-1}\cos(\omega_0)+ a^2 z^{-2} } \end{align}$	$\begin{align} \frac{ z^2-a z \cos( \omega_0) }{ z^2-2az\cos(\omega_0)+ a } \end{align}$
$a^n\sin(\omega_0)$	$\begin{align} \frac{ az^{-1} \sin(\omega_0) }{ 1-2az^{-1}\cos(\omega_0)+ a^2 z^{-2} } \end{align}$	$\begin{align} \frac{ az \sin(\omega_0) }{ z^2-2az\cos(\omega_0)+ a^2 } \end{align}$
$\begin{align} \frac{1}{n!} \end{align}$	$e^z$
ポアソン分布 $\begin{align} \frac{\lambda^n}{n!}e^{-\lambda} \end{align}$	$e^{\lambda(z-1)}$

線形性

次に、重要な性質を見てみましょう。

ラプラス変換と同様、Z変換にも線形性があります。

\begin{align}
      \mathcal{Z}[ax_n+by_n] &= a\mathcal{Z}[x_n] + b\mathcal{Z}[y_n] \\
      &= aX(z)+bY(z)
\end{align}

シフト性

Z変換において、ラプラス変換の微分と積分に相当するのはシフト性という性質です。

$x_n$のZ変換が$\mathcal{Z}[x_n] = X(z)$であるとき、その前の段階の$x_{n-1}$のZ変換は$X(z)$との関係で次のように表現できます。

\mathcal{Z}[x_{n-1}] = z^{-1}X(z) + x_{-1}

ここで$x_{-1}$は初期値です。$n \ge 0$を考慮しているのに$-1$があるのは違和感があるかもしれませんが、これはあくまで初期値とされている数値であり、単なる定数です。

そして2段階前のZ変換は次のようになります。

\mathcal{Z}[x_{n-2}] = z^{-2}X(z) + z^{-1}x_{-1} + x_{-2}

これは微分のラプラス変換と似た形であることがわかりますね。

そして逆に、次の段階のZ変換は次のようになります。

\mathcal{Z}[x_{n+1}] = zX(z) - x_0z

その次は、

\mathcal{Z}[x_{n+2}] = z^2X(z) - x_0z^2 - x_1z

ここで初期値が出てきますが、伝達関数として使う場合など、初期値が常に0である場合が多いです。その場合、シフトの関係は次のようになります。

\mathcal{Z}[x_{n+k}] = z^kX(z)

ラプラス変換では「微分は$s$を掛ける、積分は$s^{-1}$を掛ける」と同じように、Z変換では「前は$z^{-1}$を掛ける、後ろは$z$を掛ける」です。

運動方程式の離散化とZ変換による解

Z変換の定義を学んだ後、実際の物理学の問題に使う例を挙げます。まず最も簡単な例として、空気抵抗の中で走っている車などの運動方程式を考えましょう。

m\ddot{x}(t) + c \dot{x}(t) = 0

ここで$m$は車の質量、$c$は空気抵抗の係数、$x(t)$は時間の関数である位置です。ただし今回は位置ではなく、その微分である速度$v(t) = \dot{x}(t)$を考えます。すると、次のような簡単な微分方程式になります。

m\dot{v}(t) + c v(t) = 0

ここで上述のオイラー法を用いて漸化式にします。

\begin{align}
m\frac{v_{n+1} - v_{n}}{\Delta_T} + cv_n = 0
\end{align}

そしてこの漸化式を解くために、Z変換を行います。

\begin{align}
\mathcal{Z}\left[m\frac{v_{n+1} - v_{n}}{\Delta_T} + cv_n\right] = 0
\end{align}

線形性とシフト性を使って解いていきます。

\begin{align}
0 &= m\frac{\mathcal{Z}[v_{n+1}] - \mathcal{Z}[v_{n}]}{\Delta_T} + c\mathcal{Z}[v_n] \\
0 &= m\frac{(zV(z) - v_0z) - V(z)}{\Delta_T} + cV(z) \\
0 &= mzV(z)  - mv_0z - mV(z) + c\Delta_TV(z) \\ 
0 &= (m + c\Delta_Tz^{-1} - mz^{-1})V(z) - mv_0 \\
V(z) &= \frac{v_0}{1 + (c\Delta_T/m - 1)z^{-1}}
\end{align}

このように$V(z)$が得られました。しかし、私たちが欲しいのは$z$空間の$V(z)$ではなく、元の$n$で表現された$v_n$です。そのために逆Z変換を行う必要があります。

逆Z変換

ここはラプラス変換と同じ流れですね。ラプラス変換と同様に、直接計算するよりも変換表と照らし合わせる方法が一般的です。

変換表を見ると、次のものがあります。

\mathcal{Z}^{-1}\left[\frac{1}{1-az^{-1}}\right] = a^n

ここで$a = 1 - c\Delta_T/m$です。したがって、

\begin{align}
v_n &= \mathcal{Z}^{-1}\left[V(z)\right] \\
&= v_0\left(1 - \frac{c\Delta_T}{m}\right)^n
\end{align}

これで$v_n$の簡単な表現が得られました。

結果を見ると、これは時間$n$の進行に伴う指数的減衰であることがすぐに読み取れると思います。

確認するために、初期値を設定してPythonでグラフを描いてみましょう。

import numpy as np
import matplotlib.pyplot as plt

v_0 = 30
m = 1000
c = 200
Delta_T = 0.4
n = np.arange(55)
t = n*Delta_T
v = v_0*(1-c*Delta_T/m)**n

plt.figure(figsize=[5,3],dpi=100,facecolor='#f0e6f8')
plt.plot(t,v,'o-',ms=3,lw=0.5)
plt.ylabel(r'$v(t)$',size=14)
plt.xlabel(r'$t$',size=14)
plt.grid(ls=':')
plt.show()

微分方程式のラプラス変換との結果比較

上述の離散化の結果は、あくまで元の連続時間の近似です。理論上、$\Delta_T$が小さいほど本来の連続時間の結果に近づくはずです。

結果の違いを比較するために、元の微分方程式をラプラス変換で解いてみましょう。

\begin{align}
m\dot{v}(t) + cv(t) &= 0 \\
m\mathcal{L}[\dot{v}(t)] + c \mathcal{L}[v(t)] &= 0 \\
msV(s) - mv_0 + cV(s) &= 0 \\
V(s) &= \frac{v_0}{s + c/m} \\
v(t) &= v_0e^{-ct/m}
\end{align}

結果を見ると、これも指数的減衰であることは明らかです。しかし、表現は少し異なり、底は$e$になっています。

実は、指数関数の定義

\begin{align}e^a=\lim_{k \to \infty} \left(1 + \frac{a}{k}\right)^k\end{align}

から、この漸化式の結果も$\Delta_T \to 0$の場合、次のように表現できます。

\begin{align}
\lim_{\Delta_T \to 0} v_n &= \lim_{\Delta_T \to 0} v_0\left(1 - \frac{c\Delta_T}{m}\right)^n \\
&= v_0\lim_{\Delta_T^{-1} \to \infty} \left(1 - \frac{c}{m\Delta_T^{-1}}\right)^{t\Delta_T^{-1}} \\
&= v_0 (e^{-c/m})^t \\
&= v_0 e^{-ct/m}
\end{align}

これで、$\Delta_T$が十分に小さい場合に結果が一致することがわかりますね。

確認するために、連続時間での結果と各$\Delta_T$値での離散時間の結果を描いて比較してみましょう。

v_0 = 30
m = 1000
c = 200

plt.figure(figsize=[6,5],dpi=100,facecolor='#f0e6f8')

for Delta_T in [0.2,1,2,4,8]:
    t = np.arange(0,16+Delta_T,Delta_T)
    n = np.arange(len(t))
    v = v_0*(1-c*Delta_T/m)**n
    plt.plot(t,v,'o-',ms=2,lw=0.5,label=r'$\Delta_T=%.1f$'%Delta_T)

t = np.linspace(0,16,801)
v = v_0*np.exp(-c*t/m)
plt.plot(t,v,label='$v_0e^{-ct/m}$')

plt.ylabel(r'$v(t)$',size=14)
plt.xlabel(r'$t$',size=14)
plt.legend(fontsize=14)
plt.grid(ls=':')
plt.show()

PythonのlcapyとsympyでのZ変換の実装

最後に、PythonでのZ変換の実装も紹介しておきます。ラプラス変換とは異なり、Z変換はsympyに含まれていませんが、lcapyという線形回路を解析するためのモジュールがあります。この中にはsympyをベースにしたZ変換と逆Z変換があるため、これを使う方法について少し説明します。

インストールはpipで簡単に行えます。

pip install lcapy

Z変換機能を使うにはlcapyモジュールからimportする必要がありますが、使い方はsympyのラプラス変換とほぼ同じなので、詳細は省略します。

from lcapy.ztransform import ZT
import sympy as sp

a,b,x,n,z = sp.symbols('a,b,x,n,z')

x = a**n + b*n**2
X = ZT(x,n,z)
print(X) # b*z*(z + 1)/(z - 1)**3 + z/(-a + z)

latexで表示することもできます。

print(sp.latex(X))

\frac{b z \left(z + 1\right)}{\left(z - 1\right)^{3}} + \frac{z}{- a + z}

逆Z変換の例。

from lcapy.inverse_ztransform import IZT

inv = IZT(X,z,n)
print(inv) # Piecewise((a**n + b*n**2, n >= 0))
print(sp.latex(inv))

latexの結果。

\begin{cases} a^{n} + b n^{2} & \text{for}\: n \geq 0 \end{cases}

ラプラス変換と同様に、0以上の領域で定義されているため、結果にも条件が付いていますね。

その他にも、controlモジュールはZ変換の伝達関数を扱うことができます。これについては今後紹介する予定です。

ラプラス変換とZ変換の対比

以上の説明で、様々な側面からラプラス変換とZ変換を比較してきましたので、ここにまとめておきます。

ラプラス変換	Z変換
$\begin{align}\mathcal{L}[x(t)] = \int_{0}^{\infty}x(t)e^{-st}dt\end{align}$	$\begin{align}\mathcal{Z}[x_n] = \sum_{n=0}^{\infty} x_n z^{-n}\end{align}$
$X(s) = \mathcal{L}[x(t)]$	$X(z) = \mathcal{Z}[x_n]$
連続	離散

微分方程式の解析	漸化式の解析
$t$	$n\Delta_T$
$t \in [0,\infty)$	$n \in \{0, 1, 2, \cdots\}$
$e^{s\Delta_T}$	$z$
$s$	$\begin{align}\frac{\ln(z)}{\Delta_T}\end{align}$
$\begin{align}\mathcal{L}[1] = \frac{1}{s}\end{align}$	$\begin{align}\mathcal{Z}[1] = \frac{1}{1-z^{-1}}\end{align}$
$\begin{align}\mathcal{L}[e^{-\alpha t}] = \frac{1}{s+\alpha}\end{align}$	$\begin{align}\mathcal{Z}[a^n] = \frac{1}{1-az^{-1}}\end{align}$
微分	前の値
$\mathcal{L}[\dot{x}(t)] = sX(s) - x(0) $	$\mathcal{Z}[x_{n-1}] = z^{-1}X(z) + x_{-1} $
積分	次の値
$\begin{align}\mathcal{L}\left[\int_{0}^{t}x(\tau)d\tau\right] = \frac{1}{s} X(s)\end{align} $	$\mathcal{Z}[x_{n+1}] = zX(z) - x_0z $

纏め

この記事では、連続時間システムを記述するラプラス変換に対応する、離散時間システムのためのZ変換について解説しました。

Z変換は、コンピューターで扱う離散化された時間データ（時系列）に対する変換であり、「離散化されたラプラス変換」と捉えることができます。
連続時間の微分方程式が離散化されると漸化式（差分方程式） に変換され、Z変換はこの漸化式を解くための強力なツールとなります。
ラプラス変換とZ変換は、線形性や、微分/積分に対応するシフト性など、多くの類似した性質を持っています。
実際の物理問題（例：空気抵抗を受ける物体の運動）を離散化し、Z変換を用いて解く手順を示しました。得られた離散解は、離散化間隔ΔTを小さくするほど、ラプラス変換で得られた連続解に収束することが確認できました。
Pythonのlcapyモジュールを用いることで、Z変換と逆Z変換をシンボリックに実行できることも紹介しました。

Z変換は、デジタル制御やディジタル信号処理など、現代の工学分野において重要な数学的基礎となっています。

参考

動画