連続系を離散化する際の誤差の話
ラプラス変換で表された連続なシステムをz変換で離散化する時には，$ z=e^{sT}を近似する必要がある。その近似法について考察してみる。（数学的な正確さはここでは追求しない） 以下では，簡単のため，サンプリング周期$ T=1とする。これは即ち，サンプリング周波数$ f=1であり，サンプリング角周波数$ \omega_{s}=2\piである。さらに，サンプリング定理（標本化定理）より，ここで考慮すべき角周波数は，$ -\omega_{s}/2 \leq \omega \leq \omega_{s}/2，すなわち$ -\pi \leq \omega \leq \piである。この範囲外の角周波数の信号は，近似方法に関わらず再現できないので，無視する。 ここで，角周波数$ \omegaの正弦波のみを考慮すると，$ s=j\omegaの関係がある。これを$ z=e^{s}並びにその近似式に代入し，$ -\pi \leq \omega \leq \piの間で$ \omegaを動かすと，$ zの複素平面上で軌跡を描く。その軌跡により，近似を可視化しようと考える。 正確な変換$ z=e^{s} \longrightarrow z=e^{j\omega}
前進オイラー近似 $ z\simeq 1+s \longrightarrow z=1+j\omega 後退オイラー近似 $ z\simeq \frac{1}{1-s} \longrightarrow z=\frac{1}{1-j\omega} パデ近似（双一次変換） $ z\simeq \frac{1+\frac{s}{2}}{1-\frac{s}{2}} \longrightarrow z=\frac{1+j\frac{\omega}{2}}{1-j\frac{\omega}{2}} $ \omega=-\pi, -0.9\pi, \cdots, 0, 0.1\pi \cdots \piについてプロットした結果が下の図である。左上が正確な軌跡，右上が前進オイラー，左下が後退オイラー，右下がパデ近似である。
https://scrapbox.io/files/63c61265ac2245001d342c83.png
$ \omega=0のところは全ての場合で$ z=1+j0の点にあるが，$ \omegaが0から離れるにつれて，近似の軌跡は厳密な軌跡から離れていくことが分かる。しかし，パデ近似の場合は軌跡自体が正確な軌跡と同じ単位円であり，言い換えれば$ |z|が保存されていることが分かり，パデ近似が優れた近似であることが一目瞭然である。（そもそも，パデ近似の式の形は，虚軸を単位円に写像する関数である）
さらに，安定性について議論する。s平面では虚軸を境界とした左半面が安定である一方，z平面では単位円を境界としたその内側が安定である。上の図は，$ s=j\omegaを与えたことから，安定限界である虚軸をz平面に写像した線という見方ができる。この観点では，パデ近似は安定限界を良く再現できている一方，前進オイラーはより不安定な方向に，後退オイラーはより安定な方向にずれていると言える。 なお，正確な変換において，$ -\pi \leq \omega \leq \piの範囲の軌跡はちょうど円1周を描く。これは，サンプリング定理が成り立つことを表しているとも言える。つまり，この範囲を超えると，写像が重なって1対1の写像ではなくなるため，原信号を復元できない。
次に2次近似でやってみる。
前進オイラー近似 $ z\simeq 1+s+\frac{s^{2}}{2} \longrightarrow z=1-\frac{\omega^{2}}{2}+j\omega
後退オイラー近似 $ z\simeq \frac{1}{1-s+\frac{s^{2}}{2}} \longrightarrow z=\frac{1}{1-\frac{\omega^{2}}{2}-j\omega}
パデ近似 $ z\simeq \frac{1+\frac{s}{2}+\frac{s^{2}}{12}}{1-\frac{s}{2}+\frac{s^{2}}{12}} \longrightarrow z=\frac{1-\frac{\omega^{2}}{12}+j\frac{\omega}{2}}{1-\frac{\omega^{2}}{12}-j\frac{\omega}{2}}
それぞれの近似の次数を2次に上げてプロットした結果が以下の図である。いずれの場合も明らかに厳密な場合の単位円により近付いたことが分かる。元々単位円に写像されるパデ近似においても，周波数の高い所でより$ z=-1に近付いている。
https://scrapbox.io/files/63c6127c56cd70001df560d2.png
軌跡を描かせるMATLABプログラムを以下に示す。この状態では，1次近似の式はコメントアウトされており，そのまま実行すると2次近似の式で計算される。 code:matlab
omega = linspace(-pi,pi,21);
tiledlayout(2,2);
nexttile
plot(exp(1i*omega),'-o','MarkerSize',4);
title('Correct trajectory')
xlabel('Re(e^j^\omega)')
ylabel('Im(e^j^\omega)')
nexttile
%plot(1+1i*omega,'-o','MarkerSize',4); %1st
plot(1+1i*omega+1/2*(1i*omega).^2,'-o','MarkerSize',4); % 2nd
title('Forward Euler trajectory')
xlabel('Re(e^j^\omega)')
ylabel('Im(e^j^\omega)')
nexttile
%plot(1./(1-1i*omega),'-o','MarkerSize',4); %1st
plot(1./(1-1i*omega+1/2*(1i*omega).^2),'-o','MarkerSize',4); % 2nd
title('Backward Euler trajectory')
xlabel('Re(e^j^\omega)')
ylabel('Im(e^j^\omega)')
nexttile
%plot((1+1i*omega./2)./(1-1i*omega./2),'-o','MarkerSize',4); %1st
plot((1+1i*omega./2+1/12*(1i*omega).^2)./(1-1i*omega./2+1/12*(1i*omega).^2),'-o','MarkerSize',4); % 2nd
title('Pade trajectory')
xlabel('Re(e^j^\omega)')
ylabel('Im(e^j^\omega)')