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 þ 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 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 ŘEŠENÍ PRO EXPERIMENTÁLNÍ DATA 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 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 levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz PODMÍNĚNOST SOUSTAVY LINEÁRNÍCH ROVNIC þ2x + 6y = 8 þ2x + 6,00001y = 8,00001 þ þŘešení: ? þ levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz PODMÍNĚNOST SOUSTAVY LINEÁRNÍCH ROVNIC þ2x + 6y = 8 þ2x + 6,00001y = 8,00001 þ þŘešení: x = 1; y = 1 þ þ levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz PODMÍNĚNOST SOUSTAVY LINEÁRNÍCH ROVNIC þ2x + 6y = 8 þ2x + 6,00001 = 8,00001 þ þŘešení: x = 1; y = 1 þ þ2x + 6y = 8 þ2x + 5,99999y = 8,00002 þ þŘešení: ? levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz PODMÍNĚNOST SOUSTAVY LINEÁRNÍCH ROVNIC þ2x + 6y = 8 þ2x + 6,00001 = 8,00001 þ þŘešení: x = 1; y = 1 þ þ2x + 6y = 8 þ2x + 5,99999y = 8,00002 þ þŘešení: x = 10; y = -2 levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz þšpatně podmíněná soustava: þnormalizace matice soustavy tak, že všechny prvky rozšířené matice soustavy jsou menší (ale ne moc) než 1 þšpatně podmíněná soustava lineárních rovnic je taková, která po této normalizaci má inverzní matici soustavy A-1 s některými prvky o „velké“ absolutní hodnotě. PODMÍNĚNOST SOUSTAVY LINEÁRNÍCH ROVNIC levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz PODMÍNĚNOST SOUSTAVY LINEÁRNÍCH ROVNIC levy-panel-IBA-se-zavojem logo-IBA-transparent logo-MU © Institut biostatistiky a analýz PODMÍNĚNOST SOUSTAVY LINEÁRNÍCH ROVNIC 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 þ* 18. srpna 1919, V17. července 2002, White Plains, N.Y., U.S.A. þamerický teoretický fyzik, jaderný inženýr, informatik; þVzdělání: þBSc. (strojní inženýrství) & MSc. (fyzika) , Stevens Inst. Technology, Hoboken, New Jersey, U.S.A. þPh.D. (teoretická fyzika) New York Univ. þZaměstnání: þManhattan projekt (Nuclear Development Corporation of America) þzakladatel Sage Instruments, Kalifornie (testo- cí přístroje a va- þ vací přístroje a přístroje pro lékařský výzkum) þIBM, Advanced Systems Development Division (28 let) – systémy automatického řízení, komprese dat, digitální tiskárny þVýsledky: þGoertzelův algoritmus, 13 U.S.patentů (komprese dat, zpracování 2D dat, dětský inkubátor, …),… þ GERALD GOERTZEL Some Mathematical Methods of Physics (Dover Books on Physics) 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ě