1. 線形シフト不変システムと z 変換 ここで言う システム とは? 入力数列 T[ ] 出力数列 一意変換 ( 演算子 ) 概念的には,, x 2, x 1, x 0, x 1, x 2, を入力すると, y 2, y 1, y 0, y 1, y 2, が出力される. 線形システム : 線形シ

１．線形シフト不変システムとｚ変換

ここで言う「システム｣とは？ T[ ] 入力数列 出力数列 一意変換（演算子） 概念的には，  , x −2 , x −1 , x 0 , x 1 , x  2 , を入力すると , y −2 , y −1 , y 0 , y 1 , y  2 , が出力される． 線形システム： 線形システムの例 y n=x  nx n−1 2 線形でないシステムの例 y n= x n2 y n=

{

x n x n0 0 otherwise なぜ線形システム？ ● 簡単だから ● いろいろな性質がわかっている ● 設計しやすさ 実際に使われる非線形システム ● 整流 ● ニューラルネット シフト不変システム y n=T [ x n] のとき yn−k=T [ x n−k] ※ここで言う x  n , y n は {, x −1 , x 0 , x 1 , }, { , y −1 , y 0 , y 1 ,} のことをさすことに注意 シフト不変システムの例 y n=T [ x n]= x n− x n−1 x ' n=x n−1 , y ' n=T [ x ' n] とすると， y ' n=T [ x ' n]= x ' n−x ' n−1=x n−1−x  n−2 y n−1=x  n−1−x  n−2 より T [ xn−1]= y n−1 シフト不変でないシステムの例 y n=T [ x n]=x  n−x 0 x ' n=x n−1 , y ' n=T [ x ' n] とすると， y ' n=T [ x ' n]= x '  n−x '  0=x n−1− x −1 だが， y n−1=x n−1−x 0 より T [ xn−1]≠ yn−1

なぜシフト不変システム？ ● 簡単だから ● 入力の時間ずれをあまり気にしなくて良い 線形でシフト不変なシステムの場合，次の式でシステムが表現できる y n=T [ x n]=

∑

k=−∞ ∞ h k  x  n−k =

∑

k=−∞ ∞ hn−k  x k  h:インパルス応答 単位インパルス  n : n=0 で 1，それ以外で 0 T [ n]=

∑

k=−∞ ∞ h k  x n−k =h n 因果性： y n は  , x n−1 , x n のみに依存 ● n が時間を表す時，実世界で起きる現象に対する制約 ● 時系列の処理には重要 ● 画像などの場合はあまり重要でない（n は時間ではないので） このとき，システムは次の形 T [ x n ]=

∑

k=0 ∞ h k  x n−k  安定性： x  n が有界の時， y n も有界

∑

k =−∞ ∞ ∣_{h k ∣∞ ⇔システムが安定} 証明：h(k)の無限和が有限ならばシステムが安定であることの証明 x  n が有界なので ∣x n∣M となる M が存在する．したがって

∣

∑

k=−∞ ∞ h k  x  n−k 

∣

≤

∑

k=−∞ ∞ ∣h k  x  n−k ∣ 

∑

k=−∞ ∞ ∣h k  M∣=M

∑

k=−∞ ∞ ∣h k ∣∞ システムが安定ならば h(k)の無限和が有限であることの証明 x  n が有界の時， y n も有界である．ここで x  n=∣h −n∣ h −n h n≠0 , 0 h n=0 とおけば，明らかに x  n は有界である．ここで， h n≠0 のとき， h n x−n=h n⋅∣h n∣_hn =∣h n∣ h n=0 のとき， h n x−n=0=∣hn∣ となる．ここで， ∣y 0∣=

∣

∑

k=−∞ ∞ hk  x −k 

∣

=

∣

∑

k=−∞ ∞ ∣h k ∣

∣

∞ より，

∑

k =−∞ ∞ ∣h k ∣∞ である．

安定でないシステムの例： y n= x n y n−1 入力 x  n=u n （単位ステップ）のとき， y n=n1 であり， yn は発散する． ｚ変換：離散系列でのラプラス変換のようなもの 形式的には次のとおり（z は複素数） X  z =Z [ x  n]=

∑

n=−∞ ∞ x n z−n （本によっては zn を掛ける定義もある） 無限級数なので収束しないかもしれない：収束性が問題 例１： x  n=an_{u n} X  z =

∑

n=−∞ ∞ anu n z−n =

∑

n=0 ∞ az−1n= 1 1−az−1= z z −a 

∣

az −1

_∣

1 z =a に極を持つ ∣a∣∣z∣ で収束 例２： x  n=−anu−n−1 X  z =

_∑

n=−∞ ∞ −an_{u −n−1 z}−n =

∑

n=1 ∞ −a−n_zn = −a −1_z 1−a−1z= z z−a 

∣

a −1_z

_∣

_1 z =a に極を持つ ∣z∣∣a∣ で収束 ｚ変換の関数の形と収束領域を合わせないと元の x(n)が特定できない 例３： x  n=an_{u n−b}n_{u −n−1} X  z = z z −a z z−b ∣a∣∣z∣∣b∣ a a 収束領域は 極を含まない 因果性信号の 収束領域は 無限遠点を含む 反因果性信号の 収束領域は 原点を含む

収束領域と因果性 因果性をｚ領域で見ると？ ⇒原点から最も遠い極の外側が収束領域 因果的な信号・システムだけを問題にする場合はこのタイプだけ考えればよい 収束領域と安定性 安定性をｚ領域で見ると？ ⇒収束領域が単位円を含む ・単位円は z =ej  すなわち周波数領域 収束領域が単位円を含む＝離散フーリエ変換が可能 安定でないシステムの周波数応答例 y n= x n y n−1 ｚ領域では Y  z = 1 1− z−1 X  z 極は z=1 なので収束領域は単位円を含まない z=ej  とおくと，

∣

H  ej 

∣

=

_

_{21−cos}1 a b a

スペクトルを見ればわかるように、このシステムは直流信号に対して発散する それ以外では安定 演習： x  n=−1n_{u n の場合， y n はどうなる？} X  z =

∑

n=0 ∞ −z −n = 1 1z−1 Y  z = 1 1− z−1 1 1z−1= 1 21−z−1 1 21 z−1 y n=1−1 n_{u n} 2 極が単位円の外側にある例 Y  z = 1 1−2z−1 X  z  時間領域では y n=2y n−1 x n 過去に向かえば収束，未来に向かえば発散 安定で因果的⇒すべての極が単位円の内側にある 演習：次のシステムを因果的と仮定したとき，安定かどうかを調べよ． (a) y n=0.5y n−1−0.5y n−2 x n

z 変換すると Y  z =z −1 2 Y  z − z−2 2 Y  z  X  z  2− z−1z−2 2 Y  z= X  z Y  z = 2 2−z−1z−2 X  z = 2z2 2z2−z1 X  z  極は 1± j



7 4 その絶対値は 1



2 より，安定．

(b) y n=−3₂ y  n−2−3₂ y n−11₂x n z 変換すると Y  z =−3 2z −2_{Y  z −}3 2z −1_{Y  z }1 2X  z  23z−1_3z−2_{Y  z = X  z} Y  z = 1 23z−13z−2 X  z = z2 2z23z3 X  z  極は −3± j



15 4 絶対値は



3 2 なので安定ではない． 単位ステップ応答 -2.5 0 2.5 5 7.5 10 12.5 15 17.5 20 22.5 25 27.5 30 32.5 35 37.5 40 42.5 45 -1500 -1250 -1000 -750 -500 -250 0 250 500 750 1000 1250 1500 1750 2000 z2 2z23z3 z z −1 = z z−z₁ z z −z₂ z z−1

z₁=−3 j



15 4 z2= −3− j



15 4 =15 j



15 80 =15− j



15 80 =1 8 y n= z₁nu  n z₂nu n1 8u n =6 n / 2 40



15 sin n15 cos n 



u n 1 8u n ただし =tan−1





15 3 −



ｚ変換とラプラス変換 ラプラス変換 X s=

∫

−∞ ∞ x t e−st dt ここで y t=

∑

k=−∞ ∞ x t t−k  のラプラス変換は， Y s =

_∫

−∞ ∞

∑

k =−∞ ∞ x t t−k  e−st_dt =

∑

k=−∞ ∞ x  k  e−sk es_{=z とおけば} Y  z =

∑

k=−∞ ∞ x k  z−k ｓ平面とｚ平面の関係 s= j  とおくと， es =z より z =e⋅ej 

２．FIR フィルタと IIR フィルタ

FIR フィルタ y n=

_∑

k=0 N −1 h k  x n−k  Y  z =



∑

k=0 N −1 h k  z−k



X  z  FIR フィルタの特徴 ● 必ず安定 ● 直線位相にできる ● 極は z=0 のみ，それ以外は零点　→　零点だけで周波数特性を作る

直線位相のFIRフィルタ

z

₄

z

₄

z

₃

z*

₃

z

₃

z*

₃

z

₂

1/z

₂

z

₂

1/z

₂

z

₁

z*

₁

1/z

₁

1/z*

₁

z

₁

z*

₁

1/z

₁

1/z*

₁

単位円および軸に対して対称

FIR フィルタの構成 インパルス応答からフィルタを構成する y n=

∑

k=0 N −1 h k  x  n−k  N 点のインパルス応答からフィルタを構成する 周波数特性からフィルタを構成する

h n の離散フーリエ変換 H



e 2 jk N



_{がわかっているとする} W =e −2  j N _とすると h n= 1 N

∑

k=0 N −1 H



W−k

_

W−kn H  z =

_∑

n=0 N −1 h n z−n =1 N

∑

n=0 N −1

∑

k=0 N −1 H W−k W−kn_z−n = 1 N

∑

k=0 N −1 H W−k 

∑

n=0 N −1 W−k z−1n = 1 N

∑

k=0 N −1 H W−k 1−W −k z−1N 1−W−k z−1 =1−z −N N

∑

k=0 N −1 H W−k  1−W−k z−1 つまり，周波数軸上( 0~2  )での N 点における周波数特性がわかっていれば， N 点 FIR フィルタのシステム関数が構成できる IIR フィルタ y n=

∑

k=0 ∞ h k  x n−k  Y  z =



_∑

k=0 ∞ h k  z−k



X  z  無限を有限に：次の形のシステム関数が用いられることが多い Y  z =

∑

k =0 M −1 a_kz−k 1−

∑

k =1 N −1 b_kz−k X  z  特徴 ● うまく設計しないと不安定 ● うまく設計すれば少ない次数で鋭い特性 ● 非直線位相 IIR フィルタの設計 ● インパルス応答は無限に長い⇒インパルス応答からの設計は不可能 ● アナログフィルタからの設計 適当なアナログフィルタ（ｓ領域）→ディジタルフィルタ（ｚ領域） ● 直接設計 ｚ領域で極配置を決定する アナログフィルタからの設計の例 バタワース低域フィルタ 振幅特性

∣

H  ej 

∣

= 1



1

_

/_c

_

2N = 1



1−1N



j /_c

_

2N （アナログ）バタワースフィルタの伝達関数 H  s の極＝ 1−1N_s/ c2N の根 N が偶数の時　 s=_ce 2k1 2N j  N が奇数の時　 s=_ce k N j  2N 個の極のうち，ｓ平面の左半分（実数部が負）を選ぶ（安定性のため） ⇒ p1,, pN 伝達関数は H  s=∏ k=1 N _−p k s− p_k バタワースフィルタの伝達特性の分母：バタワース多項式 1: s1 2: s2 



2 s1 3:  s2s1 s1

4: s4As3Bs2As1 A=



21sin3₈ cos3 ₈ B=



22 （アナログ）バタワースフィルタの設計 1. 遮断周波数 c [rad/s]を決める 2. 阻止域端周波数 s [rad/s]を決める 3. 阻止域減衰量 As [dB]を決める 4. −10 log10

∣

H  e j s_

∣

_≥A s となるように次数 N を決める ⇒ N ≥1 2 log1010As/10−1 log₁₀_s/_c アナログからディジタルへ 単純な方法： アナログバタワースフィルタ Has =

∏

k =1 N 1 s− p_k=

∑

k=1 N _A k s− p_k

そのインパルス応答 hat =

∑

k Ake −pktu t インパルス応答をサンプリング h n=

∑

k Ake−pknu n ｚ変換 H  z =∑ k A_k 1−e−pkz−1 問題点： アナログフィルタの特性はナイキスト周波数で切れているわけではない →本当は hat  をサンプリングしてはいけない（折り返し歪みが発生） 双１次ｚ変換 周波数 −∞ , ∞ でのアナログフィルタの周波数特性 →周波数 [−, ] に写像してしまう LPF の場合，どうせ端の値は小さいので大して問題にならない 変換式 s=2⋅1−z−1 1z−1 導出： まず s= j  次に =2 tan−1_{2 とおく．} =2 tan_{2 である．} z=ej  とすれば， −≤≤ となる． 21−z −1 1z−1=2 1−e−j  1e−j =2 ej /2−e−j /2 ej /2 e−j /2=2 2 jsin / 2 2 cos/ 2 =2 j tan  2=j =s 例：２次のバタワース低域フィルタ,遮断周波数 /2 [rad/s] 1. アナログフィルタでの遮断周波数は 2 tan/ 4=2 2. バタワースフィルタの極は 2×−



21± j 2 =−



2 1± j 3. 伝達関数は H  s=



21 j s



2 1 j 



21− j s



2 1− j= 4 s22



2 s4

4. _H



₂1−z−1 1z−1



= z12 



22 z2−



22 5. 最終的に H  z = 12z−1z−2



22−



2−2 z−2= 1



22 12z−1z−2 1−



2−2



22 z −2

x(n)

y(n)

z

-1

+

z

-1 0.2929 0.2929 0.5858

+

z

-1

z

-1 -0.1716

３．線形予測

システムの出力だけを見てシステムが推定できるか？ （例）音声信号 　　　株価・為替 推定手順 ● システムの「形｣を設定する（モデル化） ● 出力系列 y(n) を最もよく説明するようにシステムパラメータを推定 y n= f  , y n−1e n e n が入力 x n に相当（するが，本当のところはわからない） ● 何が「最も良い｣のか？ 自明な解： y n=e n でもあまり意味はない ｢最も良い｣モデル化のための仮定 ● 現在の出力 y n は過去の出力値に依存して決まる ● e n をできるだけ小さくする すなわち E [ y n− f  , y  n−1] を最小化 推定すると何が良いのか？ ● 未来の予測ができる（＝儲かる） ● システムがマクロレベルで変わるとき，それを見つけることができる （音声のモデル化，音声認識） →音声信号の圧縮（符号化），認識の特徴量など f として線形関数を仮定→線形予測モデル y n=−

_∑

k =1 p a_ky n−k e n ak :線形予測係数(LPC)→観測データから推定 p:予測次数→あらかじめ決めておく 注：必ずしも線形である必要はない（ニューラル予測モデル，カーネル予測モデルなど） 非線形の場合，e(n)を小さくするのは容易だが，パラメータ推定が難しく，また汎化能力 が低い（過学習） 線形予測で表されるシステムはどんなものか？ 上式をｚ変換 Y  z =−

_∑

k a_kz−k_{Y  z E  z } Y  z = 1 1

∑

k=1 p a_kz−k E  z IIR フィルタの特殊なケース（分子が１) ※分子が１なら全極モデル，自己回帰モデルまたは AR モデル ※分母が１なら全零モデル，移動平均モデルまたは MA モデル

