1. 線形シフト不変システムと z 変換 ここで言う システム とは? 入力数列 T[ ] 出力数列 一意変換 ( 演算子 ) 概念的には,, x 2, x 1, x, x1, x2, を入力すると, y 2, y 1, y, y1, y2, が出力される. 線形システム : 線形システムの例 x nxn 1 yn= 2 線形でないシステムの例 xn yn={ 2 xn xn othewise なぜ線形システム? 簡単だから いろいろな性質がわかっている 設計しやすさ実際に使われる非線形システム 整流 ニューラルネット シフト不変システム yn=t [ xn] のとき yn k=t [ x n k] ここで言う x n, y n は {, x 1, x, x1,}, {, y 1, y, y1,} のことをさすことに注意シフト不変システムの例 yn=t [ xn]= xn xn 1 x ' n=x n 1, y' n=t [ x ' n] とすると, y' n=t [ x ' n]= x' n x ' n 1=xn 1 xn 2 yn 1=x n 1 x n 2 より T [ xn 1]= y n 1 シフト不変でないシステムの例 yn=t [ xn]=x n x x ' n=x n 1, y' n=t [ x ' n] とすると, y' n=t [ x ' n]= x' n x' =xn 1 x 1 だが, yn 1=xn 1 x より T [ xn 1] yn 1
なぜシフト不変システム? 簡単だから 入力の時間ずれをあまり気にしなくて良い 線形でシフト不変なシステムの場合, 次の式でシステムが表現できる yn=t [ xn]= hk x n k= hn kx k k= k= h: インパルス応答単位インパルス n : n= で 1, それ以外で T [ n]= hk xn k =hn k= 因果性 : yn は, xn 1, xn のみに依存 n が時間を表す時, 実世界で起きる現象に対する制約 時系列の処理には重要 画像などの場合はあまり重要でない (n は時間ではないので ) このとき, システムは次の形 T [ xn]= h k xn k 安定性 : x n が有界の時, yn も有界 k = hk システムが安定 証明 :h(k) の無限和が有限ならばシステムが安定であることの証明 x n が有界なので xn M となる M が存在する. したがって hk xn k hk xn k k= k= hk M =M hk k= k= システムが安定ならば h(k) の無限和が有限であることの証明 x n が有界の時, yn も有界である. ここで h n x n= h n, h n= h n とおけば, 明らかに xn は有界である. ここで, h n hn hn のとき, hn x n= = hn hn hn= のとき, hn x n== hn となる. ここで, より, である. y = k= k = hk hk x k = k= h k
安定でないシステムの例 : yn= xn yn 1 入力 xn=un ( 単位ステップ ) のとき, yn=n1 であり, yn は発散する. z 変換 : 離散系列でのラプラス変換のようなもの形式的には次のとおり (z は複素数 ) X z=z [ xn]= xn z n n= ( 本によっては z n を掛ける定義もある ) 無限級数なので収束しないかもしれない : 収束性が問題 例 1: x n=a n un X z= n= z=a に極を持つ a n un z n = az 1 n = 1 n= 1 az = z 1 z a a z で収束 az 1 1 a 収束領域は極を含まない 因果性信号の収束領域は無限遠点を含む 例 2: x n= a n u n 1 X z= n= z=a に極を持つ a n z n = a 1 z n=1 1 a 1 z = z z a z a で収束 a n u n 1 z n = a 1 z 1 a 反因果性信号の収束領域は原点を含む z 変換の関数の形と収束領域を合わせないと元の x(n) が特定できない 例 3: x n=a n u n b n u n 1 X z= z z a z a z b z b
a b 収束領域と因果性因果性を z 領域で見ると? 原点から最も遠い極の外側が収束領域因果的な信号 システムだけを問題にする場合はこのタイプだけ考えればよい a 収束領域と安定性安定性を z 領域で見ると? 収束領域が単位円を含む 単位円は z=e j すなわち周波数領域収束領域が単位円を含む = 離散フーリエ変換が可能 安定でないシステムの周波数応答例 yn= xn yn 1 z 領域では Y z= 1 1 z X z 1 極は z=1 なので収束領域は単位円を含まない z=e j とおくと, H e j 1 = 21 cos
スペクトルを見ればわかるように このシステムは直流信号に対して発散するそれ以外では安定 演習 : xn= 1 n un の場合, yn はどうなる? X z= n= z n = 1 1z 1 Y z= 1 1 z 1 1 1z 1 = 1 21 z 1 1 21 z 1 yn= 1 1n un 2 極が単位円の外側にある例 Y z= 1 X z 1 1 2z 時間領域では yn=2y n 1 xn 過去に向かえば収束, 未来に向かえば発散 安定で因果的 すべての極が単位円の内側にある 演習 : 次のシステムを因果的と仮定したとき, 安定かどうかを調べよ. (a) yn=.5y n 1.5y n 2 xn z 変換すると Y z= z 1 z 2 Y z Y z X z 2 2 2 z 1 z 2 Y z= X z 2 2 Y z= X z= 2z2 2 z 1 2 z 2z 2 z1 X z 1± j 7 極はその絶対値は 1 4 2 より, 安定.
(b) yn= 3 2 yn 2 3 2 yn 1 1 2 x n z 変換すると Y z= 3 2 z 2 Y z 3 2 z 1 Y z 1 2 X z 23z 1 3z 2 Y z= X z 1 Y z= 23z 1 3z X z= z 2 2 2z 2 3z3 X z 3± j 15 極は絶対値は 3 4 2 なので安定ではない. 単位ステップ応答 2 175 15 125 1 75 5 25-25 -5-75 -1-125 -15-2.5 2.5 5 7.5 1 12.5 15 17.5 2 22.5 25 27.5 3 32.5 35 37.5 4 42.5 45 z 2 2z 2 3z3 z z 1 = z z z 1 z z z 2 z z 1
3 j z 1 = 15 4 3 j 15 z 2 = 4 15 j 15 = 8 15 j = 15 8 = 6n / 2 = 1 8 yn= z 1 n u n z 2 n u n 1 8 u n 4 15 sin n15 cos n u n 1 8 u n 15 ただし =tan 1 3 z 変換とラプラス変換 ラプラス変換 X s= x t e st dt ここで yt= xt t k のラプラス変換は, Y s= k = = xk e sk k= e s =z とおけば Y z= k= k= xt t ke st dt xk z k s 平面と z 平面の関係 s= j とおくと, e s =z より z=e e j
2.FIR フィルタと IIR フィルタ FIR フィルタ FIR フィルタの特徴 yn= h k xn k N Y z= 1 hk z k X z 必ず安定 直線位相にできる 極は z= のみ, それ以外は零点 零点だけで周波数特性を作る 直線位相の FIR フィルタ z 3 z 1 1/z 1 z 4 z 2 1/z 2 z* 1 z* 3 1/z* 1 単位円および軸に対して対称 FIR フィルタの構成 インパルス応答からフィルタを構成する yn= h k xn k N 点のインパルス応答からフィルタを構成する 周波数特性からフィルタを構成する
hn の離散フーリエ変換 2 j N W =e hn= 1 N とすると H W k W kn 2 jk H e N がわかっているとする H z= hnz n = 1 n= N n= = 1 H W k W k z 1 n N n= = 1 N 1 H W k 1 W k z 1 N N 1 W k z 1 = 1 z N H W k N 1 W k z 1 H W k W kn z n つまり, 周波数軸上 ( ~2 ) での N 点における周波数特性がわかっていれば, N 点 FIR フィルタのシステム関数が構成できる IIR フィルタ yn= hk xn k hk z X z k Y z= 特徴 無限を有限に : 次の形のシステム関数が用いられることが多い Y z= M 1 k = 1 k =1 a k z k b k z k X z うまく設計しないと不安定 うまく設計すれば少ない次数で鋭い特性 非直線位相 IIR フィルタの設計 インパルス応答は無限に長い インパルス応答からの設計は不可能 アナログフィルタからの設計適当なアナログフィルタ (s 領域 ) ディジタルフィルタ (z 領域 ) 直接設計 z 領域で極配置を決定する アナログフィルタからの設計の例 バタワース低域フィルタ振幅特性
H e j 1 = 1 / c = 1 2N 1 1N j / c 2N ( アナログ ) バタワースフィルタの伝達関数 H s の極 = 1 1 N s/ c 2N N が偶数の時 N が奇数の時 2k1 2N s= c e j s= c e k N j の根 2N 個の極のうち,s 平面の左半分 ( 実数部が負 ) を選ぶ ( 安定性のため ) p, p 1, N N p 伝達関数は H s= k k=1 s p k バタワースフィルタの伝達特性の分母 : バタワース多項式 1: s1 2: s 2 2 s1 3: s 2 s1s1 4: s 4 As 3 Bs 2 As1 A=21sin 3 8 cos 3 8 B=22 ( アナログ ) バタワースフィルタの設計 1. 遮断周波数 c [ad/s] を決める 2. 阻止域端周波数 s [ad/s] を決める 3. 阻止域減衰量 A s [db] を決める 4. 1 log 1 H e j s As となるように次数 N を決める N 1 log 1 1 As/1 1 2 log 1 s / c アナログからディジタルへ単純な方法 : アナログバタワースフィルタ N H a s= k =1 1 N s p k = k=1 A k s p k
そのインパルス応答 インパルス応答をサンプリング A z 変換 H z= k k 1 e p k 1 z h a t = A k e p t k ut k hn= A k e p n k u n k 問題点 : アナログフィルタの特性はナイキスト周波数で切れているわけではない 本当は h a t をサンプリングしてはいけない ( 折り返し歪みが発生 ) 双 1 次 z 変換周波数, でのアナログフィルタの周波数特性 周波数 [, ] に写像してしまう LPF の場合, どうせ端の値は小さいので大して問題にならない変換式 s=2 1 z 1 1z 1 導出 : まず s= j 次に =2 tan 1 2 とおく. =2 tan 2 である. z=e j とすれば, となる. 2 1 z 1 1z =2 1 e j 1 1e =2 e j /2 e j /2 2 jsin / 2 =2 j e j /2 j /2 e 2cos/ 2 =2 j tan 2 = j =s 例 :2 次のバタワース低域フィルタ, 遮断周波数 /2 [ad/s] 1. アナログフィルタでの遮断周波数は 2 tan/ 4=2 2. バタワースフィルタの極は 2 21± j = 21± j 2 3. 伝達関数は H s= 21 j 21 j s21 j s21 j = 4 s 2 22 s4
4. H 2 1 z 1 1z 1 = z1 2 22 z 2 22 5. 最終的に H z= 12z 1 z 2 22 2 2 z 2 = 1 22 12z 1 z 2 1 2 2 22 z 2 x(n).2929 + + y(n) z -1.5858 z -1 z -1.2929 -.1716 z -1
3. 線形予測 システムの出力だけを見てシステムが推定できるか? ( 例 ) 音声信号株価 為替 推定手順 システムの 形 を設定する ( モデル化 ) 出力系列 y(n) を最もよく説明するようにシステムパラメータを推定 yn= f, yn 1en en が入力 xn に相当 ( するが, 本当のところはわからない ) 何が 最も良い のか? 自明な解 : yn=en でもあまり意味はない 最も良い モデル化のための仮定 現在の出力 yn は過去の出力値に依存して決まる en をできるだけ小さくするすなわち E [ yn f, yn 1] を最小化 推定すると何が良いのか? 未来の予測ができる (= 儲かる ) システムがマクロレベルで変わるとき, それを見つけることができる ( 音声のモデル化, 音声認識 ) 音声信号の圧縮 ( 符号化 ), 認識の特徴量など f として線形関数を仮定 線形予測モデル p yn= a k yn k en a k k =1 : 線形予測係数 (LPC) 観測データから推定 p: 予測次数 あらかじめ決めておく 注 : 必ずしも線形である必要はない ( ニューラル予測モデル, カーネル予測モデルなど ) 非線形の場合,e(n) を小さくするのは容易だが, パラメータ推定が難しく, また汎化能力が低い ( 過学習 ) 線形予測で表されるシステムはどんなものか? 上式をz 変換 Y z= a k z k Y z E z k 1 Y z= E z p 1 a k z k k=1 IIR フィルタの特殊なケース ( 分子が 1) 分子が 1 なら全極モデル, 自己回帰モデルまたは AR モデル 分母が 1 なら全零モデル, 移動平均モデルまたは MA モデル
どちらも 1 でなければ極零 (ARMA) モデル ちなみに経済学ではわざと安定でないモデルを使うことがある (ARIMA モデル ) 線形予測係数の推定 e(n) の2 乗和が最小になるように係数を決定する p 2 e n 2 = y n a k yn k min n n k =1 行列表現線形予測係数の p 次元ベクトル A=a 1,,a p T 出力系列の (N-p) 次元ベクトル V = y p, y p1,, y T 出力系列の N p p 次元行列 このとき = y F p 1 y p 2 y y p y p 1 y1 y N 2 y N 3 y N p 1 V FA 2 を最小化 d da V FAT V FA=O d da V T V A T F T V V T FA A T F T F A=O 2 F T V 2 F T F A=O F T F A= F T V A= F T F 1 F T V 実際の解き方 F T F A= F T V を分解してみる まず = y p 1 y p yn 2 y p 1 y p 2 y F T y p 2 y p 1 yn 3 y p y p 1 y1 F y y1 y N p 1y N 2 yn 3 y N p 1 = N p 1 N p 1 y p 1k y p 2k y p 1k yk N p 1 N p 1 y p 2k 2 y p 2k y k k = N p 1 N p 1 yk y p 1k yk 2 N p 1 y p 1k 2 N p 1 y p 2k y p 1k N p 1 y k y p 1k そこで F T F = ij とおけば 1 i, j p, また ij = n= p yn i y n j
より = F T V N p 1 k = N p1 k = F T V = 1 2 p y pk y p 1k y pk y p 2k N p1 y pk yk k = したがって全体ではこうなる ( ij = ji ) 11 12 1p 21 22 2p p1 p2 ppa1 a 2 2 a p= 1 p これを Yule-Walke 方程式 ( 正規方程式 ) という 定義どおり ij = yn i y n j として係数を求めるやり方を共分散法という n= p より高速な近似解法として自己相関法がある まず, N i 1 ij = yn i yn j= yn yn j i n= p n= p i である ここで,y が定常的であり, N p ならば, 上の式は i と j の差だけに依存すると仮定することができる 自己相関法では N m 1 m= yn y nm n= として = i j ij とする
すると解くべき方程式は 1 p 2 p 1 1 p 3 p 2 1 2 1 p 4 p 3 a 2 2 p 2 p 3 1 a1 a p 1 p 2 1 p= p となる A の左の行列は対称 Toeplitz 行列 このとき,Levinson-Dubin のアルゴリズムを使って, 上の方程式は O p 2 で解くことができる ( 通常の行列だと O p 3 ) Levinson-Dubin のアルゴリズム 1 2 p 2 p 1 1 1 p 3 p 2 1 2 1 p 4 p 3 a 2 p 2 p 3 p 4 1 a1 a p= p 2 p 1 p 2 p 3 1 = 対称 Toeplitz 行列の重要な性質 1 2 p 2 p 1 1 1 p 3 p 2 T 2 1 p 4 p 3 p とすると, p 2 p 3 p 4 = 1 p 1 p 2 p 3 = 1 1 p 1 p 1 T 1 T p 1 p T p 1 1 p 1 p 1 1 ( その 1) p=k のときの解を a k =a 1 k,,a k k T とする. p=1 のとき a 1 1 = 1 より a 1 1 = 1 /
1 p=k のとき, T k = 2 ak k 1 k ここで, T k 1 1 k 1 1 ( その 2) = T k b k k 1 k 2 = 1 k 1 2 a k 1 k k 2 k 1 1 k 3 k 2 k 2 1 k 1 1 k 2 1 k 1 1 k b1 k b 2 k b k 1 k b k ただし = 1 k 1 k 1 k= k i a i i =1 となるベクトル b k k =b 1 b k k T を考える このとき, 対称性より, b k の要素を縦に反転したベクトル = f k について, k k 2 k bk =1 1 1 k 3 k 2 b k 1 k T k f k である k b 2 k b 1 k=1 のとき, b 1 1 = 1 である k>1 について, k 2 k 1 1 k 2 T k 1 k 1 k 1 T k 1 f 1 k 1 1 したがって b k k =b 1 1 k 1 = 1 ただし ただし k f k 1 k b = i=1 k 1 k f = i=1 i b i k 1 k 1 i 1 b k i
k[ 1 T 1 k k f b k 1 b 1 k b k k 1 b 1 f ( その3) その1, その2より したがって T k ak 1 k b 1 k k f b k b k 1 k f b 1 T k k k = 2 bk a = k ak 1 k k bk f k 1 = ] f k 1 =bk k 1 k これにしたがって,k=1 から p まで順次 b k, a k 1 より 1 2 k k= k 1 k を計算していく