趣味人のブログ

Cogito, ergo sum. 我思う故に我あり.

9. ディジタル信号処理と Z 変換

本章では,先ず9.1において,サンプリング定理を復習する.次に,9.2において,ディジタル信号処理におけるサンプラーの論理モデルとサンプル値の表現を,9.3 において Z 変換の導出を示す.そして 9.4 以降において Z 変換の定義,公式と定理,ラプラス変換の像関数からの Z 変換,差分方程式の解法,逆 Z 変換,システム伝達関数と周波数応答を復習する.これらは教科書に十分な説明が記載されているため,本資料では数式の証明や導出等は省略し,留意すべき点を Note に示す.最後に 9.10 において,以上のディジタル信号処理の教科書的説明に基づいて設計したローパスフィルタの例を述べ,教科書通りの設計では期待した特性のフィルタは得られない事を示す.

9.1. サンプリング定理

サンプリング (標本化) とは連続時間信号を離散時間信号に変換する方法の事である.多くの場合,連続時間信号の瞬時値を周期的に計測して,これを離散時間信号とする.離散時間信号を計測する周期はサンプリング周期,その逆数はサンプリング周波数と呼ばれる.

変換された離散時間信号から,連続時間信号に含まれる f0 までの周波数成分を再現するには,サンプリング周波数 fsfs ≥ 2 f0 でなければならない.サンプリング周波数 fs の 1/2 の周波数をナイキスト周波数と言う.

離散時間信号では,連続時間信号に含まれるナイキスト周波数 fn = fs / 2 以下の周波数成分 f0 と,ナイキスト周波数を超える周波数成分 N fs ± f0, (N = 1, 2, 3, ...) を区別する事はできない.これは,異名現象 (エイリアシング) と呼ばれる.従って,連続時間信号をサンプリングする場合は,サンプラーの前にアナログフィルタ等を設けて,ナイキスト周波数を超える周波数成分を十分に減衰させる,若しくはナイキスト周波数を超える周波数成分によって生じるノイズを許容する事が必要となる.

9.2. サンプル値の表現

ディジタル信号処理やサンプル値制御におけるサンプラーの論理モデルを図9‑1 に示す.論理モデルにおけるサンプラーは,連続時間信号の入力と,サンプリング周期 T = 1 / fs で生成される単位インパルス関数 (櫛型関数) を乗算する事によって,連続時間信号の入力の瞬時値をインパルス列として抽出する.

f:id:suzumushi0:20170611122150p:plain
図9-1: サンプラーの論理モデル.

即ち,連続時間信号の入力を x (t), (t ≥ 0),インパルス列を p (t), (t ≥ 0) とすると,p (t) は以下の式で記述される.ここで,δ (t) は,5.2 に述べた単位インパルス関数 (デルタ関数) である.

(9-1) pt =xtδt+δt-T+δt-2T+δt-3T+ =xtn=0δt-nT =n=0xnTδt-nT,t0

Note: 単位インパルス関数は理論上,想像上の関数であるため,サンプラーの論理モデルと現実の実装の間には相当な乖離がある事に留意せよ.また,この論理モデルは,連続時間信号の入力 x (t) と周期 T で生成される単位インパルス関数を乗算しているため,連続時間信号の入力 x (t) とインパルス列 p (t) の関係は定数係数線形常微分方程式では記述できない.従って,これらの間には伝達関数や周波数伝達関数は定義されない事にも注意せよ.

9.3. Z 変換の導出

次に,連続時間信号の入力とインパルス列の関係から,Z 変換を導出する.式 (9-1) のラプラス変換を以下に示す.変換には,表5-2 に示した時間軸上の平行移動定理を利用している.

(9-2) Ps=n=0xnTe-nTs

上記の式において z = eTs に置き換えると,以下の通り Z 変換が得られる.

(9-3) Pz=n=0xnTz-n

ここで,連続時間信号の入力 x (t) に対して,x [n] = x (nT) となる数列を離散時間信号とする.x [n] は離散時間信号の n 番目の標本となる.この場合の p (t) 及び P (z) を以下に示す.

