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 IV. ROZKLAD POMOCÍ OPTIMÁLNÍ APROXIMACE levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz ROZKLAD POMOCÍ OPTIMÁLNÍ APROXIMACE skenování0003.jpg •měsíční průměry celkové koncentrace ozónu, 65°S – 65°N levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz HRÁTKY S POČÁTEČNÍ FÁZÍ •originál •φ01=φ02=π/2 •φ01= π/4; φ02=π/2 •φ01=φ02=π levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz ROZKLAD POMOCÍ OPTIMÁLNÍ APROXIMACE skenování0003.jpg •měsíční průměry celkové koncentrace ozónu, 65°S – 65°N levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz PRINCIP þdominantní část časové řady lze vyjádřit ve tvaru þs(n) = A0 + A.cos(2pfn + j0), þkde f=1/12 cyklů/měsíc = 1 cyklus/rok þ(Tvz je 1 měsíc) þneznámé jsou A0, A a j0 þexperimentální data x(n) lze vyjádřit þx(n) = s(n) + e(n), þe(n) je reziduum n-tého vzorku þ(model je dobrý, pokud jsou rezidua malá) levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz PRINCIP þmetoda nejmenších čtverců þ þ þnejlépe když je splněn požadavek na linearitu vzhledem k parametrům modelu þs(n) = A0 + Cc.cos(2pfn) + Cs.sin(2pfn), þkde Cc = A.cos(j0) a Cs = -A.sin(j0) þnebo také þ levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz PRINCIP levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz PRINCIP þrovnice pro určení A0, Cc a Cs metodou nejmenších čtverců: levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz skenování0004.jpg ŘEŠENÍ PRO EXPERIMENTÁLNÍ DATA levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz VÍCE HARMONICKÝCH SLOŽEK þx(n) = A0 + Cc1.cos(2pfn) + Cs1.sin(2pfn) + þ+ Cc2.cos(2.2pfn) + Cs2.sin(2.2pfn) + e(n) þ levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz DVĚ HARMONICKÉ SLOŽKY - ŘEŠENÍ •? levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz DVĚ HARMONICKÉ SLOŽKY - ŘEŠENÍ skenování0005.jpg levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þCO KDYŽ NEZNÁME FREKVENCE? þ þ þ þ þ þ þ þ þ21 vrcholů během 600 dní Þ perioda 600/21 » 28,6 dní Þ Þ frekvence f = 1/28,6 » 0,035 cyklů/den ROZKLAD POMOCÍ OPTIMÁLNÍ APROXIMACE skenování0006.jpg •půlnoční magnituda proměnné hvězdy v 600 následných nocích (podle Whittaker,E.T., Robinson G.: The Calculus of Observations, London, Blackie & Son 1944, str.349) •http://www.york.ac.uk/depts/maths/data/ts/, No.26 levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þCO KDYŽ NEZNÁME FREKVENCE? þto je ale jenom odhad ROZKLAD POMOCÍ OPTIMÁLNÍ APROXIMACE skenování0007.jpg •součet čtverců pro data z proměnné hvězdy pro f Îá0,03; 0,04ñ cyklů/den •maximum pro = 0,034442 = 1/29,0343 cyklů/den •vedlejší laloky v grafu nesignalizují přítomnost dalších periodických komponent levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þCO KDYŽ NEZNÁME FREKVENCE? ROZKLAD POMOCÍ OPTIMÁLNÍ APROXIMACE skenování0008.jpg •určená kosinusovka pro periodu 29,034 dní a její zbytková funkce 25 vrcholů během 600 dní Þ perioda 600/25 = 24 dní Þ Þ frekvence f = 1/24 » » 0,042 cyklů/den levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þCO KDYŽ NEZNÁME FREKVENCE? þanalýzou originální časové řady lze ekvivalentně þ(založeno na modelu þx(n)=μ+A2cos2πf2n+B2sin2πf2n+ε(n) ) þurčit frekvenci blízkou 0,042 cyklů/den þ= 0,041724 = 1/23,967 þ ROZKLAD POMOCÍ OPTIMÁLNÍ APROXIMACE levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þCO KDYŽ NEZNÁME FREKVENCE? þobecně model pro dvě frekvenční složky je þx(n) = μ + A1cos2πf1n + B1sin2πf1n + A2cos2πf2n + B2sin2πf2n + ε(n) þkriteriální funkce pro minimalizaci þ þ ROZKLAD POMOCÍ OPTIMÁLNÍ APROXIMACE skenování0009.jpg levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þCO KDYŽ NEZNÁME FREKVENCE? þ þ ROZKLAD POMOCÍ OPTIMÁLNÍ APROXIMACE skenování0010.jpg levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz ÚŽASNÁ DATABÁZE ČASOVÝCH ŘAD þhttp://www.york.ac.uk/depts/maths/data/ts/welcome.htm logo-IBA logo-MU © Institut biostatistiky a analýz V. GOERTZELŮV ALGORITMUS levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz GOERTZELŮV ALGORITMUS (GA) þefektivní algoritmus pro výpočet jedné hodnoty (vzorku) diskrétní Fourierovy transformace; þ(pro celé spektrum je GA složitější než FFT, ale pro výpočet omezeného počtu vzorků je efektivnější) þpodobně jako definiční vztah DFT počítá GA parametry jedné určité frekvenční složky analyzované časové řady (diskrétního signálu); þna rozdíl od DFT pro posloupnost reálných čísel používá pouze reálnou aritmetiku þ levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz GOERTZELŮV ALGORITMUS (GA) þALGORITMUS þ(protože je podstatou číslicovým filtrem, nazývá se také často Goerzelovým filtrem) þmá dva sériově zapojené stupně: þ(1): s(n) = x(n) + 2cos(2pf)s(n-1) – s(n-2) þx(n) = 0 pro n<0; s(-2) = s(-1) = 0 þ(2): y(n) = s(n) – e-2pjf .s(n-1) þhodnota f určuje hodnotu normalizované (počet cyklů na vzorek) frekvence analyzované harmonické složky þ!!! výpočet v reálném čase !!! levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz GOERTZELŮV ALGORITMUS (GA) þALGORITMUS þpřenosové funkce obou dílčích stupňů: þ þ þpóly tohoto filtru leží na e+j2pf a e-j2pf, tj. na jednotkové kružnici na frekvenci odpovídající f, tedy je na mezi stability a jeho stabilita může být tak závislá na numerických chybách, zejména při výpočtech s dlouhou vstupní posloupností a s aritmetikou s nižší přesností levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz GOERTZELŮV ALGORITMUS (GA) þALGORITMUS þcelková přenosová funkce pak je þ þ þ þpo zpětném vyjádření v časové doméně je (po rozvinutí rekurzivního vztahu pro kauzální posloupnost) þ þ þ (@) levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz GOERTZELŮV ALGORIMUS - VÝPOČET þPŘEDPOKLADY þfiltrace končí posledním N-tým vzorkem; þhodnoty frekvence, pro které se realizuje výpočet (DFT) jsou dány vztahem þ þ kde k je celé číslo kÎ{0, 1, …, N-1}. þPo dosazení do (@) je þ þ þ þcož je vztah lišící se od definice DFT pouze horní mezí součtu. levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz GOERTZELŮV ALGORIMUS - VÝPOČET þPokud vložíme x(N)=0, lze beztrestně snížit horní mez na N-1. þZ definiční rekurzivní diferenční rovnice (1) je při nulovém posledním vzorku x(N)=0 þs(N) = 2cos(2pf)s(N-1) – s(N-2). (*) þZ toho je výsledný algoritmus: þukončit výpočet podle vztahu (1) pro x(N-1); þpro výpočet s(N) z hodnot s(N-1) a s(N-2) použít vztah (*); þvýslednou hodnotu y(N) z s(N) a s(N-1) určit ze vztahu (2). þ levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz GOERTZELŮV ALGORIMUS - VÝPOČET þPoslední dva kroky algoritmu lze zjednodušit þ þ þ þ þ þ levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz GOERTZELŮV ALGORIMUS - APLIKACE þVZOREK VÝKONOVÉHO SPEKTRA þ þ þ þ þ þ levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz GOERTZELŮV ALGORIMUS - APLIKACE þVZOREK DFT POMOCÍ REÁLNÉ ARITMETIKY þprvní část algoritmu pracuje pouze s reálnými čísly (pokud jsou hodnoty vstupní posloupnosti reálné), komplexní výpočet reprezentuje pouze poslední krok algoritmu þ þcož můžeme rozepsat podle Eulerova vztahu þ þ þpotom lze snadno určit i modul a fázi komplexního výsledku þje-li vstup komplexní, lze jej rozložit na reálnou a imaginární část a GA zpracuje obě části separátně