Ćwiczenia 5

Model AR

Proces AR

Proces autoregresyjny rzędu [math]p[/math] jest zdefiniowany jako:

[math]x_t = \sum _{i=1}^p a_i x_{t-i} +\epsilon _t[/math]

Poniższy kod demonstruje generację pięciu realizacji procesu AR drugiego rzędu, każda o długości 500 punktów:

import numpy as np
import pylab as py
#wspolczynniki modelu AR

a = np.array([0.9, -0.6])

# chcemy badac N punktow procesu

for j in range(1,6):
    N=500
    x=np.zeros(N);
    #generujemy realizacje procesu
    for i in range(2,N):
        x[i]=a[0]*x[i-1]+a[1]*x[i-2]+np.random.randn()
    py.subplot(5,1,j)
    py.plot(x)
    py.title('realizacja'+str(j))
py.show()

Estymacja parametrów

Algorytmów służących do estymacji parametrów modelu AR jest kilka. Tu przedstawimy algorytm Yule-Walker'a:

• mnożymy stronami równania opisujące proces dla póbki [math]t[/math] i [math]t-m[/math]

[math]x_t x_{t-m} = \sum _{i=1}^p a_i x_{t-i} x_{t-m} +\epsilon _t x_{t-m} [/math]

• bierzemy wartość oczekiwaną lewej i prawej strony. Wartości oczekiwane [math]E\lbrace x_t x_{t-m}\rbrace [/math] to funkcja autokorelacji [math]R(m)[/math] więc:

[math]R(m) = \sum _{i=1}^p a_i R(m-i)+ \sigma _\epsilon ^2 \delta (m)[/math]

gdzie [math]m=0,\dots ,p.[/math]

• Dla [math]m\gt 0[/math] możemy zapisać stąd układ równań:

[math]\left[\begin{array}{c} R(1)\\ R(2)\\ \vdots \\ R(p) \end{array}\right]= \left[\begin{array}{cccc} R(0)& R(-1) &\dots &\\ R(1)& R(0) &R(-1) \dots &\\ \vdots & & &\\ R(p-1) & &\dots &R(0) \end{array}\right] \left[\begin{array}{c} a_1\\ a_2\\ \vdots \\ a_p \end{array} \right] [/math]

• stąd wyliczamy współczynniki [math]a[/math],

dla [math]m=0[/math] mamy

[math]R(0) = \sum _{k=1}^p a_k R(-k) + \sigma _\epsilon ^2[/math]

• można stąd wyliczyć [math]\sigma _\epsilon ^2 [/math]

Ćwiczenie estymacja parametrów procesu AR

• Wygeneruj 2000 próbek sygnału z modelu AR o parametrach a = {0.9, -0.6}, epsilon=2
• Oblicz funkcję autokorelacji tego sygnału: ak = np.correlate(x,x,mode='full')
• Znormalizuj tą funkcję dzieląc ją przez długość sygnału.
• Oblicz parametry zgodnie ze wzorami z poprzedniego paragrafu dla modelu rzędu 2. (wypisz konkretną postać wzorów analitycznie a następnie zaimplementuje je)

wskazówka: R[0]=ak[N-1]

• Wypisz parametry prawdziwe i estymowane.
• Sprawdź jak wpływa długość sygnału na dokładność estymaty (uruchom program kilka razy dla każdej z badanych długości sygnału)

*

W przypadku modelu rzędu [math]p[/math] estymację parametrów metodą Y-W można zaimplementować np. tak:

def parametryAR(x,p):
    '''funkcja estymująca parametry modelu AR 
    argumenty:
    x- sygnał
    p - rząd modelu
    f. zwraca:
    a - wektor współczynników modelu
    epsilon - estymowana wariancja szumu

    funkcja wymaga zaimportowania modułu numpy as np
    '''
    N = len(x)
    ak = np.correlate(x,x,mode='full')
    ak=ak/N
    R=ak[N-1:]
    RL  = R[1:1+p]
    RP = np.zeros((p,p))
    for i in range(p):
        aa = ak[N-1-i:N-1-i+p]
        RP[i,:] = aa
    a=np.linalg.solve(RP,RL)
    epsilon = (ak[N-1] - np.sum(a*ak[N:N+p]))**0.5
    return a,epsilon