(9-4) pt=n=0xnδt-nT,t0 Pz=n=0xnz-n

以上から明らかな通り,Z 変換とはインパルス列のラプラス変換の別表現である.また,インパルス列における各単位インパルス関数の係数が,離散時間信号となる事に注意せよ.

9.4. Z 変換の定義

n ≥ 0 で定義された数列 (離散時間信号) f [n] の Z 変換 F (z) を以下に示す.

(9-5) Fz=n=0fnz-n

例えば f [0] = 3, f [1] = 1, f [2] = 4, f [n] = 0, (n > 2) となる数列 f [n] の Z 変換は以下の通りとなる.

(9-6) Fz=n=0fnz-n=3+z-1+4z-2

ここで z 複素数z = eTs = eT(σ + ) = re jωT と表す事が多い.F (z) が収束する z の範囲を収束域と言う.f [n] と F (z) の関係を以下の様に表示する場合もある.

(9-7) Zfn=Fz

同様に,f [n] は F (z) の逆 Z 逆変換とも呼ばれ,以下の様に表示される.

(9-8) Z-1Fz=fn

逆 Z 変換の公式を以下に示す.ここで,CF (z) の収束領域内にあり,原点の周りを反時計方向に回る閉路である.

(9-9) fn=1j2πCFzzn-1dz

Z 変換は差分方程式の解法やディジタル信号処理の解析に用いられるが,フーリエ解析の様に観測した離散時間信号を Z 変換する事は稀である.このため,ラプラス変換と同様に,上記の逆 Z 変換の公式を実際に解く必要は殆ど無く,多くの場合逆変換には以下の Z 変換の公式や定理が用いられる.

9.5. Z 変換の公式と定理

Z 変換の主要な公式を以下に示す.

表9-1: Z変換の主要な公式.
項番 fn,n0 Fz 収束域
1 δn=1,n=00,n0 1  
2 un=1,n00,n<0 11-z-1 z>1
3 1 11-z-1 z>1
4 n z-11-z-12 z>1
5 αn 11-αz-1 z>|α|
6 nαn αz-11-αz-12 z>|α|

項番 1 は単位インパルス数列と呼ばれる.あの難解な単位インパルス関数が離散時間信号ではこんなに簡単になると喜んでいる者は,9.2 からの説明を再度読む事.離散時間信号とはインパルス列における各単位インパルス関数の係数である.項番 1 の数列を係数とするインパルス列は,単位インパルス関数となるため,これを単位インパルス数列と称しているだけの事である.

項番 2 は単位ステップ数列と呼ばれる.f [n] の定義域は n ≥ 0 であるから,これを項番 3 の様に表す場合もある.また,項番 5, 6 の α には複素数が含まれる.尚,項番 2 以降は,何れも高校数学で学習した以下の等比数列の和の公式から求める事ができる.

(9-10) n=0arn=a1-r,r<1

最後に,Z 変換の主要な定理を以下に示す.

表9-2: Z変換の主要な定理.
線 形 性 Zafn+bgn=aZfn+bZgn=aFz+bGz
畳み込み Zk=0fn-kgk=ZfnZgn=FzGz
時間軸のシフト Zfn-k=z-kZfn=z-kFz,k0
指数数列の乗算 Zanfn=Fa-1z
F (z) の微分 Znfn=-zddzFz
終値定理 f=limz11-z-1Fz
初期値定理 f0=limzFz

9.6. ラプラス変換の像関数からの Z 変換

ラプラス変換の原関数とその像関数,及び原関数を周期 T でサンプリングした数列とその Z 変換の対応関係を表9‑3 に示す.この表から,時間関数を周期 T でサンプリングした数列の Z 変換を,その時間関数のラプラス変換から求める事ができる.

表9 3: ラプラス変換と Z 変換の対応関係.
原関数 ft , t 0 像関数 Fs fn,n0 Fz
ut=1,t00,t<0 1s un=1,n00,n<0 11-z-1
t 1s2 nT Tz-11-z-12
e-at 1s+a e-naT 11-e-aTz-1

例として,連続時間関数 eat cos (ωt) を周期 T でサンプリングした数列の Z 変換を,以下のラプラス変換の像関数から求める.