※どちらも１でなければ極零(ARMA)モデル ※ちなみに経済学ではわざと安定でないモデルを使うことがある(ARIMA モデル） 線形予測係数の推定 e(n)の２乗和が最小になるように係数を決定する

∑

n e  n2=

∑

n

∣

y  n

∑

k =1 p a_ky n−k 

∣

2 min 行列表現 線形予測係数の p 次元ベクトル A=a1,, apT 出力系列の(N-p)次元ベクトル　 V = y  p , y  p1 , , y  N −1T 出力系列の  N − p× p 次元行列 F =



y  p−1 y  p−2 ⋯ y 0 y  p y  p−1 ⋯ y 1 ⋮ ⋯ ⋮ y  N −2 y  N −3 ⋯ y  N − p−1



このとき

∣V FA∣2 を最小化→ _dAd V FATV FA=O d

dAV

T_{V A}T_FT_{V V}T _{FA A}T_FT_{F A=O}

2 FT_{V 2 F}T_{F A=O} FTF A=−FTV A=− FT_{F }−1_FT_V 実際の解き方 FT_{F A=−F}T_{V を分解してみる。まず} FT_{F =}



y  p−1 y  p ⋯ y  N −2 y  p−2 y  p−1 ⋯ y  N −3 ⋮ ⋯ ⋮ y 0 y 1 ⋯ y  N − p−1



y  p−1 y  p−2 ⋯ y 0 y  p y  p−1 ⋯ y 1 ⋮ ⋯ ⋮ y  N −2 y  N −3 ⋯ y  N − p−1



=



∑

k=0 N − p−1 y  p−1k 2

∑

k=0 N − p−1 y  p−1k  y  p−2k  ⋯

∑

k=0 N − p−1 y  p−1k  y k 

∑

k=0 N − p−1 y  p−2k  y  p−1k 

_∑

k =0 N − p −1 y  p−2k 2 ⋯

∑

k=0 N− p−1 y  p−2k  y  k  ⋮ ⋮ ⋱ ⋮

∑

k=0 N − p−1 y  k  y  p−1k 

_∑

k=0 N − p−1 y k  y  p−1k  ⋯

∑

k=0 N − p−1 y k 2



そこで FTF =ij とおけば 1≤i , j≤ p ， _ij=

∑

n= p N −1 y n−i y  n− j また

FTV =



∑

k =0 N − p 1 y  pk  y  p−1k 

∑

k =0 N − p1 y  pk  y  p−2k  ⋮

∑

k =0 N − p1 y  pk  y k 



より FTV =



₀₁ ₀₂ ⋮ _0p



したがって全体ではこうなる( ij=ji )



₁₁ ₁₂ ⋯ _1p ₂₁ ₂₂ ⋯ _2p ⋮ ⋮ ⋱ ⋮ _p1 _p2 ⋯ _pp



a₁ a₂ ⋮ a_p



=−



₀₁ ₀₂ ⋮ _0p



これを Yule-Walker 方程式（正規方程式）という。 定義どおり ij=

∑

n= p N −1 y n−i y  n− j として係数を求めるやり方を共分散法という。 より高速な近似解法として自己相関法がある。まず， _ij=

∑

n= p N −1 y n−i y  n− j=

∑

n= p−i N −i −1 y n y  n j−i である。ここで，y が定常的であり， N ≫ p ならば，上の式は i と j の差だけに依存す ると仮定することができる。 自己相関法では r m=

∑

n=0 N − m−1 y n y  nm として ij=r ∣i− j∣ とする。

すると解くべき方程式は



r  0 r 1 ⋯ r  p−2 r  p−1 r 1 r 0 ⋯ r  p−3 r  p−2 r  2 r 1 ⋯ r  p−4 r  p−3 ⋮ ⋮ ⋱ ⋮ r  p−2 r  p−3 ⋯ r 0 r 1 r  p−1 r  p−2 ⋯ r 1 r 0





a₁ a₂ ⋮ ap



=−



r 1 r 2 ⋮ r  p



となる。A の左の行列は対称 Toeplitz 行列。このとき，Levinson-Durbin のアルゴリズム

を使って，上の方程式は O p2 で解くことができる。（通常の行列だと O p3 ） Levinson-Durbin のアルゴリズム



r₀ r₁ r₂ ⋯ r_{p −2} r_{p −1} r₁ r₀ r₁ ⋯ r_p−3 r_p−2 r₂ r₁ r₀ ⋯ r_{p −4} r_{p −3} ⋮ ⋮ ⋮ ⋱ ⋮ rp−2 rp −3 rp−4 ⋯ r0 r1 r_p−1 r_p−2 r_p−3 ⋯ r₁ r₀





a₁ a₂ ⋮ ap



=−



r₁ r₂ ⋮ rp



※対称 Toeplitz 行列の重要な性質 T_p=



r₀ r₁ r₂ ⋯ r_p−2 r_p−1 r₁ r₀ r₁ ⋯ r_p−3 r_p−2 r₂ r₁ r₀ ⋯ r_p−4 r_p−3 ⋮ ⋮ ⋮ ⋱ ⋮ rp−2 rp−3 rp −4 ⋯ r0 r1 r_p−1 r_p−2 r_p−3 ⋯ r₁ r₀



とすると， T_p=



r₀ r₁ ⋯ r_p−1 r₁ ⋮ T_p−1 r_p−1



=



r_p−1 T_p−1 ⋮ r₁ r_p−1 ⋯ r₁ r₀



（その１） p=k のときの解を ak=a1 k ,  , ak k T とする． p=1 のとき　 r0a1 1_=−r 1 より a1 1_=−r 1/r0

p=k のとき， Tka k =−



r1 r₂ ⋮ r_k−1 rk



ここで，



r_k−1 rk−2 T_{k −1} ⋮ r1 rk −1 ⋯ r1 r0





ak−1 0



=−



r₁ r2 ⋮ r_k−1 _k



ただし k=−

∑

i =1 k −1 r_k−ia_ik−1 （その２） T_kbk  =



r₀ ⋯ r_k−2 r_k−1 r1 ⋯ rk−3 rk−2 ⋮ ⋮ ⋮ r_k−2 ⋯ r₀ r₁ rk−1 ⋯ r1 r0





b1 k  b₂k  ⋮ b_k−1k b_kk 



=



0 0 ⋮ 0 1



となるベクトル bk =b₁k _⋯b k k_T を考える。このとき，対称性より， bk  の 要素を縦に反転したベクトル f k  について， T_k fk =



r₀ ⋯ r_{k −2} r_{k −1} r1 ⋯ rk−3 rk −2 ⋮ ⋮ ⋮ r_{k −2} ⋯ r₀ r₁ rk −1 ⋯ r1 r0





bkk bk−1 k  ⋮ b₂k b₁k



=



1 0 ⋮ 0 0



である。 k=1 のとき， b1 1 =1 r₀ である。 k>1 について，



r₀ ⋯ r_{k −2} r_k−1 r₁ ⋮ r_{k −2} T_k−1 rk−1





0 bk −1



=



_bk  0 ⋮ 0 1



ただし b k  =

∑

i=1 k −1 r_ib_ik−1



r_k−1 T_{k −1} ⋮ r₁ r_{k −1} ⋯ r₁ r₀





fk−1 0



=



1 0 ⋮ 0 _fk 



ただし f k  =

∑

i=1 k −1 r_i−1b_{k −i}k−1 したがって

T_k

[

1 1−_fk _bk 



0 bk −1 



− _bk 1−_fk _bk 



fk−1 0



]

=



0 0 ⋮ 0 1



より 1 1−f k  _bk 



0 bk−1



− bk  1−f k  _bk 



fk−1 0



=b k  （その３） その１，その２より T_k



ak−1 0



−rk−kTkb k  =−



r₁ r₂ ⋮ rk−1 k



−



0 0 ⋮ 0 r_k−_k



=−



r₁ r₂ ⋮ rk −1 rk



したがって ak =



a k−1 0



−rk−kb k  これにしたがって，k=1 から p まで順次 bk , ak  を計算していく。