logo-IBA logo-MU © Institut biostatistiky a analýz SPEKTRÁLNÍ ANALÝZA ČASOVÝCH ŘAD prof. Ing. Jiří Holčík, CSc. holcik@iba.muni.cz logo-IBA logo-MU © Institut biostatistiky a analýz VI. SPEKTRÁLNÍ ANALÝZA POMOCÍ METODY VLASTNÍCH ČÍSEL levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz ZAČÍNÁME þAR(p) proces znehodnocený aditivním bílým šumem º ARMA(p,p) proces; þteď bude … periodický signál + bílý šum þx(nTvz) = 2.cos(2πfkTvz).x(nTvz-Tvz) – x(nTvz-2Tvz) þtento systém generuje signál þx(nTvz) = 2.cos(2πfknTvz) pro n ≥ 0 þpokud jsou počáteční podmínky x(-1)=-1 a x(-2)=0 levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz ZAČÍNÁME þobecně, signál skládající se z p harmonických složek splňuje diferenční rovnici þ þ ([) þ þcož odpovídá systému s přenosovou funkcí þ þ (N) þ þ þ(Polynom A(z)=1+∑amz-m má 2p kořenů na jednotkovém kruhu v místech, která odpovídají frekvencím harmonického signálu.) levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz PRINCIP þpřepokládejme periodický signál + bílý šum w(nT) {E(|w(nTvz)2|)=σw2} þy(nTvz) = x(nTvz) + w(nTvz) þpo dosazení za x(nT) z tohoto vztahu do ([) máme þ þ þ þ þcož představuje ARMA proces s identickými AR i MA parametry þyT.a = wT.a (õ) þyT =[y(nTvz), y(nTvz-Tvz),…,y(nTvz-2pTvz)], þwT =[w(nTvz),w(nTvz-Tvz),…,w(nTvz-2pTvz)], þaT =[1, a1,…, a2p-1,a2p], þ þ levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz PRINCIP þvynásobením obou stran (õ) vektorem y a určením střední hodnoty þE(y.yT).a = E(y.wT).a = E((x+w).wT).a þΓyy.a = 0 + σw2.a þ (Γyy – σw2.I).a = 0 … vlastní (charakteristická) rovnice þσw2 je vlastní číslo autokorelační matice Γyy; þa je vlastní vektor Γyy spojený s vlastním číslem σw2; levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz PISARENKOVA HARMONICKÁ DEKOMPOZICE þMějme p náhodně fázově posunutých harmonických signálů s aditivním bílým šumem. þHodnoty autokorelační funkce jsou þ þ þ je þ průměrný výkon i-té sinusovky, Ai je její amplituda levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þmaticově þ þ þ þ þznáme-li frekvence fi, 1≤i≤p, můžeme spočítat výkon jednotlivých harmonických složek, místo hodnot gyy(mTvz) použijeme odhady ryy(mTvz), známe-li výkony, určíme rozptyl šumu þ PISARENKOVA HARMONICKÁ DEKOMPOZICE levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þpodle Pisarenka platí pro ARMA proces obsahující p harmonických složek v aditivním bílém šumu, že rozptyl σw2 odpovídá minimálnímu vlastnímu číslu autokorelační matice, pokud je rozměr autokorelační matice větší nebo roven þ(2p+1) x (2p+1). þPotom požadovaný vektor koeficientů ARMA modelu je dán vlastním vektorem náležejícím minimálnímu vlastnímu číslu. þFrekvence fi, i=1,…,p se určí řešením rovnice, dané položením jmenovatele ve vztahu (N) rovno nule, kde koeficienty am jsou určeny vlastním vektorem spojeným s minimálním vlastním číslem. þ PISARENKOVA HARMONICKÁ DEKOMPOZICE levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz VLADIMIR FEDOROVICH PISARENKO þ? þPh.D. þMoskevská státní univerzita 1963 þdisertační práce: þMatematická klasifikace objektů þobor: Teorie pravděpodobnosti a stochastické procesy þškolitel: Roland Lvovich Dobrushin Vladimer Pisarenko, V. F. The retrieval of harmonics from a covariance function Geophysics, J. Roy. Astron. Soc., vol. 33, pp. 347-366, 1973. View works by V. F Pisarenko 1960 - 1960 View works by V. F Pisarenko 1961 - 1961 View works by V. F Pisarenko 1962 - 1962 View works by V. F Pisarenko 1963 - 1963 View works by V. F Pisarenko 1964 - 1964 View works by V. F Pisarenko 1965 - 1965 View works by V. F Pisarenko 1966 - 1966 View works by V. F Pisarenko 1967 - 1967 View works by V. F Pisarenko 1968 - 1968 on View works by V. F Pisarenko 1970 - 1970 View works by V. F Pisarenko 1971 - 1971 View works by V. F Pisarenko 1972 - 1972 View works by V. F Pisarenko 1973 - 1973 View works by V. F Pisarenko 1974 - 1974 View works by V. F Pisarenko 1975 - 1975 View works by V. F Pisarenko 1976 - 1976 View works by V. F Pisarenko 1977 - 1977 View works by V. F Pisarenko 1978 - 1978 View works by V. F Pisarenko 1979 - 1979 View works by V. F Pisarenko 1980 - 1980 on View works by V. F Pisarenko 1982 - 1982 View works by V. F Pisarenko 1983 - 1983 View works by V. F Pisarenko 1984 - 1984 on on on on on View works by V. F Pisarenko 1990 - 1990 on View works by V. F Pisarenko 1992 - 1992 View works by V. F Pisarenko 1993 - 1993 on View works by V. F Pisarenko 1995 - 1995 View works by V. F Pisarenko 1996 - 1996 View works by V. F Pisarenko 1997 - 1997 View works by V. F Pisarenko 1998 - 1998 View works by V. F Pisarenko 1999 - 1999 View works by V. F Pisarenko 2000 - 2000 View works by V. F Pisarenko 2001 - 2001 View works by V. F Pisarenko 2002 - 2002 View works by V. F Pisarenko 2003 - 2003 View works by V. F Pisarenko 2004 - 2004 View works by V. F Pisarenko 2005 - 2005 View works by V. F Pisarenko 2006 - 2006 on View works by V. F Pisarenko 2008 - 2008 View works by V. F Pisarenko 2009 - 2009 on View works by V. F Pisarenko 2011 - 2011 View works by V. F Pisarenko 2012 - 2012 View works by V. F Pisarenko 2013 - 2013 View works by V. F Pisarenko 2014 - 2014 View works by V. F Pisarenko 2015 - 2015 View works by V. F Pisarenko 2016 - 2016 View works by V. F Pisarenko 2017 - 2017 View works by V. F Pisarenko 2018 - 2018 View works by V. F Pisarenko 2019 - 2019 View works by V. F Pisarenko 2020 - 2020 levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þPředpokládejme hodnoty AKF gyy(0)=3, gyy(1)=1 a gyy(2)=0. Proces obsahuje jeden harmonický signál v bílém šumu. Určete jeho frekvenci, výkon a rozptyl, tj. výkon šumu. þŘešení: þautokorelační matice je PISARENKOVA HARMONICKÁ DEKOMPOZICE PŘÍKLAD levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þjejí minimální vlastní číslo je rovno nejmenšímu kořenu charakteristického polynomu þ þ þ þl1 = 3, l2 = 3+Ö2, l3 = 3-Ö2 Þ sw2 = lmin = 3-Ö2 þodpovídající charakteristický vektor má složky a0=1, a1, a2, pro které platí þ PISARENKOVA HARMONICKÁ DEKOMPOZICE PŘÍKLAD levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þřešením získáme a1 = -Ö2 a a2 = 1 þz2 - Ö2 + 1 = 0 þ þ tj. leží na jednotkové kružnici þ þ þvýkon P1.cos(2пf1Tvz) = gyy(1)=1 Þ P1 = Ö2, þproto þ þkontrola: sw2 = gyy(0) - P1 = 3-Ö2, což souhlasí s lmin. PISARENKOVA HARMONICKÁ DEKOMPOZICE PŘÍKLAD levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þzbývá určit frekvence fi, 1≤i≤p, - to vyžaduje znalost vlastního vektoru a pro vlastní číslo sw2; þprotože jsme si již uvedli, že: þpro ARMA proces skládající se z p sinusovek v aditivním bílém šumu je rozptyl sw2 roven nejmenšímu vlastnímu číslu matice Γyy za předpokladu, že rozměr Γyy ≥ (2p+1)´(2p+1) þpožadovaný vektor koeficientů ARMA systému je dán hodnotami vlastního vektoru odpovídajícímu nejmenšímu vlastnímu číslu – frekvence fi, 1≤i≤p, pak odpovídají kořenům polynomu, jehož koeficienty PISARENKOVA HARMONICKÁ DEKOMPOZICE PŘÍKLAD levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz PRONYHO METODY þGaspard Clair François Marie Riche, Baron de Prony þ(22.7.1755 – 29.7.1839) þzákony (rozumějme funkce) popisující expanzi plynů lze vyjádřit součtem exponenciál þß þnavrhnul metodu na interpolaci naměřených hodnot na základě exponenciálního modelu interpolační funkce Gaspard de Prony Essai éxperimental et analytique: sur les lois dilatabilité de fludes élastique et sur celles de la force expansive de la vapuer de l’alkool, à différentes températures. Journal de l’École Polytechnique, Floréal et Plairial, an III (1795), vol 1, cahier 22, 24-76 levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz PRONYHO METODY þpředpokládejme, že vzorky signálu neobsahují šum a skládají se z p exponenciálních signálů þ þ þproblém je určit Ak a sk z daných N vzorků. þPotřebujeme nejméně 2p hodnot x(nTvz) – problém je bohužel nelineární. þTo co Prony vymyslel je, že výpočet Ak a sk může být vzájemně oddělen a realizován řešením dvou soustav lineárních rovnic + řešení polynomiální rovnice. levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz PRONYHO METODY þAk i sk jsou obecně komplexní a je-li x(nTvz) reálné, þpak þ þ þ þprotože součet komplexních exponenciál je homogenním řešením lineární diferenční rovnice, tak musí taková diferenční rovnice existovat levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þpro koeficienty diferenční rovnice je þ þ þ þa to všechno vede k definici tří kroků Pronyho metody: þ1. konstrukce a řešení lineární soustavy rovnic pro koeficienty bk; PRONYHO METODY levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þ2. po určení koeficientů bk stanovení kořenů rovnice þ þ þ k-tý kořen je roven þ3. z definiční rovnice þ þ pak můžeme definovat lineární soustavu p rovnic, ze které vypočítáme amplitudy Ak; þ þA TO JE ZÁKLADNÍ MYŠLÉNKA BARONA PRONYHO þ PRONYHO METODY levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz PRONYHO METODY þz-transformace dané posloupnosti x(nTvz) je þ þ þvýměnou součtů a sečtením geometrické řady je þ þ þ þ þ (\1) levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz PRONYHO METODY þX(z) lze tedy vyjádřit pomocí poměru dvou polynomů þ þ (\2) þ þsrovnáním obou (\) vztahů lze říci, že B(z) je polynom p-tého řádu s kořeny exp(s1Tvz), exp(s2Tvz), …, exp(spTvz); þ þ þpo zjednodušení čitatele C(z) je polynom (N+p-1)-ho řádu definovaného levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz PRONYHO METODY þzatímco některé koefcienty polynomu C(z) (c0, …, cp-1 a cN, …, cN+p-1) závisí na neznámých amplitudách a pólech (nebo frekvencích), jiné koeficienty cp, …, cN-1 º 0; þprotože X(z).B(z) = C(z), je v časové oblasti þx(nTvz)Äbn = cn, þcož znamená levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þprotože dále jsou cp, …, cN-1 nulové, středních N-p rovnic pro nulové pravé strany lze vyjádřit þ þ (ä) þ þ þXf.b = 0 þkde i,j-tý prvek matice Xf je xf i,j = x(pTvz+iTvz-jTvz), i=1,2,…, N-p; j = 1,2,…, p+1 a b=[1,b1,b2,…, bp]T þJe-li N=2p, pak máme právě dost rovnic k tomu určit b1,b2,…, bp. PRONYHO METODY levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz PRONYHO METODY þJakmile určíme bk, známe polynom þB(z)=1+ b1z-1+ b2z-2+…+ bpz-p þa řešením rovnice B(z)=0 dostaneme kořeny exp(s1Tvz), exp(s2Tvz), …, exp(spTvz);. þKoeficienty Ak jsou pak určeny řešením následujících p lineárních rovnic þ þ (å) þ levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz PRONYHO METODY þSUMARIZACE ALGORITMU þ 1.vypočítat koeficienty bk z (ä); 2.spočítat kořeny polynomu B(z) a tak určit þ 3.určit amplitudy A1,…, Ap z rovnice (å); levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz PRONYHO METODA NEJMENŠÍCH ČTVERCŮ þpředpokládáme signál + aditivní šum þNechť diskrétní funkce þ þ (†) þ þje modelem signálu, reprezentovaného naměřenými hodnotami x0,…, xN-1. þHodnoty bk, zk jsou obecně komplexní þ þ þAk je amplituda, Φk počáteční fáze, σk koeficient tlumení a fk frekvence oscilací, Tvz je vzorkovací perioda levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þcílem je nalézt {Ak,Φk,σk,fk} a p takové, že je minimalizována chyba þ þ þje to těžce nelineární problém nejmenších čtverců, který lze řešit iteračně postupným zlepšováním počátečního odhadu. þpokusme se aplikovat Pronyho postup, poskytující suboptimální řešení, které sice nezaručuje nalezení globálního minima, ale poskytuje dostatečně přijatelné řešení PRONYHO METODA NEJMENŠÍCH ČTVERCŮ levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þtedy ještě jednou a trochu jinak: þvztah (†) představuje homogenní řešení lineární diferenční rovnice s konstantními parametry. Jakými? To určíme! þdefinujme polynom þ þ þz (†) jeden způsob jak vyjádřit odhad je PRONYHO METODA NEJMENŠÍCH ČTVERCŮ levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þvynásobením am a sečtením přes posledních p+1 součinů þ þ þpo substituci je þ þ þ þTa nula plyne z toho, že poslední suma je právě hodnota polynomu Y(zs), tj. pro jeden jeho kořen. þ ß PRONYHO METODA NEJMENŠÍCH ČTVERCŮ levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þ þpo dosazení za je þ þ þ þ þ þ þto nám poskytuje alternativní model (součet exponenciál + aditivní šum) pomocí ARMA systému AR a MA parametry þna rozdíl od Pisarenka nejsou koeficienty ai omezeny tak, aby kořeny měly jen jednotkový modul (automaticky jen harmonické složky) PRONYHO METODA NEJMENŠÍCH ČTVERCŮ a co takhle Pisarenko?! levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þ þminimalizace vede k soustavě nelineárních þ þrovnic, proto alternativa (rozšířená Pronyho metoda): þkdyž þ þpak þ þ þto pak vede k minimalizaci en Þ pouze AR model Þ linearita (en je určena MA procesem řízeným chybou aproximace en) PRONYHO METODA NEJMENŠÍCH ČTVERCŮ levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þen je rozdíl mezi xn a jeho lineární predikcí p-tého řádu, zatímco þen je rozdíl mezi xn a jeho exponenciální aproximací þřád p je určen nějakým způsobem pro určení řádu AR procesu (modelu) þ PRONYHO METODA NEJMENŠÍCH ČTVERCŮ levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þurčíme koeficienty am AR procesu Þ najdeme kořeny AR polynomu a problém se redukuje na řešení soustavy lineárních rovnic s neznámými koeficienty bm þ þkde þ þ þ þ þ þminimalizace nejmenšími oo vede na řešení þ (T) PRONYHO METODA NEJMENŠÍCH ČTVERCŮ levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þpři řešení pomůže, když víme, že þ þ þ þ þ þparametry Ai, Φi, αi a fi pak určíme þAi =|bi | þΦ i = arctg(Imbi/Rebi) þαi =ln|zi |/Tvz þfi = arctg(Imzi/Rezi)/2πTvz PRONYHO METODA NEJMENŠÍCH ČTVERCŮ levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þPronyho metoda užitečná při analýze přechodných dějů þomezíme-li se na tlumené reálné sinusovky, pak þ þ (') þ þpro reálné x(t) požadujeme komplexně sdružené þ a ; þ PRONYHO METODA NEJMENŠÍCH ČTVERCŮ dále předpokládáme, že koeficienty tlumení jsou záporné, tj. exponenciály jsou tlumené (je-li α=0, jsou sinusovky tak jak mají být) levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þprotože vztah (') má v tom případě konečnou energii, jeho spektrální hustota je rovna FT tohoto vztahu þ þkde þ þ þspektrum je lineárně úměrné energii, nikoliv jako u AR modelů, kde je úměrné (nelineárně) výkonu PRONYHO METODA NEJMENŠÍCH ČTVERCŮ lze produkovat spektra s úzkými či širokými laloky – závisí na koeficientech tlumení (šířka pásma pro -3dB je a/p [Hz]) levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þPROBLÉMY þpočet exponenciál ~ řád AR systému – protože je třeba určit 2p parametrů, max. řád by měl být pmax£N/2, zatímco u AR modelů je možné p>N/2; þpřítomnost šumu ovlivňuje přesnost Pronyho odhadů; þšum může způsobit, že koeficienty tlumení jsou moc velké; þ PRONYHO METODA NEJMENŠÍCH ČTVERCŮ levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz skenování0004.jpg levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þVÝHODY OPROTI PISARENKOVI þnejsou třeba odhady autokorelačních funkcí; þvyskytuje se méně falešných spektrálních čar díky možnému lepšímu odhadu řádu (?); þodhady frekvencí a výkonů mají menší chyby; þjednodušší výpočet (dvě lineární soustavy + nalezení kořenů polynomu x výpočet komplexních vlastních čísel) þ þ A JE TO! þ PRONYHO METODA PRO VÝPOČET SPEKTRÁLNÍCH ČAR http://1.bp.blogspot.com/_5eb1cI127Ig/SutdOHvbd4I/AAAAAAAAAPY/BueUE7ykBSQ/s320/pat+a+mat.gif