SPEKTRÁLNÍ ANALÝZA ČASOVÝCH ŘAD prof. Ing. Jiří Holčík, CSc. holcik@iba.muni.cz © Institut biostatistiky a analýz V. PARAMETRICKÉ METODY ODHADU VÝKONOVÉHO SPEKTRA pokračování MA A ARMA MODELY © Institut biostatistiky a analýz ZÁKLADNÍ VZTAHY MEZI PARAMETRY MODELU A AUTOKORELAČNÍ POSLOUPNOSTÍ ARMA: y(nTvz ) = −∑ ak y(nTvz − kTvz ) + ∑ bk x(nTvz − kTvz ) k =1 p k =0 p q .y(nTvz − mTvz ), E E[ y(nTvz )y(nTvz − mTvz )] = −∑ akE[ y(nTvz − kTvz )y(nTvz − mTvz )] + k =1 + ∑ bkE[ x(nTvz − kTvz )y(nTvz − mTvz )] k =0 q γ yy (mTvz ) = −∑ ak γ yy (mTvz − kTvz ) + ∑ bk γ xy (mTvz − kTvz ) k =1 k =0 p q γ xy(mTvz) … vzájemná korelační posloupnost mezi x(nTvz) a y(nTvz) ∞  γ xy (mTvz ) = E[ y(nTvz )x(nTvz + mTvz )] = E ∑ h(kTvz )x(nTvz − kTvz ).x(nTvz + mTvz ) =  k =0  = ∑ h(kTvz )E[ x(nTvz ).x(nTvz + mTvz + kTvz )] = ∑ h(kTvz )γ xx (mTvz + kTvz ) = k =0 k =0 ∞ ∞ předpokládáme kauzální filtr 0, m>0  = h( −mTvz ).σ2 =  2 x σ x .h( −mTvz ) m ≤ 0 © Institut biostatistiky a analýz ZÁKLADNÍ VZTAHY MEZI PARAMETRY MODELU A AUTOKORELAČNÍ POSLOUPNOSTÍ ARMA: p  − ∑ ak .γ yy (mTvz − kTvz )  k =1  p q −m  2 γ yy (mTvz ) = − ∑ ak .γ yy (mTvz − kTvz ) + σ x ∑ h(kTvz ).bk +m k =0  k =1 γ yy ( −mTvz )    m>q 0≤m≤q m<0 nelineární vztah mezi γ yy(mTvz) a parametry ak a bk, protože impulzní charakteristika h(kTvz) = f(ak,bk) © Institut biostatistiky a analýz ZÁKLADNÍ VZTAHY MEZI PARAMETRY MODELU A AUTOKORELAČNÍ POSLOUPNOSTÍ  opakování: AR: y(nTvz ) = −∑ ak y((n − k )Tvz ) + ∑ bk x((n − k )Tvz ) k =1 k =0 p q σ 2 b 0 = 1 a γ xx (kTvz ) =  x 0 p   − ∑ ak .γ yy (mTvz − kTvz )  p k =1  γ yy (mTvz ) = − ∑ ak .γ yy (mTvz − kTvz ) + σ 2 x  k =1 γ yy ( −mTvz )    k=0 k≠0 m > 0, (q = 0) m = 0, (0 ≤ m ≤ q, q = 0) m<0 © Institut biostatistiky a analýz MA MODELY  opakování:  0  q −m  γ yy (mTvz ) = σ 2 ∑ bk bk +m = σ 2 dm x x  k =0 γ yy ( −mTvz )   m>q 0≤m≤q m<0 yn = m =0 ∑b q m x n −m h(k ) = {bk } δ 2 E { x n +m x n } =  x 0 m=0 jindy E{ x n } = 0 © Institut biostatistiky a analýz MA MODELY B( z ).B( z ) = D( z ) = −1 m= − q dm z −m ∑ q (b0z0+b1z1+…+bqzq).(b0z0+b1z-1+…+bqz-q) = = b02+b0b1z+…+b0bqzq+b0b1z-1+b12+b2b1z+…+bqb1zq-1+…+ +b0bqz-q+b1bqz-q+1+b2bqz-q+2+…+bq2 = = (b0b0+b1b1+…+bqbq).z0 +(b0b1+b1b2+…+bq-1bq).z-1 + d0 d-1 +(b1b0+b2b1+…+bqbq-1).z1 +…+(…).zq ⇓ dm = q− m k =0 ∑b b k d1 k +m pro m ≤ q © Institut biostatistiky a analýz MA MODELY Tedy  0 γ yy (mT ) =  2 σ x dm q m >q m ≤q metoda momentů metoda momentů !!! nepočítáme parametry MA modelu, ale jen výkonovou spektrální hustotu !!! a výkonové spektrum MA Γyy ( f ) = m= − q γ yy (mTvz ).e − j2 πfmTvz ∑ odhad spektra: Pyy ( f ) = MA m= − q ryy (mTvz ).e − j2 πfmTvz ∑ aneb aneb metoda Blackmanova- Tukeyho metoda Blackmanova- Tukeyho © Institut biostatistiky a analýz q !!! A TO JE KLASIKA !!! !!! A TO JE KLASIKA !!! © Institut biostatistiky a analýz MA MODELY ALTERNATIVA: stanovení {bk} založené na aproximaci MA procesu AR procesem vysokého řádu, tj. p » q  dekompoziční  teorém (Wold 1938) jakýkoliv ARMA nebo MA proces může být jednoznačně reprezentován AR modelem max. ∞ řádu; jakýkoliv ARMA nebo AR proces lze reprezentovat MA modelem max. ∞ řádu; © Institut biostatistiky a analýz  MA MODELY ALTERNATIVA: stanovení {bk} založené na aproximaci MA procesu AR procesem vysokého řádu, tj. p » q v tom případě platí, že B(z) = 1/A(z), resp. B(z).A(z)=1 a proto q   1 n = 0 an + ∑ bk .an−k =  k =1 0 n ≠ 0 frekvenční charakteristiky Hammingovy MA DP s délkou impuzní odezvy 322 vzorků (modrá) a její AR ekvivalent řádu 274 (černá) © Institut biostatistiky a analýz MA MODELY ALTERNATIVA: p odhad {bk} můžeme zpřesnit pomocí metody nejmenších čtverců určíme kvadratickou chybu      e = ∑ an + ∑ bk .an−k  − 1, a0 = 1; ak = 0 pro k < 0 n =0  k =0  q 2 a tu minimalizujeme výběrem parametrů {bk}  −1 b = −R aa .raa p − i− j n =0 p −i kde   R aa ( i − j ) = ∑ an .an+ i− j , pro i, j = 1, 2, ..., q;   raa (i) = ∑ an .an+i , pro i = 1, 2, ..., q. n =0 vymyslel to Durbin a ukázalo se, že je to přibližně odhad vymyslel to Durbin a ukázalo se, že je to přibližně odhad s maximální věrohodností, je-li proces normální s maximální věrohodností, je-li proces normální © Institut biostatistiky a analýz URČENÍ ŘÁDU MA MODELU  pro vyjádření úzkých pásem je třeba hodně vysoký řád; hodnoty odhadu autokorelační funkce musí rychle klesat k nule, protože γ yy(mTvz)=0 pro | m| >q test, zda ryy(mTvz)→0 je založený na srovnávání ryy(qTvz) s rozptylem hodnot ryy(mTvz) pro mq   k =1 γ yy (mTvz ) =  p q 2 − ∑ ak γ yy (mTvz − kTvz ) + σ x ∑ bkh(kTvz − mTvz ) 0 ≤ m ≤ q  k =1 k =m  σ 2 (.Tvz ).1 + ∑ bk exp( − j2πfkTvz ) x odhad odhad  ARMA k =1 výkonového Pyy ( f ) = výkonového 2 p spektra spektra 1 + a exp( − j2πfkT ) q 2 ∑ k =1 k vz © Institut biostatistiky a analýz ARMA MODELY  citlivost na šum stejnosměrná složka a lineární trend znehodnocuje spektrum na nízkých kmitočtech – odstranit předem !!! zobecnění problému: yn = xn + wn a nechť je wn bílý šum s rozptylem σw2 a je nekorelovaný s xn; pak výkonové spektrum  Py ( f ) = Tvz .σ 2 x 2 i 1 + ∑ ap (i).e − j2 πfiTvz + σ 2 Tvz w  Py ( f ) = Tvz .σ 2 σ 2 + σ 2 A( z ).A * (1/ z*) .Tvz 2 x w + σ w Tvz = x A( z ).A * (1/ z*) A( z ).A * (1/ z*) !!! tohle už ale není přenosová funkce AR, alébrž ARMA !!! © Institut biostatistiky a analýz [ ] ARMA MODELY mnoho teoreticky optimálních postupů, ovšem bez praktického významu, protože:   jsou výpočetně příliš náročné (maticové operace, iterační optimalizační postupy); iterační optimalizace nezaručuje konvergenci řešení nebo konvergenci ke skutečně optimálnímu řešení; používají metody nejmenších čtverců ⇒ řešení soustavy lineárních rovnic; oddělený výpočet AR a MA parametrů proto suboptimální postupy :   © Institut biostatistiky a analýz ARMA MODELY METODA PRVNÍ  pro m>q je γ yy (mTvz ) = f (ak ), resp. f (ak ) dosadíme γ yy a řešíme p lineárních rovnic (pro modely vyšších řádů horší výsledky díky slabším odhadům γ yy(mTvz) pro větší m); přeurčená soustava lineárních rovnic (pro m>q) její řešení opět optimalizací metodou nejmenších čtverců © Institut biostatistiky a analýz ARMA MODELY METODA PRVNÍ předpokládejme, že známe odhady autokorelační funkce až do zpoždění M > p+q ryy (q − 1)  ryy (q)  r (q + 1) ryy (q)  yy     ryy (M − 1) ryy (M − 2)   ryy (q − p + 1)   a1   ryy (q + 1)  r (q + 2)  ryy (q − p + 2) a 2  .  = −  yy               ryy (M − p)  ap  ryy (M)       R yy a = −ryy protože Ryy je rozměru (M-q) x p a M-q > p, lze použít metodu nejmenších ;  a = −(R T .R yy )−1.R T .ryy výsledkem minimalizace je yy yy (může být použito i váhování k potlačení méně spolehlivých odhadů AK funkcí) © Institut biostatistiky a analýz ARMA MODELY METODA PRVNÍ p   −k A( z) = 1 + ∑ ak .z k =1  analyzovanou sekvenci vyfiltrujeme FIR filtrem A( z ) a dostaneme p  u(nTvz ) = y(nTvz ) + ∑ ak .y(nTvz − kTvz ), n = 0, 1, ..., N − 1; k =1 kaskádní zapojení ARMA(p,q) [HARMA(z)=B(z)/A(z)] s  AR(p) HAR ( z ) = A( z ) reprezentuje přibližně MA(q) proces s HMA(z)=B(z) © Institut biostatistiky a analýz ARMA MODELY METODA PRVNÍ Y( z ) B( z ) 1 H( z ) = = = B( z ) ⋅ X( z ) A( z ) A( z ) 1 Y( z ) = ⋅ U( z ) ⇒ U( z ) = A( z ) ⋅ Y( z ) A( z) © Institut biostatistiky a analýz ARMA MODELY METODA PRVNÍ z filtrované sekvence u(nTvz) pro p≤ n≤ N-1 je ruu(mTvz), odhad výkonového spektra potom je q  MA Pyy ( f ) = ∑ ruu (mTvz ) exp( − j2πfmTvz ) m= − q (opět můžeme „woknovat“ k potlačení méně spolehlivých odhadů autokorelačních funkcí nebo filtrace AR(q) oběma směry ⇒ dvě různé posloupnosti autokorelační funkce)  ARMA Pyy ( f ) = MA Puu ( f ) 1 + ∑ ak exp( − j2πfkTvz ) k =1 p 2 © Institut biostatistiky a analýz ARMA MODELY METODA DRUHÁ vstupní/výstupní identifikace metodou nejmenších  opakování: γ yy (mTvz ) = −∑ ak γ yy (mTvz − kTvz ) + ∑ bk γ xy (mTvz − kTvz ) k =1 k =m p q nelinearita vztahu pro výpočet γ yy(mTvz) plyne z toho, že neznáme γ xy(mTvz-kTvz), protože neznáme x(nTvz) y(nTvz ) = −∑ ak y(nTvz − kTvz ) + ∑ bk x(nTvz − kTvz ) k =1 k =0 p q © Institut biostatistiky a analýz ARMA MODELY METODA DRUHÁ q p y(nTvz ) = −∑ ak y(nTvz − kTvz ) + ∑ bk x(nTvz − kTvz ) + x(nTvz ), n = 0, 1, ...,N − 1 k =1 k =1 y = H.θ + ν y = [y(0) y(Tvz) … y(NTvz-Tvz)]T ν = [x(0) x(Tvz) … x(NTvz-Tvz)]T θ = [-a1 –a2 … -ap b1 b2 … bq ]T  y −1  y 0 H=      y N− 2  y −2 y −1     y −p y −p +1  x −1 x0  x −2 x −1  x N− 3 y N−3  y N−p −1 x N−2 x −q   x −q+1        x N−q−1    matice vstupních/výstupních dat o rozměru N x (p+q) © Institut biostatistiky a analýz ARMA MODELY METODA DRUHÁ minimalizace nejmenšími   Θ = (HT .H)−1.HT .y počáteční podmínky y-p, …, y-1, x-q, …, x-1 buď specifikovat nebo nulové odhad hodnot vstupního šumového signálu z analyzované posloupnosti AR modelem vysokého řádu © Institut biostatistiky a analýz ARMA MODELY METODA TŘETÍ (už tu částečně byla) ∞ B( z ) 1 = , kde C( z ) = 1 + ∑ c k z −k A( z ) C( z ) k =1 H( z ) = B( z ) 1 1 = = , A( z ) A( z ) C( z ) B( z ) koeficienty {ck} se určí nějakým AR algoritmem a z nich {ak} a {bk} M   −k C( z ) = 1 + ∑ c k z , M ≥ p + q k =1     B( z ) 1   je-li p>q, pak A( z) = C( z) ⇒ B( z) ⋅ C( z) = A( z) nebo v časové oblasti     ∑ bk .c n−k = an, n = 1,2,...; b0 = 1 q k =0 © Institut biostatistiky a analýz ARMA MODELY METODA TŘETÍ protože by mělo platit an = 0 pro n>p, je q   ∑ bk .c n−k = 0, pro n = p + 1, p + 2,..., p + q k =0       po určení {b1, b 2 ,..., b q } se spočítají {a1, a 2 ,..., ap } z   c p−1  cp   c cp  p +1      c p+ q−1 c p+ q−2      b1   c p+1   c p+1−q        c   c p+ 2−q  b 2   p+ 2  . =                c p  b q  c p+ q          ∑ bk .c n−k = an , pro n = 1, 2,..., p; c 0 = 1 q k =0 © Institut biostatistiky a analýz ARMA MODELY METODA TŘETÍ což je maticově   1 0  a1   c1    a  c 1  2  =  2 c1            ap  c p c p−1 c p −2      0  1    b   0  1           c p −q  b q    © Institut biostatistiky a analýz URČENÍ ŘÁDŮ ARMA MODELU 2(p + q) 2 AIC(p, q) = ln σ wpq + , N dodatečné kritérium: posouzení bělosti posloupnosti po inverzní filtraci analyzované posloupnosti navrženou ARMA soustavou  σ wpq je odhad rozptylu vstupní chybové posloupnosti odhad p  ryy (q) ryy (q − 1)  r (q + 1) ryy (q)  yy R' yy =     ryy (q + p − 1) ryy (q + p − 2)   ryy (q − p + 1)   ryy (q − p + 2)       ryy (q)   rozšířené modifikova né Y.-W. rovnice det(R’yy)=0, pokud je rozměr modelu větší než řád analyzovaného procesu © Institut biostatistiky a analýz