Jak znaleźć rząd modelu?

Kryterium Akaike (AIC):

[math]\mathrm{AIC}(p)= \frac{2p}{N} +\ln(V) [/math]

[math]p[/math] - ilość parametrów modelu,

[math]N[/math] - ilość próbek sygnału,

[math]V[/math] - wariancja szumu.

Kryterium to karze za zwiększanie ilości parametrów i nagradza za zmniejszanie niewytłumaczonej wariancji.

Poniższy kod jest przykładową implementacją kryterium AIC:

def kryterium_AIC(x,maxymalnyRzad):
    zakres_rzedow = range(1,maxymalnyRzad)
    N = len(x)
    AIC = np.zeros(len(zakres_rzedow))
    for p in zakres_rzedow:
        a,epsilon = parametryAR(x,p)
        AIC[p-1] = (2.0*p)/N + np.log(np.sqrt(epsilon))
        print 'p:', p, ' a:',a,' epsilon: ',epsilon
    return AIC

Zobaczmy jak działa to na przykładowym syganle AR:

import numpy as np
import pylab as py
from numpy.fft import fft, fftshift
def parametryAR(x,p):
    '''funkcja estymująca parametry modelu AR 
    argumenty:
    x- sygnał
    p - rząd modelu
    f. zwraca:
    a - wektor współczynników modelu
    epsilon - estymowana wariancja szumu

    funkcja wymaga zaimportowania modułu numpy as np
    '''
    N = len(x)
    ak = np.correlate(x,x,mode='full')
    ak=ak/N
    R=ak[N-1:]
    RL  = R[1:1+p]
    RP = np.zeros((p,p))
    for i in range(p):
        aa = ak[N-1-i:N-1-i+p]
        RP[i,:] = aa
    a=np.linalg.solve(RP,RL)
    epsilon = (ak[N-1] - np.sum(a*ak[N:N+p]))**0.5
    return a,epsilon


def kryterium_AIC(x,maxymalnyRzad):
    zakres_rzedow = range(1,maxymalnyRzad)
    AIC = np.zeros(len(zakres_rzedow))
    for p in zakres_rzedow:
        a,epsilon = parametryAR(x,p)
        AIC[p-1] = (2.0*p)/N + np.log(np.sqrt(epsilon))
        print 'p:', p, ' a:',a,' epsilon: ',epsilon
    return AIC

#wspolczynniki modelu AR 
a = np.array([0.9, -0.7])
epsilon=2
N=600
x=np.zeros(N);
#generujemy realizacje procesu
for i in range(2,N):
    x[i]=a[0]*x[i-1]+a[1]*x[i-2] +epsilon*np.random.randn()

py.subplot(2,1,1)
py.plot(x)
py.title('wygenerowany sygnal')
py.subplot(2,1,2)

AIC = kryterium_AIC(x,6)
py.plot(range(1,len(AIC)+1),AIC)
py.title('Kryterium AIC')
py.show()

Widmo modelu AR

Widmo modelu można wyliczyć analitycznie znając jego współczynniki: Przepisujemy równanie modelu:

[math]x_t = \sum _{i=1}^p a_i x_{t-i} +\epsilon _t[/math]

[math]\sum _{i=0}^p a_i x_{t-i} =\epsilon _t[/math]

biorąc transformaty [math]Z[/math] obu stron mamy równanie algebraiczne:

[math]A(f)X(f) =E(f)[/math]

[math]X(f)=A^{-1}(f) E(f)=H(f) E(f)[/math]

Stąd widmo:

[math]S(f) = X(f)X^*(f)=H(f)VH^*(f)[/math]

Jak znaleźć [math]A[/math] — transformata Z

