Metody Analiz Widmowych, Filtracji i Prognozowania

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:

  1. Interpolacja
  2. Filtracja
  3. Analizy widmowe
  4. Prognozowanie lub predykcja
  5. Analizy czasowo częstotliwościowe
Z interpolacją mamy do czynienia wówczas gdy chcemy uzupełnić proces stochastyczny w przedziale czasowym, w którym jest on określony. Estymacja danych poza przedział czasowy, w którym te dane są określone nazywana jest prognozowaniem lub predykcją. Prognozowanie to również określanie wartości danych przed początkiem przedziału czasowego, w których są one określone. Interpolacja i prognozowanie powinny być tak przeprowadzone aby charakterystyka częstotliwościowa szeregu czasowego interpolowanego lub przedłużonego wstecz lub w przód była taka sama jak szeregu czasowego oryginalnego.

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:

  1. wartość oczekiwana


  2. Wartość oczekiwana lub średnia dla przebiegów ciągłych zdefiniowana jest następująco:

    W praktyce dysponujemy skończonym czasem trwania procesu ciągłego zatem estymator wartości średniej dla przebiegów ciągłych zdefiniowany jest następująco:

    Wartość oczekiwana nazywana jest czasem momentem pierwszego rzędu. Estymator wartości oczekiwanej szeregu czasowego jest średnią arytmetyczną wartości tego szeregu

     
  3. autokowariancja


  4. Funkcja autokowariancji sygnału jest to czasowa funkcja korelacji sygnału z samym sobą przedstawiona w funkcji różnych przesunięć czasowych . Jeżeli sygnał jest okresowy wówczas jego funkcja autokowariancji jest również okresowa. Badanie funkcji autokowariancji sygnału pozwala określić po jakim czasie proces stochastyczny traci całkowicie korelacje ze swoimi poprzednimi wartościami. Czas ten określony jest poprzez wartość argumentu tej funkcji, dla którego po raz pierwszy osiąga ona wartość zero

    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.
     

  5. widmo mocy

Widmo mocy jest wariancją składowej procesu stochastycznego o częstotliwości zawartej w pewnym wąskim paśmie częstotliwości . Definicja widma mocy jest następująca:

gdzie jest teoretyczną składową lub oscylacją procesu o centralnej częstotliwości .

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:

  1. kowariancja wzajemna (cross-kowariancja)


  2. Funkcja kowariancji wzajemnej jest podobna do funkcji autokowariancji z tą różnicą, że opisuje korelację pomiędzy dwoma sygnałami w czasie. Funkcja kowariancji wzajemnej dwóch sygnałów losowych charakteryzuje wzajemne zależności wartości jednego sygnału losowego od wartości drugiego sygnału losowego.
    Dla stacjonarnych przebiegów ciągłych o zerowych wartościach oczekiwanych funkcja kowariancji wzajemnej określona jest następującym wzorem:

    Funkcja kowariancji wzajemnej niekoniecznie musi być funkcją parzystą i charakteryzuje się cechą asymetrii to znaczy gdy zamienimy miejscami to . Jeżeli to procesy stochastyczne są nieskorelowane.

    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 w tych szeregach są ze sobą nieskorelowane.
     

  3. widmo mocy wzajemnej (cross-spectrum)


  4. widmo mocy wzajemnej jest transformatą Fouriera kowariancji wzajemnej i wyrażone jest przez następujące równanie:

    Z własności asymetrii funkcji kowariancji wzajemnej wynika, że . Ponieważ funkcja kowariancji wzajemnej jest nieparzysta to widmo mocy wzajemnej jest funkcją zespoloną

    gdzie jest składową synfazową natomiast jest składową kwadraturową.
     
  5. koherencja

Koherencja jest funkcją częstotliwości pokazującą jaki jest współczynnik korelacji pomiędzy dwoma procesami stochastycznymi zależny od częstotliwości. Wysoka koherencja dla danej częstotliwości wskazuje na to, że dwa sygnały mają wysoką koncentrację mocy dla tej częstotliwości.
Następną własnością koherencji jest możliwość pomiaru wzajemnego kąta fazowego pomiędzy dwoma sygnałami w funkcji częstotliwości. Sygnały mogą mieć różnice faz równą zero i wówczas są one synchroniczne lub różnice 180 stopni i wówczas są asynchroniczne.
Funkcja koherencji pomiędzy dwoma przebiegami ciągłymi wyrażona jest następującym równaniem:

gdzie jest widmem mocy szeregu czasowego .

Koherencja zmienia się w przedziale od 0 do 1. Jeżeli dla wszystkich częstotliwości , wówczas szeregi czasowe są niekoherentne lub niezależne, nieskorelowane dla wszystkich częstotliwości. Jeżeli dla wszystkich częstotliwości, wówczas szeregi czasowe są całkowicie koherentne.
 
 
Wartość oczekiwana

-
Autokowariancja

Kowariancja wzajemna

Wariancja

Kowariancja

Autokorelacja

Korelacja wzajemna

Odchylenie standardowe

?
1
Współczynnik korelacji

Widmo mocy

Widmo mocy wzajemnej

1
Koherencja

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ą:

gdziejest 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 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, natomiastnazywamy 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 ostatnichpunktó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 prognozyjego 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: są oknami wagowym o długościach oraz .

Powyższa funkcja osiąga minimum jeżeli:
,
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 otrzymujemyestymatoró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 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ł.