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(nT) = 2.cos(2πfkT).x(nT-T) – x(nT-2T) þtento systém generuje signál þx(nT) = 2.cos(2πfknT) pro n ≥ 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í þ þ þ þ þ(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(nT)2|)=σw2} þy(nT) = x(nT) + w(nT) þ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(nT), y(nT-T),…,y(nT-2pT)], þwT =[w(nT),w(nT-T),…,w(nT-2pT)], þ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 – σ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(mT) použijeme odhady ryy(mT), 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 þ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 - Ö2z + 1 = 0 þ þ tj. leží na jednotkové kružnici þ þ þvýkon P1.cos(2пf1T) = gyy(1)=1 Þ P1 = Ö2, þprotože þ þ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 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 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(nT) – 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 þz-transformace dané posloupnosti x(nT) 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 es1.T, es2.T,…, esp.T; þ þ þ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(nT)Ä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(pT+iT-jT), 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 es1.T, es2.T,…, esp.T. þKoeficienty Ak jsou pak určeny řešením následujících p lineárních rovnic þ þ (å) þ þA TO JE ZÁKLADNÍ MYŠLÉNKA BARONA PRONYHO 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 bm, zm jsou obecně komplexní þ þ þAm je amplituda, Θm počáteční fáze, αm koeficient tlumení a fm frekvence oscilací, T je vzorkovací perioda levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þcílem je nalézt {Am,Θm,αm,fm} 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. þProny navrhnul 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 þještě jednou a 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. þ ß þa co takhle þPisarenko? 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 PRONYHO METODA NEJMENŠÍCH ČTVERCŮ 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: þ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 |/T þfi = arctg(Imzi/Rezi)/2πT 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é; þ þ þ A JE TO! PRONYHO METODA NEJMENŠÍCH ČTVERCŮ http://1.bp.blogspot.com/_5eb1cI127Ig/SutdOHvbd4I/AAAAAAAAAPY/BueUE7ykBSQ/s320/pat+a+mat.gif