Transformata Z jest dyskretnym odpowiednikiem transformaty Laplace'a:

[math]X(z) = Z\lbrace x[n]\rbrace = \sum _{n=0}^\infty {x[n]z^{-n}}[/math]

gdzie [math]z=Ae^{i \phi }[/math] jest liczbą zespoloną. Szczególnym przypadkiem tej transformaty jest dyskretna transformata Fouriera - wystarczy podstawić [math]A=1/N[/math] i [math]\phi = - 2 \pi k/ N [/math] porównaj.

Własności transformaty Z

Transformata ta jest liniowa tzn.

[math]\mathrm{Z}\lbrace a_1x_1[n] +a_2x_2[n]\rbrace =a_1X_1(z)+a_2X_2(z)[/math]

jak ją policzyć od sygnału przesuniętego w czasie to:

[math]\mathrm{Z}\lbrace x[n-k]\rbrace = z^{-k}X(z)[/math]

dla impulsu:

[math]\mathrm{Z}\lbrace \delta [n]\rbrace =1[/math]

więc

[math]\mathrm{Z}\lbrace \delta [n-n0]\rbrace = z^{-n0} [/math]

Stosując tą transfomatę do procesu AR dostajemy:

[math]\mathrm{Z}\lbrace x[n] + a_1 x[n-1] + \dots + a_p x[n-p]\rbrace = (1 + a_1 z^{-1} + \dots + a_p z^{-p})X(z)[/math]

[math]=A(z)X(z)[/math]

Widmo procesu

Najpierw rozważmy konkretny przykład. Niech model będzie rzędu p=2 i ma współczynniki [math]a_1 = 0.9, a_2 = -0.6[/math] i [math]\varepsilon = 2[/math]. Wyliczamy wartości funkcji [math]A(z) = a_1 z^{-1}+a_2 z^{-2}[/math] dla [math]z = e^{i \omega}[/math]:

a=[0.9, -0.6]
epsilon = 2
w=np.arange(-np.pi,np.pi,0.1)
z=np.exp(1j*w)
# dla zadanego modelu
A=-1+a[0]*z**(-1)+a[1]*z**(-2);

Następnie obliczamy odwrotność A :

H=1./A

i obliczamy widmo:

Sp=H*H.conj()*epsilon**2
Sp = Sp.real

Możemy je wykreślić w funkcji częstości [math]\omega[/math].

py.plot(w,Sp )

Operacje te możemy uogólnić i zaimplementować jako funkcję do obliczania widma modelu zadanego przez parametry:

def widmoAR(parametry_a, epsilon, N_punktow):
    w = np.linspace(-np.pi,np.pi,N_punktow)
    z = np.exp(1j*w)
    A = -1 * np.ones(N_punktow) + 1j*np.zeros(N_punktow)
    for i in range(len(parametry_a)):
        A += parametry_a[i]*z**(-(i+1))
    H = 1./A
    Sp = H*H.conj()*epsilon**2 # widmo
    Sp = Sp.real
    return Sp, w

Ćwiczenie

Proszę:

• Wygenerować realizację modelu AR [math]a = \{0.6, -0.7, 0.3, -0.25\}, \quad \varepsilon = 2[/math]

def generujAR(a, epsilon, N):
    x=np.zeros(N)
    rzad = len(a)
    for i in range(rzad,N):
        for p in range(len(a)):
            x[i] += a[p]*x[i-(p+1)]
        x[i] += epsilon*np.random.randn()
    return x

• Obliczyć widmo dla tego modelu
• Wyestymować parametry modelu na podstawie sygału, zakładając, że rząd jest p = 3,4,5,6
• Obliczyć widmo dla wyestymowanego modelu
• Wykreślić widma prawdziwego modelu i modeli estymowanych

 *

Ćwiczenie

Dla modelu z poprzedniego ćwiczenia proszę wygenerować realizację sygnału długości 1000 punktów. Proszę porównać widma:

