Wiesław Kosek
Centrum Badań Kosmicznych
PAN, Bartycka 18A, 00-716 Warszawa
kosek@cbk.waw.pl
http://www.cbk.waw.pl/~kosek
tel. 8403766 wew. 306
Wprowadzenie
Od lat 60-tych wraz z rozwojem elektronicznej techniki obliczeniowej nastąpił szybki rozwój teorii analiz szeregów czasowych oraz rozwój związanych z nimi metod numerycznych. Niektóre z tych metod wcześniej nie mogły być zastosowane, ze względu na ogromną czasochłonność obliczeń spowodowaną licznymi sumowaniami i mnożeniami występującymi przy wyznaczaniu splotów, autokowariancji, transformat oraz numerycznych wartości całek. Operacje takie mają zastosowanie przy wyznaczaniu widm mocy, filtrowaniu oraz prognozowaniu szeregów czasowych. W analizie szeregów czasowych wyróżniane są następujące zagadnienia:
Zadaniem filtracji jest usunięcie z szeregu czasowego oscylacji o określonych częstotliwościach. Szczególnym przypadkiem filtracji danych jest wygładzanie, które polega na usuwaniu z danych składowych o wysokich częstotliwościach a pozostawianie tych, które odpowiadają częstotliwościom niskim. Filtr jest urządzeniem lub procesem fizycznym, który działając na szereg czasowy modyfikuje jego własności w ściśle określony sposób.
Analizy widmowe mają za zadanie wyznaczenie funkcji częstotliwości mówiącej o tym, które składowe częstotliwościowe w szeregu mają największą moc. Moc ta utożsamiana jest zwykle z kwadratem amplitudy oscylacji o danej częstotliwości.
Zadaniem analiz czasowo częstotliwościowych jest przeprowadzanie analizy widmowej ograniczonej, w krótkim przedziale czasowym gdy zachodzi podejrzenie, że amplitudy lub fazy oscylacji procesu stochastycznego zmieniają się w czasie.
W analizie szeregów czasowych głównym problemem jest znalezienie takich metod analiz widmowych, filtracji oraz prognozowania, które zminimalizują błąd wyznaczanych estymatorów widma mocy, przefiltrowanego szeregu czasowego oraz prognozy. Błąd ten wynika z ograniczonej liczby danych szeregu czasowego lub występowania szumu pomiarowego. Częściową eliminację błędu wynikającego z ograniczonej liczby danych uzyskuje się np. poprzez zastosowanie okien wagowych lub wydłużenie szeregu czasowego jeżeli jest to możliwe. Eliminacje błędów związanych z występowaniem szumu pomiarowego uzyskuje się poprzez wykorzystanie estymatorów autokowariancji.
I. PROCES STOCHASTYCZNY
Procesem stochastycznym
nazywamy szereg liczb przypisanych pewnym określonym momentom czasu
.
Z wartości tego procesu nie można jednoznacznie wyznaczyć jego wartości
poza tymi momentami czasu, w których jest określony, ze względu na występujący
w tym procesie czynnik losowy (w odróżnieniu od deterministycznej funkcji
ciągłej
, której wartości
można określić dla każdego argumentu
).
Procesy stochastyczne charakteryzowane są następującymi parametrami:
Dla stacjonarnych przebiegów ciągłych o
zerowej wartości oczekiwanej funkcja autokowariancji określona jest następującym
wzorem:
W przypadku niestacjonarnego przebiegu
ciągłego jego funkcja autokowariancji jest funkcją nie tylko przesunięcia
czasowego
ale również czasu
.
Autokowariancja jest zawsze funkcją parzystą co oznacza, że
.
W praktyce dysponujemy skończonym czasem
trwania procesu ciągłego zatem estymator autokowariancji dla przebiegów
ciągłych zdefiniowany jest następująco:
Dla stacjonarnych szeregów czasowych dyskretnych
o zerowej wartości oczekiwanej nieobciążony i obciążony estymator autokowariancji
przedstawiają się następującymi równaniami:
,
Różnica pomiędzy estymatorem autokowariancji
obciążonym i nieobciążonym jest taka, że dla nieskończenie długiego szeregu
czasowego
jego nieobciążony
estymator autokowariancji równy jest funkcji autokowariancji nieskończonego
przebiegu ciągłego
,
z którego szereg ten został próbkowany.
- Autokorelacja
Autokorelacja jest to znormalizowana autokowariancja
powstała przez podzielenie autokowariancji przez wariancję. Autokorelacja
ma taką samą charakterystykę częstotliwościowa co autokowariancja a jej
wartości zmieniają się od –1 do 1.
- Wariancja
Wariancja nazywana czasem momentem drugiego
rzędu jest szczególnym przypadkiem autokowariancji dla przesunięcia czasowego
.
Dla stacjonarnych przebiegów ciągłych
o
zerowej wartości oczekiwanej wariancja określona jest następującym wzorem:
W praktyce dysponujemy skończonym czasem
trwania procesu stochastycznego zatem estymator wariancji dla przebiegów
ciągłych zdefiniowany jest następująco:
Dla stacjonarnych szeregów czasowych dyskretnych
o zerowej wartości oczekiwanej estymator wariancji określony jest następującym
równaniem:
- odchylenie standardowe
Odchylenie standardowe definiowane jest jako pierwiastek kwadratowy
z wariancji.
Powyższe równanie dla przebiegu ciągłego
o
zerowej wartości oczekiwanej można przedstawić następująco:
Równanie to jest jedną z wielu wersji
zapisu twierdzenia Wienera-Chinczyna, które wyraża się wzorem:
Widmo mocy przebiegu ciągłego jest więc
transformatą Fouriera funkcji autokowariancji. Skończony czas trwania procesu
stochastycznego, a także występowanie szumu w postaci błędu pomiarowego
wpływają na błąd wyznaczenia estymatora widma mocy
obliczanego
za pomocą analizy widmowej. Zadaniem analizy widmowej jest wyznaczenie
estymatora
w taki sposób
aby informację o występujących częstotliwościach oscylacji podać możliwie
najdokładniej. Informacja ta wygląda w ten sposób, że częstotliwościom
procesu stochastycznego przypisane są maksima, których wysokości na ogół
są proporcjonalne do kwadratów średnich amplitud oscylacji znajdującychsię
w tym procesie.
Stacjonarność procesu stochastycznego
Proces stochastyczny jest stacjonarny
w szerszym sensie jeżeli
,
,
a także jego funkcja autokowariancji
zależy
tylko od przesunięcia czasowego
i
nie zależy od czasu
.
Niestacjonarność szeregów czasowych może być spowodowana zmianą w czasie ich charakterystyki częstotliwościowej, wariancji lub wartości średniej. W przypadku niestacjonarności spowodowanej zmienną wartością średnią przed przystąpieniem do analizy widmowej lub prognozowania należy zmiany niestacjonarne częściowo wyeliminować, poprzez odjęcie modelu deterministycznego lub obliczenie jedno lub wielokrotnych różnic szeregu czasowego.
ZWIĄZKI POMIĘDZY DWOMA
PROCESAMI STOCHATYCZNYMI
Związki pomiędzy dwoma procesami stochastycznymi
oraz
mogą
być określane przy pomocy następujących funkcji:
W praktyce dysponujemy skończonym czasem
trwania procesu ciągłego zatem estymator kowariancji wzajemnej dla przebiegów
ciągłych zdefiniowany jest następująco:
Obciążone i nieobciążone estymatory kowariancji
wzajemnej wyznaczane są analogicznie do estymatorów autokowariancji opisanych
wcześniej.
- korelacja wzajemna (cross-korelacja)
Korelacja wzajemna powstaje poprzez podzielenie kowariancji wzajemnej przez odchylenia standardowe dwóch szeregów czasowych.
- kowariancja
kowariancja jest szczególnym przypadkiem
kowariancji wzajemnej dla przesunięcia czasowego
.
- współczynnik korelacji
Współczynnik korelacji jest szczególnym
przypadkiem funkcji korelacji wzajemnej dla przesunięcia czasowego
.
Współczynnik korelacji to liczba z przedziału <-1,1>. Zerowy współczynnik
korelacji pomiędzy dwoma procesami stochastycznymi oznacza , że są one
nieskorelowane, ale nie oznacza wcale, że poszczególne składowe częstotliwościowe
i
w
tych szeregach są ze sobą nieskorelowane.
Koherencja zmienia się w przedziale od
0 do 1. Jeżeli
dla wszystkich
częstotliwości
,
wówczas szeregi czasowe
i
są
niekoherentne lub niezależne, nieskorelowane dla wszystkich częstotliwości.
Jeżeli
dla wszystkich częstotliwości,
wówczas szeregi czasowe
i
są
całkowicie koherentne.
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
SZUM BIAŁY I KOLOROWY
Przyjmuje się, że element losowy reprezentowany jest przez szum biały
(w
praktyce może to oznaczać błąd pomiarowy), który ma:
- wartość oczekiwaną równą zero,
- funkcję autokowariancji równą zero dla niezerowego przesunięcia czasowego
,
gdzie
jest
deltą Diraca,
- widmo mocy stałe dla wszystkich częstotliwości
znormalizowanych
.
Widmo szumu nieskorelowanego jest stałe
co oznacza, że moc szumu białego rozłożona jest równomiernie w całym paśmie
częstotliwości. Dzięki tej własności przebieg przypomina światło białe.
Ze względu na to, że autokowariancja szumu białego dla zerowego przesunięcia
czasowego (wariancja) jest deltą Diraca szum
biały nie daje się wytworzyć w układach fizycznych. W związku z tym, że
szum biały jest niemożliwy do wygenerowania, to każdy fizyczny szum jest
kolorowy. Widmo mocy szumu kolorowego jest stałe dla przedziału częstotliwości
takiego,
że
.
II. FILTRY CYFROWE
Równanie filtru cyfrowego ma następującą
postać ogólną:
gdzie
jest
długością filtru,
jest
liczba biegunów filtru,
jest
sygnałem wyjściowym,
jest
sygnałem wejściowym.
Jeżeli
dla
wówczas
filtr nazywamy nierekursywnym, jeżeli
wówczas
filtr jest typu rekursywnego. Transformując
transformatą
(tzn.
$Y(z)
), w szczególnym przypadku gdy
transformata
sprowadza
się do transformaty Fouriera) obie strony równania filtru otrzymamy:
Stosunek transformaty sygnału na wejściu
do transformaty sygnału wyjściowego określa tzw. funkcję przenoszenia
filtru:
gdzie
nazywamy
odpowiedzią impulsową filtru.
Taka odpowiedź impulsowa spełnia następujące
równanie splotu:
Jeżeli
nie
może być zredukowane do wielomianu skończonego rzędu od
wówczas
ma
nieskończoną długość, i filtr taki nazywamy filtrem nieskończonej odpowiedzi
impulsowej (infinite impulse response - IIR). W przeciwnym wypadku filtr
jest filtrem skończonej odpowiedzi impulsowej (finite impulse response
- FIR). Każdy filtr nierekursywny jest filtrem
typu FIR. Filtry rekursywne charakteryzują się tym, że sygnał wyjściowy
wprowadzany jest zwrotnie na wejście. Do filtrów tego typu należą filtry:
Butterwortha, Czebyszewa oraz Eliptyczne pierwszego i drugiego rodzaju.
Do filtrów nierekursywnych zaliczyć można filtr Ormsby, McLellan'a, Gauss'a
oraz filtry projektowane z wykorzystaniem okien wagowych.
Projektowanie filtru IIR lub FIR polega na znalezieniu takiej odpowiedzi
impulsowej
, dla której
funkcja
zwana
odpowiedzią amplitudową lub funkcją wzmocnienia będzie miała określony
kształt. Odpowiedzią fazową filtru
jest
funkcja mówiąca o tym jak zostaje przekształcona faza z wejścia
do wyjścia filtru. Ze względu na kształt funkcji odpowiedzi amplitudowej filtru mamy następujące rodzaje filtrów:
- dolnoprzepustowy (low pass filter - LPF)
- górnoprzepustowy (high pass filter - HPF)
- środkowoprzepustowy (band pass
filter - BPF)
- środkowozaporowy (band stop filter -
BSF)
Idealny LPF (tzw filtr podstawowy) powinien
przenosić informację w paśmie częstotliwośći od
do
,
a usuwać informację w paśmie od
do
,
gdzie
nazywamy
częstotliwością obcięcia.
W praktyce okazuje się, że niemożliwa jest
realizacja filtrów o charakterystyce bliskiej filtrom idealnym. Chcąc otrzymać
filtr HPF lub BPF wystarczy wykonać pewne operacje na odpowiedzi impulsowej
filtru podstawowego, np. odpowiedź impulsową filtru środkowoprzepustowego
uzyskuje
się poprzez obliczenie różnicy pomiędzy dwiema odpowiedziami impulsowymi
filtru LPF dla dwóch różnych częstotliwości obcięcia
,
.
Filtr Gaussa
LPF Gaussa jest niesymetrycznym i nierekursywnym
filtrem typu FIR o odpowiedzi impulsowej będącej funkcją Gaussa:
gdzie
jest
parametrem zmieniającym nachylenie zboczy funkcji,
jest długością okna wygładzającego oraz
.
Długość okna wygładzającego wyrażona jest przez tzw FWHM (full width at
half maximum), co oznacza odległość pomiędzy dwoma punktami na funkcji
Gaussa znajdującymi się na przeciwległych jej zboczach w połowie jej wysokości.
Sygnał wygładzony obliczany jest wzorem:
Filtr Ormsby
Filtr Ormsby (1961) jest symetrycznym filtrem nierekursywnym typu FIR
realizowanym równaniem:
gdzie
jest
odpowiedzią impulsową LPF Ormsby,
-
jest
częstotliwością obcięcia,
-
jest częstotliwością końcową (roll off termination frequency),
jest
długością filtru.
- środkowoprzepustowy filtr Ormsby
Odpowiedź impulsowa środkowoprzepustwego filtru Ormsby określona jest następującym wzorem:
gdzie
jest to tzw. pasmo
przenoszenia (transition band) określający
nachylenie zbocza funkcji przenoszenia,
jest
połową sumy pasma przepuszczania oraz pasma przenoszenia.
Funkcja przenoszenia środkowoprzepustowego
filtru Ormsby moduł transformaty Fouriera funkcji
jest
trapezem, którego górna podstawa jest równa
,
natomiast dolna podstawa jest równa
.
Dokładność tego filtru zależy od nachylenia zbocza filtru
.
Zmniejszanie parametru
powoduje
zmniejszenie dokładności filtru. Jednakową wartość błędu filtracji uzyskuje
się dla warunku
.
Filtr Butterworth'a
Kwadrat odpowiedzi amplitudowej LPF Butterwortha
wyrażony jest wzorem:
gdzie,
jest
częstotliwością obcięcia, dla której
,
jest
liczbą biegunów filtru.
Wstawiając w miejsce ilorazu
w
powyższym równaniu ilorazy
lub
otrzymujemy
odpowiednio sinusowy lub tangensowy filtr Butterworth’a. Zaprojektowanie
filtru Butterworth’a polega na znalezieniu
takich współczynników
oraz
w
równaniu ogólnym filtru cyfrowego, które spełniają powyższe równanie odpowiedzi
amplitudowej filtru. Większość filtrów rekursywnych w tym filtr Butterwortha
są filtrami niesymetrycznymi i fazy oscylacji ulegają zniekształceniu z
wejścia do wyjścia filtru. Jeżeli jednak przefiltrujemy dane w kierunku
rosnącego czasu
, a następnie
w kierunku przeciwnym, wówczas powstałe po pierwszej filtracji przesunięcie
fazowe zostanie wyeliminowane. Filtracja w dwóch kierunkach powoduje jednak
podniesienie funkcji przenoszenia do kwadratu.
Filtr Transforamty Fouriera
Filtr Transformaty Fouriera przedstawia
się następującym równaniem:
gdzie
jest operatorem
dyskretnej transformaty Fouriera
,
natomiast
jest funkcja
przenoszenia filtru, w której
jest
częstotliwością obcięcia filtru,
jest
argumentem funkcji przenoszenia.
- środkowoprzepustowy filtr Transformaty Fouriera
Najprostszym środkowoprzepustowym filtrem transformaty Fouriera jest
taki filtr, w którym funkcja przenoszenia jest funkcją prostokątną:
gdzie
jest okresem
centralnym wyznaczanej oscylacji,
jest
argumentem funkcji,
jest
połową szerokości pasma przenoszenia filtru.
Dokładność wyznaczenia oscylacje z szeregu
czasowego można poprawić jeżeli zamiast funkcji prostokątnej wykorzystamy
funkcję paraboliczą (Popiński and Kosek 1995a,b):
Wartość parametru
ma
bardzo duże znaczenie na wynik otrzymanej oscylacji. Czym większy parametr
tym
oscylacja ma szybciej zmieniająca się amplitudę w czasie. Parametr
musi
spełniać następujący warunek:
gdzie
jest
liczbą danych.
W celu skrócenia czasu obliczeń Dyskretna Transformata Fouriera, zastąpiona została Szybką Transformatą Fouriera - FFT mixed-radix Fast Fourier Transformation) (Singleton 1969). Zastosowanie transformaty Singleton'a znacznie zwiększa użyteczność FFT do obliczeń, ze względu na to, że liczba danych niekoniecznie musi być równa całkowitej potędze liczby 2, tak jak w przypadku klasycznej metody FFT (Cooley i Tuckey 1965), a jest iloczynem całkowitych potęg liczb 2, 3 i 5. Pozwala na lepsze dopasowanie długości danych do obliczeń. Najważniejszą jednak zaletą FFT jest znaczne przyspieszenie obliczeń, szczególnie wtedy, gdy wzrasta liczba danych.
III. ANALIZY WIDMOWE
Wprowadzenie
Zadaniem analizy widmowej jest wyznaczenie
widma mocy szeregu czasowego w funkcji częstotliwości lub okresu. Wartości
tego widma dla określonych częstotliwości są na ogół proporcjonalne do
odpowiadających im kwadratów amplitud oscylacji. Przejście z dziedziny
czasu do dziedziny częstotliwości realizowane jest przy pomocy Transformaty
Fouriera, która leży u podstaw każdej
analizy widmowej, gdyż przenosi w sposób pełny informacje z dziedziny czasu
do dziedziny częstotliwości. Podstawą analiz widmowych jest transformata
Fouriera pozwalająca na przeniesienie informacji z dziedziny czasu
do
dziedziny częstotliwości
:
Zadaniem analizy spektralnej jest wyznaczenie
estymatora pewnej dodatniej funkcji częstotliwości, której maksima pokazują
energię oscylacji znajdujących się w danym procesie stochastycznym. Najprostszą
formą analizy spektralnej jest obliczenie kwadratu modułu transformaty
Fouriera
, jednak
taki estymator widma mocy jest bardzo mało
dokładny. Czas obliczenia Dyskretnej Transformaty Fouriera (Discrete Fourier
Transform - DFT) można skrócić poprzez zastosowanie Szybkiej Transformaty
Fouriera (Fast Fourier Transform - FFT), lub Szybkiej Transformaty Fouriera
Singleton’a.
Inną wielkością, którą uzyskuje się przy pomocy FFT jest kąt fazowy. Jeżeli sygnał reprezentowany przez sumę poszczególnych składowych okresowych o różnych częstotliwościach, kąt fazowy pokazuje relacje pomiędzy różnymi składowymi częstotliwościowymi. Widmo fazowe mówi jaki jest kąt fazowy odpowiadający każdej częstotliwości.
Analiza Harmoniczna
Analiza harmonicza zakłada, że proces
stochastyczny
da się przedstawić
podobnie jak funkcja ciągła
szeregiem
harmonik (funkcji okresowych) czyli sumą trygonometryczną:
gdzie
,
i
nazywamy
współczynnikami Fouriera.
Gdy założymy, że szereg czasowy
jest
równoodległy, wówczas dla wystarczająco dużej liczby danych
estymatory
współczynników Fouriera będą określone następującymi wzorami:
,
Przedstawiona metoda polega na korelowaniu
różnych fal testowych o częstotliwościach
z
układem, w którym występują częstotliwości
.
Wystąpienie korelacji oznacza zdudnianie fali testowej o częstotliwości
z
falą o częstotliwości
,
która znajduje się w danym procesie stochastycznym.
nazywamy periodogramem,
natomiast
nazywamy estymatorem
widma fazowego.
Metoda Blackman'a-Tuckey'a
Własność funkcji autokowariancji wykorzystali
Blackaman i Tuckey (1965 ) wprowadzając do analiz szeregów czasowych metodę
analizy widmowej będącej Transformatą Fouriera estymatora funkcji autokowariancji
procesu, który wyrażony jest następującym wzorem:
gdzie
jest oknem wagowym,
jest
liczbą punktów procesu,
jest
przesunięciem czasowym.
Własność funkcji autokowariancji najlepiej
pokazać na przykładzie. Jeżeli założymy że:
wówczas
dla liczby danych procesu
jego
funkcja autokowariancji wyrażona jest następującym wzorem:
gdzie
,
,
oznaczają
amplitudy, częstotliwości i fazy.
Łatwo zauważyć, że funkcja autokowariancji
ma
tą samą charakterystykę częstotliwościową co sam proces stochastyczny
,
z którego jest wyznaczona. W praktyce liczbą
danych jest ograniczona czasem obserwacji
i
obcięta funkcja autokowariancji zawiera dodatkowe harmoniki, które nie
występują w szeregu czasowym.
Wzór:
stanowi
podstawę metody Blackmana-Tuckeya. Ze względu na to, że funkcja autokowariancji
jest parzysta to wyznaczony tą metodą estymator widma mocy
będzie miał mniejszy błąd kwadrat modułu
Transformaty Fouriera ze względu na to, że nieobciążony estymator autokowariancji
nie
zawiera elementu losowego tak jak proces stochastyczny
.
W metodzie Blackman'a-Tuckey'a występują dwa
efekty powodujące, że estymator widma mocy nie jest równy teoretycznemu
widmu mocy. Po pierwsze estymator funkcji autokowariancji
wyznaczony
z ograniczonego czasem trwania procesu
generuje
pewne harmoniki nie występujące w tym procesie. Po drugie występuje tzw.
efekt przecieku widma mocy objawiający się tym, że w estymatorze widma
mocy pojawiają się dodatkowe maksima nie odpowiadające rzeczywistym oscylacjom.
Maksima te jednak zanikają wraz ze wzrostem ilości punktów
procesu.
Pomnożenie nieobciążonego estymatora autokowariancji przez okno wagowe
częściowo
eliminuje efekt przecieku widma mocy, natomiast pierwszy efekt można jedynie
wyeliminować poprzez wydłużenie procesu jeżeli jest to możliwe.
MESA - Maximum Entropy Spectral Analysis
MESA jest oparta na kryterium doboru takiego estymatora widma mocy,
który maksymalizuje następującą funkcję:
Została opublikowana przez Burga w 1967
roku. Jest to również transformata Fouriera autokowariancji procesu, który
wyrażony jest przez współczynniki autoregresji. W latach 70-tych zauważono,
że metoda MESA jest widmem procesu autoregresji, który opisany jest równaniem:
Wartość procesu autoregresji w momencie
czasu
wyrażona
jest poprzez jego wartości poprzednie pomnożone przez współczynniki autoregresji.
Czym rząd autoregresji jest wyższy tym proces autoregresji dokładniej aproksymuje
niższe częstotliwości procesu stochastycznego. Widmo mocy tej metody wyprowadza
się jako transformatę Fouriera funkcji
autokowariancji
, wyrażonej
przez równanie procesu autoregresji
Widmo to ma postać:
Estymator widma mocy metody MESA wyrażony
jest następującym równaniem:
gdzie
jest estymatorem
wariancji szumu.
Gdy rząd autoregresji
widmo
mocy metody MESA jest rosnące albo malejące, gdy
widmo
mocy może mieć najwyżej jedno maksimum. Aby liczba maksimów widma mocy
odpowiadała rzeczywistym oscylacjom procesu stochastycznego wartość rzędu
autoregresji
powinna być
optymalna. Jest wiele różnych sposobów określenia rzędu autoregresji. W
szczególności dla zaszumionych procesów stochastycznych (zbliżonych charakterem
do szumu kolorowego) polecane jest zmodyfikowane kryterium Rovelli - Vulpiani
(1983):
- MESA estymatora autokowariancji
Estymator funkcji autokowariancji można
potraktować jak proces stochastyczny
,
ze względu na to, że posiada on taką samą charakterystykę częstotliwościową
co ten proces. Z równania funkcji autokowariancji wyrażonej przez współczynniki
autoregresji widać, że w odróżnieniu od procesu autoregresji pozbawiona
jest elementu losowego
,
zatem wyznaczone widmo mocy metodą MESA estymatora autokowariancji będzie
dokładniejsze niż widmo mocy metody MESA procesu
.
Pamiętać należy jednak, że wysokości maksimów widma mocy będą wtedy proporcjonalne
do amplitud oscylacji nie w drugiej a w czwartej potędze.
Analiza spektralna metodą filtru środkowoprzepustowego
Definicji widma mocy jest następująca:
gdzie
jest
wyjściem idealnego filtru środkowoprzepustowego, do którego wejściem jest
proces stochastyczny
o
zerowej wartości oczekiwanej. W powyższym równaniu
jest
splotem wejścia filtru
z
odpowiedzią impulsową idealnego filtru środkowoprzepustowego
:
Filtr ten przepuszcza częstotliwości w
wąskim paśmie
dla
.
W praktyce odpowiedź impulsowa filtru środkowoprzepustowego
przepuszcza częstotliwości w paśmie, którego najmniejsza szerokość wynosi
.
- zastosowanie środkowoprzepustowegi filtru Ormsby
Jeżeli do analizy widmowej zastosowany zostanie środkowoprzepustowy filtru Ormsby, wówczas estymator widma mocy wyrażony będzie następującym równaniem:
gdzie
,
jest pasmem przenoszenia (transition band) określającym nachylenie zbocza
funkcji przenoszenia,
jest
połową sumy pasma przepuszczania oraz pasma przenoszenia,
jest
długością filtru.
- zastosowanie środkowoprzepustowego filtru transformaty Fouriera
Oscylację o częstotliwości
można
wyznaczyć filtrem Transformaty Fouriera następującym wzorem:
gdzie
jest operatorem
transformaty Fouriera,
jest
funkcją przenoszenia filtru określoną równaniem:
jest
argumentem częstotliwości,
jest
połową szerokości pasma przenoszenia.
W celu uzyskania oscylacji o mniejszym
błędzie może być zastosowana paraboliczna funkcja przenoszenia:
Wzór na estymator widma mocy ma postać
:
gdzie
jest
liczbą punktów, które należy pominąć na końcach wyznaczanej oscylacji
ze
względu na błędy końców.
Proces stochastyczny
jako
wejście do filtru transformaty Fouriera może być również procesem zespolonym.
Wówczas wyznaczona tym filtrem oscylacja jest również procesem zespolonym.
IV. PROGNOZOWANIE
Wprowadzenie
W prognozowaniu numerycznym głównym problemem jest wyznaczenie przewidywanej wartości szeregu czasowego poza przedziałem czasowym, w którym szereg ten jest określony. Wpływ zjawisk losowych powoduje, że nie znamy prawdziwej wartości szeregu czasowego poza punktami (tj. także pomiędzy punktami), w których jest on określony. Najprostszą metodą prognozowania numerycznego jest wyznaczenie składowej deterministycznej szeregu czasowego, w postaci pewnej funkcji określonej wzorem matematycznym. Funkcja ta wyznaczona poprzez dopasowanie metodą np. najmniejszych kwadratów modelu matematycznego do punktów szeregu czasowego i ekstrapolowana poza przedział czasowy, w którym szereg ten jest określony nazywana jest prognozą. Oprócz deterministycznych metod prognozowania znane są także metody stochastyczne takie jak metoda autoregresji (Box i Jenkins 1970) oraz filtr Kalmana (Gelb 1988; Morabito i in. 1988). W niniejszym opracowaniu przedstawione są tylko stochastyczne metody prognozowania. Estymator prognozy stacjonarnych szeregów czasowych obliczony dowolną metodą prognozowania jest tym dokładniej wyznaczony im większa jest liczba danych. Jeżeli szereg czasowy jest niestacjonarny, wówczas prognozowanie staje się trudne do rozwiązania. Zastosowana do niestacjonarnych danych metoda prognozy nie jest w stanie przewidzieć zmian niestacjonarnych. Jeżeli założymy, że zmiany niestacjonarne charakterystyki częstotliwościowej lub wariancji zmieniają się nieznacznie wokół ich średnich wartości, wówczas dla dostatecznie długiego szeregu czasowego eliminują się. W takim przypadku zastosowana liniowa metoda prognozy jest w stanie przewidywać część szeregu czasowego, która nie jest niestacjonarna. Stochastyczne metody prognozowania autoregresji i autokowariancji wykorzystują informację o korelacji pomiędzy oddalonymi od siebie o różne odstępy czasowe punktami szeregu czasowego (autokorelacji) i ekstrapolują tą informację poza przedział czasowy danych w postaci szeregu prognozy. W obliczonej prognozie punkty powinny być skorelowane ze sobą w podobny sposób jak punkty w samym szeregu czasowym. Charakterystyka częstotliwościowa prognozy powinna być taka sama jak szeregu czasowego, z którego została wyznaczona.
metoda autoregresji
Niektóre geofizyczne szeregi czasowe zawierają oscylacje, których amplitudy mogą być pobudzane w sposób losowy i wówczas najlepszą metodą ich reprezntacji jest model autoregresji (Akaike, 1969; Ulrych and Bishop, 1975; Ulrych and Ooe, 1979):
Stacjonarny proces stochastyczny
o
wartości oczekiwanej równej zero, określony w równoodległych momentach
czasu można modelować za pomoca procesu autoregresji rzędu M następującym
równaniem (Box i Jenkins 1970; Priestley 1981):
gdzie:
są
współczynnikami autoregresji,
jest szumem białym.
Dla procesu autoregresji jego wartości
w następnych momentach czasu
,
,...,
również
spełniają równanie autoregresji:
.
W równaniach tych przyszłe wartości szumu
białego
,
,...,
są
nieznane, zatem należy je zastąpić wartością oczekiwana szumu białego
.
Poza tym, prawdziwe wartości współczynników autoregresji
oraz
przyszłych wartości procesu
muszą
być zastąpione ich estymatorami
oraz
.
Znając rząd autoregresji
oraz
estymatory współczynników autoregresji
dla
,
prognozy obliczane są ze wzorów:
.
Prognozę w następnych momentach czasu możemy wyznaczyć tylko wówczas, gdy znane są jej wartość w punktach poprzednich. Problemem w prognozowaniu metodą autoregresji jest właściwe określenie rzędu autoregresji oraz estymatorów współczynników autoregresji.
- rząd autoregresji
Rząd autoregresji można wyznaczyć kryterium
Akaike (AIC) - Akaike's Information Criterion (Akaike, 1974):
gdzie
jest
estymatorem wariancji residuów otrzymanym metodą największego prawdopodobieństwa
opartą na modelu
parametrowym,
jest
wariancją residuów.
Innym sposobem wyznaczenia rzędu autoregresji
jest kryterium końcowego błędu prognozy FPE (final prediction error):
,
Wymienione wyżej kryteria odnoszą się
do procesów stochastycznych, w których amplituda sygnału jest na ogół większa
od odchylenia standardowego szumu. Do wyznaczenia rzędu autoregresji
zaszumionych
szeregów czasowych może być zastosowane zmodyfikowane kryterium
Rovelli-Vulpiani (Rovelli i Vulpiani 1983):
gdzie
jest estymatorem
autokowariancji, w którym
jest
oknem wagowym,
jest
liczbą danych,
jest przesunięciem
czasowym.
W teorii szeregów czasowych okno wagowe Kaiser'a posiada pewne szczególne właściwości, a mianowicie w sposób najbardziej optymalny zmniejsza obciążenie wyznaczanych estymatorów, które spowodowane jest skończonym czasem obserwacji szeregów czasowych. Okno wagowe Kaisera wyprowadzone zostało z warunku minimum sumy kwadratów różnicy pomiędzy rzeczywistym widmem mocy nieskończonego szeregu czasowego, a jego estymatorem, który obliczany jest ze skończonego szeregu czasowego (Tretter 1976).
Rząd autoregresji
,
jak również estymatory współczynników autoregresji z obciążonego estymatora
autokowariancji szeregu czasowego, który jest iloczynem estymatora nieobciążonego
pomnożonego przez współczynniki okna wagowego Kaisera
(Kosek
1992; Tretter 1976).
Niektóre geofizyczne szeregi czasowe mogą być przykładem szeregów czasowych przypominających swoim charakterem szum kolorowy. W tym przypadku zastosowanie innych kryteriów wyboru optymalnego rzędu autoregresji, np kryteriów FPE (Final Prediction Error) i AIC (Akaike's Information Criterion) (Akaike 1969; Barrodale i Erickson 1980) daje bardzo zaniżoną wartość samego rzędu autoregresji, co prowadzi do niewystarczająco dokładnej reprezentacji niskich częstotliwości przez równanie autoregresji. Rząd autoregresji w decydujący sposób wpływa na charakterystykę częstotliwościową obliczanej prognozy, a długookresowa charakterystyka częstotliwościowa jest tym lepiej wyznaczona im rząd autoregresji jest wyższy. Na ogół zmodyfikowane kryterium Rovelli-Vulpiani wyznacza wyższy rząd autoregresji niż inne kryteria. Wzrost rzędu autoregresji powoduje z drugiej strony wprowadzenie dodatkowych korelacji w zakresie wysokich częstotliwości do szeregu opisanego modelem autoregresji.
- estymatory współczynników autoregresji
Jeżeli rząd autoregresji procesu stochastycznego
jest ustalony estymatory współczynniki autoregresji można obliczyć z równań
Yule'a-Walker'a (Box i Jenkins 1970).
Prognoza autokowariancyjna szeregów czasowych jednowymiarowych
Prognoza metodą autokowariancyjną prognozująca
stacjonarne szeregi czasowe jednowymiarowe zaprezentowana została po raz
pierwszy w tomie prac konferencyjnych (Kosek 1993a; Kosek i in. 1995b).
Jak wskazuje sama nazwa tej metody, główną rolę w niej odgrywa funkcja
autokowariancji, a właściwie jej estymator. Funkcja autokowariacji nieskończonego
stacjonarnego ciągłego szeregu czasowego, którego wartość oczekiwana jest
równa zero jest określona następującym równaniem (Otnes i Enochson 1972):
gdzie
jest przesunięciem czasowym.
Funkcja ta zawiera pełną informację o charakterystyce
częstotliwościowej szeregu czasowego, z którego jest wyznaczona. Dla przykładu
jeżeli szereg czasowy
jest
sumą sinusoid o stałych amplitudach, okresach i fazach oraz szumu białego,
wówczas jego funkcja autokowariancji dla
jest
także sumą sinusoid o tych samych okresach, fazach równych zero i amplitudach,
które są kwadratami amplitud sinusoid szeregu czasowego, oraz dla
jest
wariancją tego procesu. W praktyce szeregi czasowe są skończone i określone
jedynie w punktach dyskretnych i zamiast funkcji autokowariancji posługujemy
się jej estymatorem. Obciążony estymator autokowariancji jest określony
następującym wzorem:
gdzie xt jest stacjonarnym, szeregiem
czasowym określonym w równoodległych momentach czasu, o wartości oczekiwanej
równej zero,
jest liczbą
punktów szeregu czasowego.
Jeżeli w szeregu czasowym zwiększymy liczbę
punktów o jeden dodając nieznany 'estymator pierwszego punktu prognozy'
wówczas
uzyskany nowy nieobciążony estymator autokowariancji określony jest równaniem:
Nieobciążony estymator autokowariancji
pomnożnony przez dowolne okno wagowe, różne od okna prostokątnego, nazywamy
obciążonym estymatorem autokowariancji. Ideą autokowariancyjnej metody
prognozowania jest znalezienie takiego 'estymatora pierwszego punktu
prognozy'
,
aby suma kwadratów różnic obciążonych estymatorów autokowariancji szeregu
czasowego i autokowariancji szeregu czasowego z dodanym pierwszym punktem
prognozy była jak najmniejsza, co można zapisać następująco (Kosek 1993a,
1997):
gdzie
jest oknem wagowym
dla
,
jest
oknem wagowym dla
.
Po przyrównaniu pierwszej pochodnej powyższej
funkcji do zera:
otrzymujemy wielomian trzeciego stopnia:
, w którym nieznany 'estymator
pierwszego punktu prognozy'
jest
rzeczywistym pierwiastkiem tego równania:
a współczynniki
są
określone następującymi wzorami:
Równania trzeciego stopnia ma jeden pierwiastek rzeczywisty i dwa zespolone
gdy okno wagowe
spełnia
następujący warunek:
W analizie szeregów czasowych wszystkie
okna wagowe spełniają powyższy warunek i wówczas równanie 3-go stopnia
ma jeden pierwiastek rzeczywisty i dwa zespolone (Kosek 1993a).
Badania modelowe na szeregach czasowych
będących kombinacją liniową sinusoid o stałych okresach, amplitudach i
fazach wykazały, że wyznaczona bezwzględna wartość 'estymatora pierwszego
punktu prognozy' jest zwykle mniejsza od prawdziwej bezwzględnej wartości
modelowego szeregu czasowego. Zbiór wartości obliczonych 'estymatorów pierwszych
punktów prognozy' zachowuje natomiast charakterystykę częstotliwościową
szeregu czasowego. Zadanie to rozwiązane jest w sposób polegający na tym,
że 'pierwszy punkt prognozy'
jest
iloczynem 'estymatora pierwszego punktu prognozy'
oraz
pewnego stałego współczynnika
,
który jest średnią obliczoną metodą estymacji szorstkiej (ang. robust)
(Press i in. 1992). Wartość tego współczynnika uzyskuje się poprzez porównanie
'estymatorów pierwszych punktów prognozy' z rzeczywistymi wartościami szeregu
czasowego dla ostatnich
punktów
szeregu czasowego.
Średnia 'robust'
(Priestley
1981), obliczana jest z następującego warunku:
gdzie parametry
dla
są
obliczane na podstawie znanych punktów szeregu czasowego
oraz
ich pierwszych estymatorów prognozy
,
gdzie przyjete zostało apriori
(Kosek
i in. 1995b). Jeżeli dla ![]()
wówczas
parametr musi być zastąpiony poprzednim
.
Po obliczeniu pierwszego punktu prognozy
,
jest on dołączany do szeregu czasowego, co umożliwia wyznaczenie następnego
punktu prognozy
itd.
Jeżeli szereg czasowy ma niestacjonarną
wartość oczekiwaną, np. pojawia się w nim trend liniowy, wówczas można
obliczyć jego jednokrotnie różnice
dla
a następnie obliczyć prognozę takich różnic.
Prognoza szeregu czasowego jest wówczas zsumowaną jednokrotnie prognozą
różnicy.
W tej metodzie prognozowania każde okno
wagowe może zostać zastosowane jednak najbardziej odpowiednie jest okno
wagowe Keiser'a ze wzgledu na to, że błąd estymatora autokowariancji wynikający
ze skończonej długości szeregu czasowego jest najmniejszy. Zastosowanie
tego okna wagowego wydłuża jednak czas realizacji algorytmu metody autokowariancyjnej
ze względu na występujące funkcje Bessela
w równaniu tego okna wagowego. Czas realizacji algorytmu metody autokowariancyjnej
można jednak skrócić, po pierwsze, poprzez zastosowanie metody rekurencyjnej
do obliczania następnego nieobciążonego estymatora autokowariancji, gdy
znany jest poprzedni:
dla
.
Po drugie, czas realizacji tego algorytmu
można skrócić jeszcze bardziej, jednak z nieznacznym pogorszeniem dokładności
prognozy, poprzez zastosowanie okna wagowego Bartletta (Tretter
1976):
Zastosowanie tego okna wagowego spowoduje
znaczne uproszczenie wzorów na współczynniki
w
równaniu 3-go stopnia:
Następnym krokiem do zredukowania czasu
obliczeń algorytmu metody autokowariancyjnej jest kolejne uproszczenie
wzorów, które wyprowadzone zostały z nieco zmienionego warunku, w którym
zaniedbane zostały estymatory wariancji:
Jeżeli nastepnie zastosowane jest okno
wagowe Bartlett'a 'estymator pierwszego punktu prognozy'
wyrażony
jest następującym równaniem:
gdzie
oknem wagowym
Bartletta.
Z powyższego wzoru możemy odczytać, że 'estymator pierwszego punktu prognozy' jest dyskretną realizacją splotu szeregu czasowego z jego obciążonym estymatorem autokowariancji, który podzielony jest przez (N-2) krotny estymator wariancji tego szeregu. We wzorze tym odnajdujemy także sens autokowariancyjnej metody prognozowania, w której najważniejszą rolę odgrywa estymator autokowariancji szeregu czasowego.
Prognoza autokowariancyjna szeregów czasowych dwuwymiarowych
Dwuwymiarowe stacjonarne szeregi czasowe
można prognozować metodą autokowariancyjną podobnie jak jednowymiarowe
stacjonarne szeregi czasowe. Nieobciążony estymator autokowariancji dwuwymiarowego
szeregu zasowego
,
dla
określony jest następującym równaniem:
, dla
gdzie
jest
liczbą danych.
Estymator autokowariancji zespolonego szeregu
czasowego może być wyrażony poprzez estymatory autokowariancji jego części
urojonej i rzeczywistej:
Po przedłużeniu szeregu czasowego o pierwszy
punk prognozy
jego estymator
autokowariancji jest określony następującym równaniem rekurencyjnym:
W prognozie autokowariancyjnej suma kwadratów
modułów różnicy pomiędzy estymatorem autokowariancji szeregu czasowego
oraz przedłużonego szeregu czasowego powinna być najmniejsza (Kosek 1993):
gdzie:
i
są
oknami wagowym o długościach
oraz
.
Powyższa funkcja osiąga minimum jeżeli:
i
,
Rozwiązaniem tych równań jest składowa
rzeczywista i urojona pierwszego punktu prognozy:
gdzie
,
,
,
.
Po przedłużeniu szeregu czasowego o pierwszy punktu prognozy następny punkt prognozy może być obliczony w sposób rekurencyjny.
V. CZASOWO-CZĘSTOTLIWOŚCIOWE ANALIZY WIDMOWE
Dla szeregów czasowych niestacjonarnych z usuniętą niestacjonarnością wartości średniej wynik analizy widmowej jest zależny od przedziału czasowego danych, którego ta analiza dotyczy. Przedział czasowy danych można zawsze skrócić badając widmo mocy w różnych zakresach przesuwanych w czasie. Ograniczenie liczby danych szeregu czasowego wiąże się jednak ze zmniejszeniem dokładności wyznaczanego estymatora widma mocy, ze względu na efekt przecieku widma mocy. Każda analiza widmowa daje jedynie informacje o średnich amplitudach częstotliwości znajdujących się w procesie stochastycznym. Często jednak mamy do czynienia ze zmieniającymi się w czasie amplitudami oscylacji. Informacja o czasowych zmianach spektrum może być uzyskana poprzez zastosowanie analizy widmowej nie do całości danych a do przedziału danych, który porusza się po całej długości danych. Zawężanie tego przedziału czasowego poprawia czasową rozdzielczość otrzymanego 'czasowo zmiennego widma mocy' jednak pogarsza jego dokładność w częstotliwości ze względu na tzw. efekt przecieku widma mocy. W związku z tym zaczęto poszukiwać innych metod czasowo-częstotliwościowych analiz widmowych. W przypadku analiz widmowych problem niestacjonarności może być rozwiązany poprzez jednoczesne określenie charakterystyki czasowej i częstotliwościowej szeregu czasowego. Wcześniejsze analizy czasowo częstotliwościowe wykorzystywały Transformatę Gabora (Gabor 1946), będącą modyfikacją Transormaty Fouriera polegającą na wprowadzeniu pod znak całki okna Gaussa przesuwającego się w dziedzinie czasu. Umożliwiało to określoną lokalizację czasową informacji o częstotliwości. Podstawową wadą tej Transformaty była niejednakowa rozdzielczość dla niskich i wysokich częstotliwości (Chui 1992). Na początku lat 80-tych wprowadzona została do analiz procesów nieliniowych Transformata Wavelet, której podstawowym zadaniem jest przejście do dziedziny czasowo-częstotliwościowej (Chui 1992; Gasquet i Witomski 1990; Grossman i Morlet 1984, 1988; Goupillaud i in. 1984). Cechą charakterystyczną Transformaty Wavelet jest jednakowa rozdzielczość informacji w dziedzinie czasu i częstotliwości (Farge i Rabreau 1988). Moduł Transformaty Wavelet szeregu czasowego pokazuje zmiany energii oscylacji w funkcji czasu i częstotliwości. Czasowo zmienna analiza spektralna wykorzystuje dowolny filtr środkowoprzepustowy. W odróżnieniu od transformaty Wavelet skonstrułowana ona jest w taki sposób, aby można było sterować rozdzielczością czasową i częstotliwościową widma mocy przy czym wzrost rozdzielczości czasowej zmniejsza rozdzielczość częstotliwościową i odwrotnie.
Transformata Gabor'a
W 1946 roku Gabor zmodyfikował transformatę
Fouriera wprowadzając lokalizacyjne okno (funkcję) Gaussa
,
w której parametr
przesuwał
okno w dziedzinie czasu. Okno to okazało się jednak nieefektywne do analizowania
sygnałów jednocześnie o wysokich i niskich częstotliwościach. Transformata
Gabora jest pewną modyfikacją Transformaty Fouriera, do której wprowadzane
jest pod znak całki okno wagowe Gauss'a i wyrażona jest następującym wzorem:
gdzie
jest oknem (funkcją) Gauss'a, w którym
parametr
.
Zadaniem transformaty Gabora jest lokalizacja informacji z Transformaty
Fouriera z szeregu czasowego
wokół
czasu
.
Transformata Wavelet
W 1983 roku Morlet analizując dane sejsmiczne
wprowadził dodatkowy parametr skali
w
transformacie Gabora do lokalizacyjnego okna Gaussa, i tak powstała transformata
Wavelet. Funkcja analizująca
dzięki
parametrowi skali rozszerza się lub zawęża w zależności od częstotliwości.
Transformata Wavelet jest funkcją dwuwymiarową i jest splotem sygnału
z
tą funkcją analizującą:
gdzie
jest parametrem
translacji lub przestrzeni,
jest
parametrem skali lub dylatacji.
jest funkcją analizującą Morleta, w której
parametr
.
Transformata Wavelet (dyskretnego procesu stochastycznego
charakteryzuje
się tym, że funkcja
osiąga
optymalny kompromis pomiędzy rozdzielczością skali i przestrzeni (częstotliwości
i czasu). Taka informacja o czasie i częstotliwości zawarta w
pozwala
badać widmo procesu stochastycznego
lokalnie
w czasie. Warto dodać, że transformata Fouriera daje bardzo dobrą rozdzielczość
skali (częstotliwości), podczas gdy nie daje żadnej rozdzielczości przestrzeni
(czasu). Moduł transformaty Wavelet może być również wyznaczony dla zespolonego
procesu stochastycznego. W tym przypadku można badać energię prawoskrętnych
(
) i lewoskrętnych (
)
oscylacji. W przypadku gdy
parametr
skali (dylatacji) może być interpretowany jako okres badanej oscylacji
procesu stochastycznego
.
Transformata Wavelet ma obecnie szerokie inne zastosowania
niż analiza czasowo - częstotliwościowa sygnałów, między innymi:
- Podobnie jak w przypadku rozwijania funkcji
w
szereg funkcji trygonometrycznych można sygnał przedstawić szeregami funkcji
wavelet. Funkcje te odpowiadają pewnym przedziałom częstotliwości, i mogą
być łatwo analizowane. Usuwanie defektów na starych płytach gramofonowych
może być przykładem zastosowania transformaty Wavelet, w której dokonuje
się pewnych operacji na funkcjach, które związane są z częstotliwościami
powodującymi te defekty.
- Badanie dźwięku i głosu. Poszczególne
sylaby czy litery odpowiadają pewnym charakterystycznym częstotliwościom,
które można przypisać określonym wartościom parametru skali
.
Dzięki temu łatwo jest rozpoznawać dźwięk lub mowę w postaci zdigitalizowanej,
po czym dowolnie go obrabiać np. w celu nawiązania dźwiękowego kontaktu
z komputerem.
- Zastosowanie przy kompresji danych.
Proces stochastyczny rozwinięty w szereg funkcji wavelet może być zapamiętany
w postaci współczynników rozwinięcia, które obcięte na pewnym poziomie
dokładności pozwalają poprzez odwrotną transformatę Wavelet odtworzyć go
z wymaganą dokładnością. Okazuje się że informacja przechowywana w
postaci współczynników rozwinięcia zajmuje około 20 razy mniej pamięci
dyskowej niż sam zbiór stanowiący proces stochastyczny.
- Zastosowanie w analizie obrazu do poprawiania
jego ostrości. Służy do tego transformata Haara, która należy również do
kategorii transformat Wavelet.
- Zastosowanie w astronomii do badania
spektrogramów światła gwiazd, co pomoga je klasyfikować. Poza tym gromady
gwiazd czy galaktyk udaje się grupować w pewne klasy odpowiadające pewnym
skalom (częstotliwościom). W ten sposób można badać rozmieszczenie galaktyk
w przestrzeni, a także rozmieszczenie obszarów próżni poprzez analizę minimów
w widmie mocy transformaty Wavelet.
Harmoniczna Transformata Wavelet
W analizie harmonicznej transformaty wavelet
(HWT) współczynniki transformaty wavelet dla sygnału
określone są następującym równaniem splotu (Newland, 1998):
,
gdzie
jest
zespoloną harmoniczną funkcją wavelet zlokalizowaną w dziedzinie częstotliwości
blisko częstotliwości
oraz
jest
parametrem translacji.
Transformata Fouriera obu stron równania
pozwala otrzymać obraz współczynnika transformaty wavelet w dziedzinie
częstotliwości:
W praktyce mnożymy dyskretną transformatę
Fouriera szeregu czasowego
przez
transformatę Fouriera analizującej funkcji wavlet, która ma kształt funkcji
prostokątnej. W dziedzinie częstotliwości transformata Fouriera analizującej
funkcji wavelet mnożona jest dodatkowo przez funkcję gaussa w celu otrzymania
lepszej rozdzielczości częstotliwościowej (Newland 1998). Następnie po
zastosowanniu odwrotnej transformaty Fouriera estymatory współczynników
transformaty wavelet określone są następującym równaniem:
,
gdzie
,
,
jest szeregiem czasowym,
,
is
znormalizowaną centralną częstotliwością funkcji prostokątnej pomnożonej
przez okno gaussa:
gdzie
jest połową szerokości pasma oraz
jest
parametrem wygładzającym.
Dla sygnału
o
liczbie
punktów otrzymujemy
estymatorów
współczynników transformaty wavelet, które aproksymują współczynniki określone
całką splotu.
Czasowo-częstotliwościowa analiza widmowa metodą filtru środkowoprzepustowego
Estymator widma mocy czasowo zmiennej analizy spektralnej z wykorzystaniem
filtru środkowoprzepustowego przedstawiony jest następującym równaniem:
gdzie
jest okresem centralnym
oscylacji
,
jest
liczbą punktów zawartą w całkowitej
wielokrotności
fali o okresie
,
jest
okresem próbkowania danych.
Oscylację
można
wyznaczyć dowolnym filtrem środkowoprzepustowym. W przypadku środkowoprzepustowego
filtru transformaty Fouriera oscylacje otrzymujemy wzorem:
w którym
jest paraboliczną funkcją przenoszenia,
jest
połową szerokości pasma funkcji przenoszenia w jednostkach częstotliwości
znormalizowanej (Popiński and Kosek 1995a,b).
Zmniejszanie parametru
zwiększa rozdzielczość częstotliwościową czasowo zmiennego widma mocy,
podczas gdy rozdzielczość czasowa zmniejsza się (Kosek and Kołaczek 1996).
Rozdzielczość czasowa spektrum zależy także od liczby fal
i
zmniejsza się jeżeli liczba ta jest zwiększana (Kosek 1995a). Parametru
musi spełniać następującą nierówność:
.
Rozdzielenie spektrum na lewoskrętną i prawoskrętna część możliwe jest
tylko w przypadku dwuwymiarowych szeregów czasowych. Jeżeli zastosowany
zostanie filtr Transformaty Fouriera oraz badany proces stochastyczny jest
procesem dwuwymiarowym wówczas możliwe jest badanie zmian energii lewoskrętnych
(
) i prawoskrętnych
oscylacji.
Pierwiastek kwadratowy czasowo zmiennego spektrum szeregu czasowego dwuwymiarowego 'nazywany widmem amplitudowym' jest proporcjonalny do chwilowych amplitud lewoskrętnych i prawoskrętnych oscylacji. Obliczenia estymatora widma mocy można przyspieszyć poprzez zastosowanie FFT Singletona.
Czasowo-częstotliwościowa koherencja
metodą filtru środkowoprzepustowego
Czasowo zmienna koherencja pomiędzy dwoma
szeregami czasowymi
i
mówi
o zmianach współczynnika korelacji pomiędzy nimi w funkcji częstotliwości
i określona jest następującym wzorem:
gdzie
,
,
.
LITERATURA
Akaike H, (1969) Power Spectral Estimation Through the Autoregressive Model Fitting. Ann. Inst. Stat. Math. 21 pp 69-72.
Akaike H. (1974). A New Look at the Statistical Model Identification. IEEE Trans. Automatic Control, Ac-19}, pp. 716--722.
Aardom L., 1978. "Earth Rotation and Polar Motion from Laser Ranging to the Moon and Artificial Satellites", in: 9th Research Conference on Application of Geodesy to Geodynamics, Ohio State University, I.I.Mueller (ed.), Report of the Department of Geodetic Science, No. 280.
Barnes R.T.H., Hide R., White A.A., and Wilson C.A., (1983) Atmospheric Angular Momentum Fluctuations, length-of-day changes and polar motion, Proc. R. Soc. London, A387, 31--73.
Barrodale I., and Erickson R.E., 1980. "Algorithm for least-squares linear prediction and maximum entropy spectral analysis - Part I: Theory and Part II: Fortran program, Geophysics, 45, No 3, 420-446.
Bendat J.S., and A.G. Piersol, (1966) Random data: Analysis and Measurement Procedures, John Willey and Sons, Inc., New York.
Box G.E.P., Jenkins G. M. (1976). {\sl Time series analysis, forecasting and control.} Holden Day, San Francisco.
Brosche P. Seiler U., Sundermann J., and Wunsch J. 1989. "Periodic changes in Earth's rotation due to ocean tides", Astron. Astrophys. 220, 318-320.
Brosche P.,1990. "Tidal Deceleration of the Earth" in Variations in Earth Rotation, McCarthy and Carter, eds., Geophysical Monograph 59, AGU, 47-49.
Brzeziński A., (1987) Statistical Investigations on Atmospheric Angular Momentum Functions and their Effects on Polar Motion, Manuscripta Geodetica, 12, 268--281.
Brzeziński A., (1992) Polar Motion Excitation by Variations of the Effective Angular Momentum Functions: Considerations Concerning Deconvolution Problem, Manuscripta Geodetica, 17, 3--20.
Brzeziński A., (1994) "Polar motion excitation by variations of the effective angular momentum function, II: extended model", Manuscripta Geodetica, Vol. 19, pp. 157-171.
Burg J.P., 1967. "Maximum Entropy Spectral Analysis", 37th Ann. int. meet. Society of Exploration Geophysicists, Oklahoma City.
Chandler S.C., 1892. "On the Variation of Latitude", Astr. J., 267, pp 17-22.
Chao B.F., (1989) Length of Day Variations Caused by El-Nino Southern Oscillation and quasi biennial oscillation, Science, 243, 923--925.
Chao B.F., (1993) Excitation of Earth's Polar Motion by Atmospheric Angular Momentum Variations, 1980-1990, Geophys. Res. Lett., 20, No.2, 253--256.
Chao B.F., Ray R.D., Gipson J.M., Egbert G.D., Ma C. 1996. "Diurnal/Semidiurnal polar motion excited by oceanic tidal angular momentum", J. Geophys Res. Vol 101, No B9, 20151-20163.
Carter W.E., Robertson D.S., and Abell M.D. 1979. "An Improved Polar Motion and Earth Rotation Monitoring Service Using Radio Interferometry", in: Proc. of the IAU Symposium No. 82, on Time and the Earth's Rotation. D.D. McCarthy and J.D.H., Pilkington (eds.), Reidel. Publ. Company.
Chelliah M., 1990. "The Global Climate for June-August 1989: A Season of Near Normal Conditions in the Tropical Pacific", Journal of Climate, 3, 138-162.
Cooley J.W. and Tuckey J.W., 1965. "An algorithm for the machine calculation of comlex Fourier series", Math. Comp. Vol. 19, pp. 297-301.
Chui C.K., 1992. "An Introduction to Wavelets, Wavelet Analysis and its Applications", Vol. 1. Academic Press, Boston, San Diego.
Dickey J.O., Newhall X.X., Williams J.G., 1985. "Earth Orientation from Lunar Laser Ranging and an Error Analysis of Polar Motion Services", J. Geophys. Res. 90, 9353-9362.
Dickey J.O., 1989. "Axial Rate of Spin of the Earth", Mueller and Zerbini (eds), Proc. of the Int. Conf. Interdisciplinary Role of Space Geodesy, Erice, Italy.
Dickey J.O., 1993. "Atmospheric Excitation of the Earth's Rotation: Progress and Prospects via Space Geodesy", AGU Monograph, Contributions of Space Geodesy to Geodynamics: Earth Dynamics", Smith and Turcotte eds., Geodymanics Series, Vol. 24, 55-70.
Dickey J.O., Ghil M., Marcus S.L., 1991. "Extratropical Aspects of the 40-50 day Oscillation in Length-of-Day and Atmospheric Angular Momentum", J. Geophys. Res., 96, 22643-22658.
Dickman S.R., (1981) Investigation of Controversial Polar Motion Features Using Homogeneous International Latitude Service Data, J. Geophys. Res., 86, No B6, 4904--4932.
Dickman S.R., 1989. "Polar Motion", Proc. of Int. Conf. On the Interdisciplinary Role of Space Geodesy, Erice, Italy, Mueller I.I. and Zerbini S. (eds.) Springer Verlag, New York.
Eubanks T.M., Steppe J.A., Dickey J.O., Rosen R.D., Salstein D.A., (1988) Causes of Rapid Motions of the Earth's Pole, Nature, 334, No 6178.
Eubanks T.M., 1993. "Variations in the orientation of the Earth", AGU Monograph, Contributions of Space Geodesy to Geodynamics: Earth Dynamics, Smith and Turcotte eds., Geodynamics Series, Vol. 24, 1-54.
Farge M., Rabreau G., 1988. "Transformee en ondelettes pour detecter et analyser les structures coherentes dans les ecoulements turbulents bidimensionnels", C.R. Acad. Sci. Paris, t.307, Serie II, 1479-1486.
Feissel M. and Gambis D., (1980) La mise en evidence de variations rapides de la duree du jour, Compt. Rendus Hebdomadaires Sceances Acad. Sci., Ser. B, 291, 271--273.
Feissel M. and Lewandowski W., (1984) A comparative analysis of Vondrak and Gaussian smoothing techniques, Bull. Geod., 58, 464--474.
Freedman A.P., Steppe J.A., Dickey J.O., Eubanks T.M., Sung L.-Y., 1994. "The short term prediction of universal time and length of day using atmospheric angular momentum", J. Of Geophys. Res. Vol 99. No, B4, pp 6981-6996.
Gabor, 1946. "Theory of Communication", J. IEE (London) 93, 429-457.
Gambis D. (1991). M\'ethode d'analyse de s\'eries par transformation en ondolettes. Application \`a la rotation de la Terre. {\sl Journ\'ees Syst\`eme de r\'ef\'erence, Capitaine (Ed.)}, Paris Observatory, Paris, 147--152.
Gambis D. (1992). Wavelet transform analysis of the length of the day and the El-Ni\~no/Southern Oscillation variations at interseasonal and interannual time scales. Ann Geophysicae 10, 429--437.
Gasquet C., Witomski P., 1990. "Analyze de Fourier et Applications, Filtrage, Calcul Numerique, Ondelettes", Masson, Paris.
Gelb A. (ed.). 1988. "Applied Optimal Estimation", The MIT Press, Cabridge, Massachusetts, and London, England.
Goupillaud P., Grossmann A., Morlet J., 1984. "Cycle-octave and Related Transforms in Seismic Signal Analysis", Geoexploration 23, Elsevier, Amsterdam, 85-102.
Gubanov V.S., Petrov, S.D., Rusinov Yu. L., and Titov O.A., 1995. "Estimating, Combining, Filtering, and Forecasing of EOP by Least- Squares Collocation Technique", Paper presented at the IERS Workshop 1995.
Guillemain P., R. Kronland-Martinet, B. Martens, 1990. "Estimation of Spectral Lines with the Help of the Wavelet Transform", In Wavelets and Applications. Meyer Y. (ed.), Masson Paris.
Gross R.S. and U. Lindqwister, (1992) Atmospheric Excitation of Polar Motion During the GIG 91 Measurement Campaign, Geophys. Res. Lett., 19, 849--852.
Gross R.S., 1990. "The Secular Drift of the Rotation Pole" in the Earth Rotation and Coordinate Reference Frames, IAG Symp. 105, Boucher and Wilkins eds. 146-153.
Gross R.S. and U. Lindqwister, 1992. "Atmospheric Excitation of Polar Motion During the GIG 91 Measurement Campaign," Geophys. Res. Lett., 19, 849-852.
Gross R.S., 1995. "Combinations of Earth Orientation Measurements: SPACE94, COMB94, and POLE94", (submitted to Journal of Geophys. Res.)
Grossmann A. and J. Morlet, 1984. "Decomposition of Hardy Functions Into Square Integrable Wavelets of Constant Shape", SIAM J. Math. Anal.
Grossmann A. and J. Morlet, 1988. "Decomposition of Functions Into Wavelets of Constant Shape, and Related Transforms", Mathematics and Physics, Lectures on recent results, World Scientific.
Hamdan K. and Sung Li-Yu., 1996. "Stochastic modelling of length of day and universal time" J. of Geodesy, 70, 307-320.
Herring T.A., Dong D., King R.W., 1991. "Sub-milliarcsecond Determination of Pole Position using Global Positioning System Data", Geophys Res. Lett., 18, 1893-1897.
Hide R., Birch N.T., Morrison L. V., Shea D.J., and White A., (1980) Atmospheric Angular Momentum Fluctuations and Changes in the Length of Day, Nature, 286. 114--117.
Hide R., White A.A., and Wilson C.A., 1983. "Atmospheric Angular Momentum Fluctuations, length-of-day changes and polar motion", Proc. R. Soc. London, A387, 31-73.
Hozakowski W., (1989) Polar Motion Prediction by the Least-Squares Collocation Method, Proc. of the IAG Symp. No 105: Earth Rotation and Coordinate Reference Frames, Edinburgh, C. Boucher, G.A. Wilkins (eds), Springer Verlag, 50--57.
Hamdan K. and Sung Li-Yu., (1996) Stochastic modelling of length of day and universal time, J. of Geodesy, 70, 307--320.
Hvpfner J. 1995. "Periodische Anteile in der Erdrotation und dem atmosph rischen Drehimpuls und ihre Genauigkeiten", Zeitschrift fur Vermessungswesen, Heft 1, Januar 1995, 120. Jahrgang. 8-16.
IERS (1992). International Earth Rotation Service, NEOS Earth Orientation Bulletin A.
Iijima S. 1971. "On the Chandlerian and Annual Ellipses in the Polar Motionas obtained from every 12 Year Period", extra collection of papers contributed to the IAU Symp. No. 48, Rotation of the Earth, Mizusava, Japan.
Kołaczek B., Kosek w., 1985. "Analysis od Short Periodical Variations of Pole Coordinates Determined by Different Techniques in the MERIT Campaign", Proc. of the Int. Conf. On Earth Rotation and Terrestrial Reference Frame, Columbus, Ohio, USA, 505-524.
Kołaczek B., Kosek W., Nastula J., Hozakowski W. (1989). Variations of Polar Motion in 1986-1988, {\sl Proc. of the INTERCOSMOS Symposium "Use of Artificial Satellite Observations for Geodesy and geophysics"}, Cracow, 12-17 June, 1989, pp. 139--143.
Kołaczek B., Kosek W. (1991). Systematic Differences between radiointerferometric, and laser series of pole coordinates and short periodical variations of these coordinates. {\sl Polska Akademia Nauk, Centrum Badań Kosmicznych, Warszawa, Raport No 4}, Dec. 1991.
Kołaczek, B., Nastula, J., Gambis, D., Kosek, W., and W. Hozakowski, (1991) The semiannual oscillation of polar motion and its atmospheric excitation in the period of 1979-1989, Astron. Astrophys., 243, 276--182.
Kołaczek B., Kosek W. (1992). Variations of 120 days Oscillations of Earth Rotation and Atmospheric Angular Momentum. {\sl 7th International Symposium "Geodesy and Physics of the Earth" IAG Symposium No. 112}, Potsdam, Germany, Oct. 5-10, 1992.
Kołaczek B. and Kosek W., (1993) Variations of 80-120 days Oscillations of polar motion and Atmospheric Angular Momentum, Proc. 7th International Symposium "Geodesy and Physics of the Earth" IAG Symposium No. 112, Potsdam, Germany, 1992. H. Montag and Ch. Reigber (eds.), Springer Verlag, 439--442.
Kołaczek B. and Kosek W., (1994) Systematic Differences Between Radiointerferometric and Laser Series Of Pole Coordinates and Short Periodic Variations of Pole Coordinates, in Geophysical Monograph 82, IUGG Volume 17, on - Gravimetry and Space Techniques Applied to Geodynamics and Ocean Dynamics (B.E. Schutz ed.) International Union of Geodesy and Geophysics and the American Geophysical Union, 95--101.
Kołaczek B. and Kosek W., (1996) Outline of Characteristics of Polar Motion Determined by Observations, Publications of the Institute of Geophysics, Polish Academy of Sciences, M-18(273), Teisseyre R. (ed.), Warsaw, Poland, 163--170.
Kosek W., (1987) Computations of Short Periodical Variations of Pole Coordinates of the MERIT Campaign by the Use of Ormsby Filter and the Maximum Entropy Spectral Analysis with the Optimum Filter Length, Bull. Geod., 61, No 2., 109-124.
Kosek W. (1987b). Computation of Short Periodical Oscillations of Pole Coordinates Determined by Laser Technique in the MERIT Campaign Using Maximum Entropy Spectral Analysis and an Ormsby Band Pass Filter. {\sl Technical Report of Deutsches Geodatisches Forschungsinstitut in Munich. Interner Bericht, Ref.: WK/33/87/DGFI/Abt.I}, Feb. 1987.
Kosek W., (1989) The Computation of Short Periodic Variations in the Stochastic Geophysical Processes with Special Consideration of the Earth Rotation, Ph.D. Thesis, Institute of Geophysics of the Polish Academy of Sciences, Warsaw, Poland.
Kosek W. (1990). The computation of short periodic variations in the stochastic geophysical processes with special attention to the Earth Rotation. {\sl Ph.D. Thesis., Institute of Geophysics of the Polish Academy of Sciences}, Warsaw, Poland.
Kosek W. (1991). On high Frequency Spectral Analysis of very noisy stochastic Processes and its application to the ERP series determined by VLBI and SLR. {\sl U.S. Naval Observatory Circular No
177, Washington D.C., 20392-5100}, Feb 15.
Kosek W. (1992a). The Autocovariance Prediction of the Earth Rotation Parameters. {\sl 7th International Symposium "Geodesy and Physics of the Earth" IAG Symposium No. 112}, Potsdam, Germany, Oct
5-10, 1992.
Kosek W. (1992b). The Common Short Periodic Oscillations in the Earth's Rate of Rotation, Atmospheric Angular Momentum and Solar Activity, paper submitted to Bulletin Geodesique.
Kosek W. (1992) Short periodic autoregressive prediction of the Earth Rotation Parameters, Artificial Satellites, Planetary Geodesy, 27 No 2, 9--17.
Kosek W., Kaczkowski J. (1992). Short Periodic Oscillations in x and y pole coordinates of the SLR and VLBI Techniques Detected After Filtering with the Kalman Filter. {\sl 3-rd Orlov Conference "The Study of the Earth as Planet by Methods of Astronomy, Geophysics and Geodesy"} 7-12 Sep. 1992, Odessa, Ukraine.
Kosek W., (1993a) The Autocovariance Prediction of the Earth Rotation Parameters, 7th International Symposium Geodesy and Physics of the Earth, IAG Symposium No. 112, Potsdam, Germany, 1992. H. Montag and Ch. Reigber (eds.), Springer Verlag, 44
Kosek W., (1993b) The Common Short Periodic Oscillations in the Earth's Rate of Rotation, Atmospheric Angular Momentum and Solar Activity, Bull. Geod., 67, 1--9.
Kosek W., (1994) Short Period Oscillations in GPS Pole Coordinates Data, IERS Technical Note No 16 - Results from the SEARCH'92 Campaign., J.O. Dickey and M. Feissel (eds.) September 1994. Observatoire de Paris, pp. IV-61 -- IV-65.
Kosek W., (1995a) Time Variable Band Pass Filter Spectra of Real and Complex-Valued Polar Motion Series, Artificial Satellites, Planetary Geodesy, 30, No 1, 27--43.
Kosek W., (1995b) Analysis of Short Period GPS Pole Coordinates Data, Artificial Satellites, Planetary Geodesy, 30, No 2-3, 153--159.
Kosek W., Nastula J., Kołaczek B., 1995a. "Variability of Polar Motion Oscillations with Periods from 20 to 150 Days in 1979-1991", Bull. Geod., 69, 308-319.
Kosek W., (1997) Autocovariance Prediction of Short Period Earth Rotation Parameters (accepted to Artificial Satellites).
Kosek W., Nastula J., Kołaczek B., (1995a) Variability of Polar Motion Oscillations with Periods from 20 to 150 Days in 1979-1991, Bull. Geod., 69, 308--319.
Kosek W., McCarthy D.D., Luzum B., (1995b) Possible Improvement of Polar Motion Prediction Using Autocovariance Prediction Procedures, Proc. Journees 1995, Systemes de Reference Spatio-Temporels, Warsaw, Poland, 113--116.
Kosek W. and Kołaczek B., (1995) Irregular Short Period Variations of Polar Motion, Proc. Journees 1995, Systemes de Reference Spatio-Temporels, Warsaw, Poland, 1995. 117--120.
Kosek W., and Kołaczek B., (1996) Semi-Chandler and Semiannual Oscillations of Polar Motion, Geophys. Res. Lett.
Kosek W., McCarthy D.D., Luzum B., 1995b. "Possible Improvement of Polar Motion Prediction Using Autocovariance Prediction Procedures," Proc. Journees 1995 "Systemes de Reference Spatio-Temporels," Warsaw, Poland, 113-116.
Kosek W., McCarthy D.D., Luzum B., 1998. "Possible Improvement of Earth Rotation Forecas Using the Autocovariance Prediction Procedures", Journal of Geodesy.
Kosek W., Kołaczek B., 1995. "Irregular Short Period Variations of Polar Motion. Proc. Journees 1995 "Systemes de Reference Spatio-Temporels", Warsaw, Poland, Sep. 18-20, 1995, Published by the Space Research Centre, PAS, 117-120.
Kosek W., 1997. "Autocovariance Prediction of Short Period and Subseasonal Earth Rotation Parameters", Artificial Satellites.
Koopmans L.H., 1974. "Spectral Analysis of Time Series", Academic Press, New York.
Kuehne J.S., S.Johnson, and C.R.Wilson, 1993. "Atmospheric excitation of non-seasonal polar motion", J.Geophys. Res., 98, 191973-19978.
Lambeck K., (1980) The Earth's Variable Rotation: Geophysical Causes and Consequences., Cambridge University Press, Cambridge.
Langley R.B., King R.W., Shapiro I.I., Rosen R.D., Salstein D.A., 1981. "Atmospheric Angular Momentum and the Lenght of Day: A Common Fluctuation with a Period Near 50 Days", Nature, 294, 730-733.
Lichten S.M., Marcus S.M., Dickey J.O., 1992. "Sub-Daily Resolution of Earth Rotation Variations with Global Positioning System Measurements", Geoph. Res. Lett., 19, pp 537-540, 1992.
Luzum B.J., McCarthy D.D., Kosek W., (1997) New Prediction Algorithm for Polar Motion, (submitted to J. of Geodesy).
Markowitz W., 1970. "Sudden Changes in Rotational Acceleration of the Earth and Secular Motion of the Pole", (Mansinha L., Smylie D.E. and Beck A.E. eds.), Earthquake Displacement Fields and the Rotation of the Earth, Reidel, Dordrecht, Netherlands, pp 69-81.
Max J. (1981). Methodes et techniques de traitement du signal et applications aux measures physiques (Tome 1, -Principes Generaux et methodes classiques), Mason. S. A., Paris, France.
McCarthy D.D. and Luzum B.J., (1991) Prediction of Earth Orientation, Bull. Geod., 65, 18--21.
Morabito D.D., Eubanks T.M., Steppe J.A., 1988. "Kalman Filtering of Earth orientation changes", in Babcock and Wilkins (eds.), The Earth's Rotation and Reference Frames for Geodesy and Geodynamics, Kluwer, Dordrecht, pp. 257-267.
Moritz H., and Mueller I.I., 1987. "Earth Rotation, theory and observation", Ungar. Publ. Co., New York.
Nastula J., 1988. "Variations of Pole Motion Velocity and Their Correlation with Variations of E.A.M. Function in the Period 1985-1987. Report at the 6th International Symposium Geodesy and Physics of the Earth. Potsdam, Aug. 22-27 1988, GDR.
Nastula J., 1992. "Short Periodic Variations in the Earth's Rotation in the period 1984-1990", Ann. Geophys., 10, 441-448.
Nastula J., D. Gambis and M. Feissel, (1990) Correlated High Frequency Variations in Polar Motion and of the Length of the Day in Early 1988, Ann. Geophysicae, 8, 565--570.
Nastula J., A. Korsun, B. Kolaczek, W. Kosek, W. Hozakowski, (1993) Variations of the Chandler and Annual Wobbles of Polar Motion in 1846-1988 and their Prediction, Manuscripta Geodetica, 18, 131--135.
Nastula J. (1995) Short Periodic Variations of the Polar Motion and Atmospheric Angular Momentum Excitation Functions in the Period 1984-1992, Ann. Geophysicae, 13, 217--225.
Nastula J., Kosek W., Kołaczek B., 1996. "Analyses of Zonal Atmospheric Excitation Functions and Their Correlation with Polar Motion Excitation Functions" (submitted to Annales Geophysique)
Okubo S., (1982) Is the Chandler Period Variable?, Geophys. J.R. Astr. Soc., 71, 629--646.
Ormsby J.F.A., 1961. "Design of Numerical Filters with Application to Missile Data Processing", J. Assoc. Compt. Mach., 8, 440-466.
Otnes R.K. and Enochson L., (1972) Spectral Analysis and Time Series, Academic Press, London.
Petrov S., Brzezinski A., Gubanov V., 1996. "A Stochastic Model for Polar Motion with Application to Smoothing, Prediction and Combining", Artificial Satellites, Planetary Geodesy No 26, Vol. 31 No 1.
pp 51-70.
Popiński W., Kosek W., 1994. "The Wavelet Transform and Its Application for Short Period Earth Rotation Analysis", Artificial Satellites, Planetary Geodesy No 22 Vol. 29, No 2, SRC, Warsaw, Poland, pp 75-86.
Popiński W. and Kosek W., (1995a) The Fourier Transform Band Pass Filter and Its Application to Polar Motion Analysis, Artificial Satellites, Planetary Geodesy, 30 No 1, 9--25.
Popiński W. and Kosek W., (1995b) Discrete Fourier and Wavelet Transforms in Analysis of Earth Rotation Parameters, Proc. Journees 1995, Systemes de Reference Spatio-Temporels, Warsaw, Poland,
1995, 121--124.
Priestley M.B., (1981) Spectral Analysis and Time Series, Academic Press, London.
Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P., 1992. "Numerical Recipes in Fortran, The Art of Scientific Computing", Second Edition, Cambridge University Press, New York.
Price J., 1985. "Fourier Transform and Applications", Plenum Press. New York and London.
Rosen R.D. and Salstein D.A., (1983) Variations in Atmospheric Angular Momentum on Global and Regional Scales and the Length of Day, J. Geophys. Res., 88, No. C9, 5451--5470.
Rovelli A. and Vulpiani A., 1983. "Characteristic Correlation Time as Estimate of Optimum Filter Lenght in Maximum Entropy Spectral Analysis", Geophys J. R. astr.Soc. 72, 293-306.
Salstein D.A. and Rosen R.D., (1989) Regional Contributions to the Atmospheric Excitation of Rapid Polar Motions, J. Geophys. Res., 94, No. D7, 9971--9978.
Salstein, D.A. and R.D. Rosen, (1995) Angular momentum and energetics in reanalysis products, Proceedings of the Twentieth Annual Climate Diagnostics Workshop, US Dept. Commerce, NOAA NWS/CPC/NCEP, 337--340.
Singleton R.C., (1969) An Algorithm for Computing the Mixed Radix Fast Fourier Transform, IEEE Transactions on Audio and Electroacoustics, Vol. AU.-17, No.2 June 1969.
Speed T.P., 1985. "Some Practical and Statistical Aspects of Filtering and Spectrum Estimation", In Fourier Techniques and Applications, Price J.F. (ed.), Plenum Press, New York.
Stephenson F.R., and Morrison L.V., 1984. "Long term changes in the rotation of the Earth: 700 B.C. to A.D. 1980", Phil. Trans.,R. Soc. London., A 313, 47-70.
Stoyko A., Stoyko N., 1969. "Rotation de la Terre, phenomenes geophysiques et activitre du soleyil", Bull Acad. R.Belg. ser. 5, 55, pp 279-285.
Teisseyre R. (ed.) 1989. "Gravity and Low Frequency Dynamics", PWN-Polish Scientific Publishers, Warszawa, Elsevier, Amsterdam, Oxford, New York, Tokyo.
Tapley B.D., Schutz B.E., Eanes R.J., 1985. "Station Coordinates, Baselines and Earth Rotation from LAGEOS Laser Ranging: 1976-1984", J.Geophys Res., 90, 9235-9249.
Tretter S. A., (1976) Introduction to Discrete-time Signal Processing, John Wiley and Sons, New York.
Ulrych T. J., Bishop T. N. (1975). Maximum Entropy Spectral Analysis and Autoregressive Decomposition. Reviews of Geophysics and Space Physics.,Vol. 13., No 1. February (1975).
Ulrych T. J., Ooe M. (1979). Autoregressive and Mixed Autoregressive-Moving Average Models and Spectra, (ed.) Haykin S., Nonlinear Methods of Spectral Analysis., Topics in Applied Physics, Vol. 34, Springer-Verlag Berlin Heidelberg New York.
Wilson C.R. and Haubrich R.A., (1976) Meteorological Excitation of the Earth's Wobble, Geophys. J. R. Astr. Soc., 46, 707--743.
Vondrak J., (1985) Long period behavior of polar motion between 1900 and 1984, Ann. Geophys., 3, 351--356.
Vondrak J., (1989) Prediction of Polar Motion from Air and Water Excitation. Report No 402., Department of Geodetic Science and Surveying, The Ohio State University, Columbus, Dec, 1989.
Vondrak J., Ron C. 1995. "Long Period Oscillations in Earth Orientation Parameters: Comparative Analysis of Different Solutions", Proc. Journees 1995 "Systemes de Reference Spatio-Temporels", Warsaw, Poland, Sep. 18-20, 1995, Published by the Space Research Centre, PAS, 95-102.
Wilson C.R., Vicente R.O., 1980. "An Analysis of the Homogeneous Polar Motion Series", Geophys. J.R. astr. Soc., 62, 605-616.
Yoder C.F., Williams J. G., and Parke M.E., (1981) Tidal
Variations of Earth's Rotation, J. Geophys. Res., 86, No B2, 881--891.
Pytania egzaminacyjne:
1. Podaj przykłady zastosowania estymatora autokowariancji w analizach widmowych oraz prognozowaniu.
2. Jaka jest róznica pomiędzy analizą widmową metodą MESA oraz mtodą Blackman'a - Tuckey'a ?
3. Ogólna charakterystyka filtrów cyfrowych oraz ich podział.