(9-11) Le-atcosωt=s+as+a2+ω2

ここで,右辺は以下の通り部分分数に分解できる.

(9-12) s+as+a2+ω2=121s+a+jω+1s+a-jω

表9-3 より,上記の像関数に対応する Z 変換を以下に示す.

(9-13) Zs+as+a2+ω2 =1211-e-a+jωTz-1+11-e-a-jωTz-1 =122-e-a+jωT+e-a-jωTz-11-e-a+jωTz-11-e-a-jωTz-1 =122-e-aTe-jωT+ejωTz-11-e-aTe-jωT+ejωTz-1+e-2aTz-2 =1-e-aTcosωTz-11-2e-aTcosωTz-1+e-2aTz-2

別の例として,連続時間関数 eat sin (ωt) のラプラス変換を以下に示す.

(9-14) Le-atsinωt=ωs+a2+ω2

上記の例と同様に,これを周期 T でサンプリングした数列の Z 変換は以下の通りとなる.

(9-15) Zωs+a2+ω2=e-aTsinωTz-11-2e-aTcosωTz-1+e-2aTz-2

Note: 原関数が単位インパルス関数となる場合は,これをサンプリングしても離散時間信号は無限大に発散してしまうため,Z 変換の単位インパルス数列には対応しない事に注意せよ.

9.7. Z 変換による差分方程式の解法

Z 変換による下記の 1 階の差分方程式 (定数係数線形差分方程式) の解法を示す.高階の場合も解法は原理的に同じである.

(9-16) yn-ayn-1=xn

上記の両辺の Z 変換を以下に示す.

(9-17) Yz-az-1Yz=Xz

ここで,X (z), Y (z) は各々数列 x [n], y [n] の Z 変換を表す.上記を Y (z) について解くと以下の通りとなる.

(9-18) Yz=Xz1-az-1

上記は数列 y [n] の Z 変換を与えるから,これを下記の様に逆 Z 変換する事によって y [n] が得られる.

(9-19) yn=Z-1Xz1-az-1

Note: 式 (9‑16) において,y [0] を求める際に y [−1] が必要となるが,ディジタル信号処理では,n < 0 における離散時間信号 f [n] は 0,即ち初期条件を 0 とする事を暗黙の前提としている.これは,表9‑2 の時間軸シフトの定理は f [n] = 0, (n < 0) を前提としている事からも明らかである.ディジタル信号処理ではまず必要とされないが,差分方程式の初期値問題を解く場合は,表9‑2 の時間軸シフトの定理を k < 0 とした場合が必要となる.参考までに k = −2 の場合の時間軸シフトの定理を以下に示す.

(9-20) Zfn+2=z2Zfn-z2f0-zf1=z2Fz-z2f0-zf1

9.8. 逆 Z 変換

逆 Z 変換の具体的な方法を説明する.説明を簡単にするために,Y (z) がべき級数の比で与えられる下記の関数の逆変換方法を示す.

(9-21) Yz=b0+b1z-1a0+a1z-1+a2z-2

べき級数展開による方法

Y (z) がべき級数の比で与えられる場合は,以下の様に,分子を分母で除算するべき級数展開によって,逆 Z 変換を求める事ができる.

(9-22) Yz =b0+b1z-1a0+a1z-1+a2z-2 =y0+y1z-1+y2z-2+

部分分数分解による方法

ラプラス変換と同様に,逆 Z 変換においても,式 (9‑21) を部分分数に分解して,Z 変換の公式に当てはめる事によって簡易に逆変換ができる.Y (z) の極が単極,若しくは複素数の場合,これらを各々 α, β とすると,式 (9‑21) は以下の様に部分分数に分解できる.

(9-23) Yz=k0+k11-αz-1+k21-βz-1

ここで,k1, k2 は以下の様に留数定理 (ヘビサイドの展開定理) から,k0 は係数比較から求められる.

(9-24) k1=limzα1-αz-1Yz=1-αz-1Yzz=α k2=limzβ1-βz-1Yz=1-βz-1Yzz=β k0=b0a0-k1+k2

以上と Z 変換表より,y [n] は以下の通りとなる.