• prawdziwe, obliczone z prawdziwych parametrów modelu
• obliczone z estymowanego modelu
• obliczone przez periodogram
• obliczone metodą Welcha
• obliczone metodą wielookienkową Thomsona

import numpy as np
import pylab as py
from numpy.fft import fft, fftshift,fftfreq
import gendpss as dpss

def generujAR(a, epsilon, N):
    x=np.zeros(N)
    r = len(a)
    for i in range(r,N):
        for p in range(len(a)):
            x[i] += a[p]*x[i-(p+1)]
        x[i] += epsilon*np.random.randn()
    return x
    
def parametryAR(x,p):
    N =  len(x)
    ak = np.correlate(x,x,mode='full')
    ak = ak/N
    R = ak[N-1:]
    RL = R[1:1+p]
    RP = np.zeros((p,p))
    for i in range(p):
        aa = ak[N-1-i:N-1-i+p]
        RP[i,:] = aa
    a = np.linalg.solve(RP,RL)
    epsilon = (ak[N-1] - np.sum(a*ak[N:N+p]))**0.5
    return a,epsilon

def widmoAR(parametry_a, epsilon, N_punktow):
    w = np.linspace(-np.pi,np.pi,N_punktow)
    z = np.exp(1j*w)
    A = -1 * np.ones(N_punktow) + 1j*np.zeros(N_punktow)
    for i in range(len(parametry_a)):
        A += parametry_a[i]*z**(-(i+1))
    H = 1./A
    Sp = H*H.conj()*epsilon**2 # widmo
    Sp = Sp.real
    return Sp, w

def periodogram(s, okno , Fs):
    s = s*okno
    N_fft = len(s)
    S = fft(s,N_fft)#/np.sqrt(N_fft)
    P = S*S.conj()/np.sum(okno**2)
    P = P.real # P i tak ma zerowe wartośći urojone, ale trzeba ykonać konwersję typów
    F = fftfreq(N_fft, 1/Fs)
    return (fftshift(P),fftshift(F))

def pwelch(s,okno, przesuniencie, Fs):
    N = len(s)
    N_s = len(okno)
    
    start_fragmentow = np.arange(0,N-N_s+1,przesuniencie)
    ile_fragmentow = len(start_fragmentow)
    ile_przekrycia = N_s*ile_fragmentow/float(N)
    print ile_przekrycia, ile_fragmentow
    P_sredni = np.zeros(N_s)
    for i in range(ile_fragmentow):
        s_fragment = s[start_fragmentow[i]:start_fragmentow[i]+N_s]
        (P, F) = periodogram(s_fragment,okno,Fs)
        P_sredni += P
    return (P_sredni/ile_przekrycia,F)#(P_sredni/ile_przekrycia,F)

def mtm(s, NW = 3, Fs = 128):
    K = int(2*NW-1)
    N = len(s)
    w = dpss.gendpss(N,NW,K)
    S=np.zeros(N)
    for i in range(K):
        Si = np.abs(fft(s*w.dpssarray[i]))**2
        S[:] += Si.real
    S = S/K
    F = fftfreq(N,1.0/Fs)
    return (fftshift(S),fftshift(F))


#wspolczynniki modelu AR 
a = np.array([0.3, 0.2, 0.5, -0.25 ,-0.3])
epsilon = 2
N=256
# obliczanie widma z modelu
Sp, w = widmoAR(a,epsilon,N)

#generujemy realizacje procesu

x = generujAR(a, epsilon, N)

# estymujemy wspolczynniki modelu metodą Yula-Walkera
# obliczamy widmo dla estymowanego modelu
a_est,epsilon_est = parametryAR(x,5)
Sp_est, w = widmoAR(a_est,epsilon_est,N)

okno = np.blackman(N)
Fs = 2*np.pi                    
P_periodogram,F_periodogram = periodogram(x, okno , Fs=Fs)
okno = np.blackman(N/4)
P_welch, F_welch = pwelch(x,okno, len(okno)/2, Fs=Fs)                
P_mtm, F_mtm = mtm(x, NW = 4.5, Fs =Fs)

py.plot(w,Sp)
py.plot(w,Sp_est)
py.plot(F_periodogram,P_periodogram)
py.plot(F_welch, P_welch/(len(P_periodogram)/len(P_welch)))#uwaga Welch ma inne df
py.plot(F_mtm, P_mtm)

#py.legend(('prawdziwy','estymowany z AR','periodogram','Welch','mtm'))

print 'enenrgia sygnału: ', np.sum(x**2)
print 'enenrgia spektrum AR',np.sum(Sp)
print 'enenrgia est',np.sum(Sp_est)
print 'enenrgia mtm',np.sum(P_mtm)
print 'enenrgia welch',np.sum(P_welch)
print 'enenrgia period',np.sum(P_periodogram)
py.show()

Ogólny schemat

W praktyce musimy najczęściej estymować widmo mając dany tylko sygnał. Wówczas powinniśmy pwstąpić według następującego algorytmu:

• oszacować rząd modelu np. przy pomocy kryterium Akaikego
• wyestymować parametry modelu
• obliczyć widmo dla estymowanego modelu

Wielokanałowe modele AR

Model autoregresyjny rozpatrywaliśmy do tej pory jako proces generacji pojedynczego sygnału x. Nic jednak nie stoi na przeszkodzie, aby w chwili czasu t opisywać wartość nie jednego sygnału, ale jakiejś większej ich liczby, np. k. Możemy wtedy zapisać wzór modelu k-kanałowego jako

[math]X(t) = \sum _{i=1}^p A(i) X(t-i) +E(t)[/math]

gdzie X(t) i E(t) są wektorami odpowiednio — wartości wszystkich opisywanych modelem sygnałów (x₁, x₂, ..., x_k) w chwili czasu t i wartości k niezależnych procesów czysto losowych (ϵ₁, ϵ₂, ..., ϵ_k) w chwili czasu t:

[math]X(t)=\left(\begin{array}{l} x_1(t)\\x_2(t)\\...\\x_k(t) \end{array} \right)\ \ \ E(t)=\left(\begin{array}{l} \epsilon_1(t)\\ \epsilon_2(t)\\...\\ \epsilon_k(t) \end{array} \right)[/math]

Proces opisujący jednocześnie więcej niż jeden sygnał nazywamy wielokanałowym (ang. multichannel, multivariate). Zauważmy, że (jeden) model wielokanałowy nie jest prostym powieleniem kilku modeli jednokanałowych. Tutaj współczynniki A są macierzami rozmiaru k×k, zawierającymi oprócz zależności osobno dla każdego z sygnałów (kanałów modelu) również wyrazy pozadiagonalne, opisujące zależności między sygnałami składowymi.

Wzory opisujące model w dziedzinie częstości, w szczególności wzór na transformację Fouriera sygnału X oraz na widmo S (czyli funkcję gęstości widmowej mocy) wyglądają identycznie jak w przypadku jednokanałowym z tym, że odpowiednie wielkości (A, H, S) są teraz macierzami rozmiaru k×k:

[math]A(f)X(f) =E(f)[/math]

[math]X(f)=A^{-1}(f) E(f)=H(f) E(f)[/math]

[math]S(f) = X(f)X^*(f)=H(f)VH^*(f)[/math]

Macierz H nazywamy macierzą przejścia modelu (ang. transfer matrix). Widzimy, że transformata Fouriera sygnału powstaje przez pomnożenie macierzy przejścia H przez transformatę Fouriera szumu E, która teoretycznie nie zależy od częstości. Oznacza to, że własności widmowe opisane współczynnikami modelu (z których wylicza się macierz przejścia) zawarte są w H.

Widmo S procesu wielokanałowego jest również macierzą rozmiaru k×k. Na przekątnej zawiera ona tzw. widma własne (ang. autospectra), a poza przekątną widma wzajemne (ang. cross-spectra) mówiące o wspólnych dla danej pary kanałów składowych w częstości.