(9-25) yn=k0δn+k1αn+k2βn

尚,Y (z) の極が n 重極の場合,これを α とすると,Y (z) は以下の様に部分分数に分解できる.

(9-26) Yz=k11-αz-1n+k21-αz-1n-1++kn1-αz-1

ここで,kr は以下の留数定理から求められる.

(9-27) kr =1r-1!limzαdr-1dzr-11-αz-1nYz =1r-1!dr-1dzr-11-αz-1nYzz=α

9.9. システム伝達関数と周波数応答

9.3 に述べた通り,Z 変換とはインパルス列のラプラス変換の別表現であり,かつ離散時間信号の初期条件は 0 であるから,システムに入出力されるインパルス列のラプラス変換の比,即ち Z 変換の比から伝達関数や周波数伝達関数が定義できる.システムの入力の離散時間信号 x [n] と出力の離散時間信号 y [n] が以下の差分方程式で記述されるシステムにおける,入力 x [n] と出力 y [n] の Z 変換の比を,システム伝達関数と言う.

(9-28) yn=k=1Nakyn-k+r=0Mbrxn-r

上記の差分方程式は,出力の離散時間信号 y [n] が,過去 N 回の出力と,現在から過去 M 回の入力の一次結合となるシステムを示している.上記の Z 変換を以下に示す.

(9-29) Yz=k=1NakYzz-k+r=0MbrXzz-r

ここで,X (z), Y(z) は各々入出力の離散時間信号 x [n], y [n] の Z 変換を表す.よってシステム伝達関数 H (z) は以下で与えられる.

(9-30) Hz=YzXz=r=0Mbrz-r1-k=1Nakz-k

また,Z 変換はラプラス変換において z = eTs に置き換えたものであるから,周波数伝達関数は,以下に示す様に,システム伝達関数において z = e jωT に置き換える事によって得られる.

(9-31) HejωT=r=0Mbre-jrωT1-k=1Nake-jkωT

9.10. ローパスフィルタの例

以上,教科書的なディジタル信号処理の基本を復習した.ここで,ディジタルフィルタの設計の各論に移る前に,以下の伝達関数で示される,時定数が τ,カットオフ周波数が ωc = 1 / τ のローパスフィルタをディジタル化した例を示す.

(9-32) Gs=1τs+1=ωcs+ωc

表9‑3 より,上記の伝達関数に対応するシステム伝達関数を以下に示す.サンプリング周期は T とする.

(9-33) Hz=ωc1-e-ωcTz-1

また,上記のシステム伝達関数に対応する,周波数伝達関数を以下に示す.

(9-34) HejωT=ωc1-e-ωcTe-jωT

上記の周波数伝達関数のゲイン線図を図9‑2 に示す.横軸の角周波数は ω / ωc として正規化されている.また図中の青線は正規化されたナイキスト周波数 fn を示しており,この例ではカットオフ周波数 ωc の 20 倍の角周波数としている.尚,サンプリング周波数は fs = 1 / T = 1 [MHz] としている.

f:id:suzumushi0:20170611122151p:plain
図9-2: ディジタル化されたローパスフィルタのゲイン線図の例.

上記のゲイン線図と図6-8 に示した一次遅れ要素 (ローパスフィルタ) のゲイン線図を比較すると,ナイキスト周波数以下の特性は一見似ている様に見えるが,ゲインの最大値が 120 dB 以上となっており特性が全く異なるローパスフィルタとなっている.

実は,ディジタル信号処理の論理モデルは,現実の実装モデルと相当な乖離があるため,理論モデルを解説した教科書通りに設計しても,期待した特性のディジタルフィルタは得られない.そこで,10章において,論理モデルと実装モデルの関係を述べ,論理モデルを実装モデルに適用する方法を明らかにする.

10. ディジタル信号処理の論理モデルと実装モデル - 趣味人のブログ

本章の参考文献

明石一, "制御工学 増訂版," 大学講座 機械工学12, 共立出版, 1979.
A. Oppenheim, R. Schafer, 伊達訳, "ディジタル信号処理," (上), (下), コロナ社, 1978.