Ponieważ macierz H zawiera w sobie związki między wszystkimi sygnałami wchodzącymi w skład opisywanego układu oraz jest ona niesymetryczna, może posłużyć do skonstruowania funkcji opisującej zależności pomiędzy sygnałami. Funkcją taką jest na przykład kierunkowa funkcja przejścia (ang. directed transfer function, DTF) opisana w wersji znormalizowanej wzorem.

[math]\gamma_{ij}(f)=\frac{\left|H_{ij}(f)\right|^2}{\sum_{m=1}^{k}\left|H_{im}(f)\right|^2}[/math]

Opisuje ona stosunek wpływu w częstości f sygnału transmitowanego z kanału j do kanału i w danym zestawie w stosunku do wszystkich wpływów transmitowanych w tej częstości do kanału i. Dzięki temu jest ona znormalizowana: jej wartości zawierają się w przedziale [0, 1]; wartość 0 oznacza brak transmisji z kanału j do kanału i, wartość bliska 1 oznacza silny przepływ sygnału w tym kierunku.

Ćwiczenia

Aby zapoznać się z działaniem funkcji DTF możemy użyć wtyczki DTF do programu Svarog. Oblicza ona znormalizowaną funkcję DTF dla wybranego zestawu kanałów wczytanych z pliku i wybranego rzędu modelu.

Sygnał testowy 1 (dane_dtf_1.bin, dane_dtf_1.xml) jest to sygnał 3-kanałowy, w którym kanały 1 i 3 pochodzą z zapisu EEG ludzkiego z czaszki, pierwszy z przodu głowy, drugi z tyłu głowy. Kanał 2 jest to szum o rozkładzie normalnym. Sygnał zbierany był w trakcie czuwania z zamkniętymi oczami, oczekujemy więc silnej składowej w paśmie alfa (8-12 Hz), która, jak jest to wiadomo, jest generowana z tyłu głowy. Jak zinterpretujesz otrzymany wynik? Powinniśmy więc zasadniczo zaobserwować wpływ (transmisję) sygnału 3 na sygnał 1, przy niskich wartościach transmisji z i do pozostałych kanałów (oraz pomiędzy pozostałymi kanałami – szum powinien być niezależny od kanału 1 i 3).

Sygnał testowy 2 (dane_dtf_2.bin, dane_dtf_2.xml) jest to również sygnał 3-kanałowy, demonstrujący tzw. problem wspólnego źródła. Sygnał 1 pochodzi z rejestracji ludzkiego EEG podczas czuwania z zamkniętymi oczami, sygnał 2 to częściowo ten sam sygnał opóźniony o jedną próbkę z dodanym szumem, a sygnał 3 to sygnał 1 opóźniony o dwie próbki i z dodanym (innym) szumem. Ponieważ sygnały 2 i 3 są opóźnionymi wersjami sygnału 1 oczekujemy wykrycia transmisji 1→2 i 1→3. Pytanie jest czy powinniśmy oczekiwać transmisji 2→3? W końcu w kanale 3 jest sygnał podobny do sygnału z kanału 2 i na dodatek opóźniony w stosunku do niego o 1 próbkę.

Zaobserwuj wielkość wykrytej transmisji między sygnałami 2 a 3 w dwóch sytuacjach: gdy do analizy bierzesz tylko te dwa sygnały (2 i 3) oraz gdy do analizy brane są wszystkie trzy sygnały na raz. Czy możemy zaobserwować różnicę między sytuacją użycia pary sygnałów i informacji z pełnego układu 3-kanałowego?

[1]

Ćwiczenie

Teraz, kiedy umiemu już estymować widma różnymi metodami proszę pobawić się prawdziwymi synałami. Przykładowe sygnały i sposoby ich wczytywania podane są tu. Proszę zapisać sygnał c4spin.txt w swoim bieżącym katalogu, wczytać go i oszacować widmo metodą AR i metodą Welcha.