Ćwiczenia 5

Analiza_sygnałów_-_ćwiczenia/AR_2

Procesy AR

Dla przypomnienia: proces AR generowany jest tak, że kolejna próbka jest ważoną sumą [math]p[/math] poprzednich próbek i niezależnej zmiennej losowej o średniej 0 i wariancji [math]\sigma^2[/math]:

[math]x_n = \sum_{k=1}^p a_k * x_{n-k} +\varepsilon_n[/math]

i [math]\varepsilon_n \sim N(0,\sigma^2)[/math]

Proces AR można zatem scharakteryzować podając:

• współczynniki [math]a[/math] oraz
• [math]\sigma^2[/math].

Zadanie: ilustracja realizacji procesu

Poniższy kod po uzupełnieniu będzie ilustrował jak mogą wyglądać pojedyncze realizacje procesu opisywanego przez współczynniki a=[0.9, -0.7] i wariancję 1.

...# stosowne importy

def realizacjaAR(a,sigma, N):
    x = np.zeros(N)
    for i in range(2,N): #kolejno tworzymy próbki w każdej realizacji
        x[i]=a[0]*x[i-1]+a[1]*x[i-2]+ ...
    return x
  
a = np.array([0.9, -0.7])
sigma = 1
N_realizacji = 5 # liczba realizacji
N=500 #liczba punktów w realizacji

realizacja = np.zeros((N_realizacji, N)); # macierz na wszystkie realizacje
for r in range(0,N_realizacji):    #generujemy realizacje procesu
    realizacja[r,:] = ...
    
for r in range(0,N_realizacji):   #rysujemy realizacje procesu
    py.subplot(5,1,r+1)
    py.plot( ...)
    py.title('realizacja'+str(r))
py.show()

Zadanie: ilustracja funkcji autokorelacji procesu

Dla współczynników [math]a[/math], dla których proces ten jest stacjonarny, charakteryzuje się też pewną konkretną funkcją autokorelacji. Poniższe ćwiczenie powinno nam uświadomić:

• jak może wyglądać estymowana funkcja autokorelacji dla realizacji procesu AR
• jak estymata funkcji autokorelacji zależy od długości realizacji (czyli od ilości dostępnych informacji). Co dzieje się z funkcją autokorelacji dla poszczególnych realizacji, gdy zwiększamy liczbę punktów w realizacji N od 50 do 5000?

Poniższy kod powinien stanowić kontynuację poprzedniego.

for r in range(0,N_realizacji):   #rysujemy funkcję autokorelacji poszczególnych realizacji
    f_corr =...
    tau = np.arange(-N+1,N,1)
    ind = range(...-10,...+10) # tu szykujemy indeksy, dzięki którym będziemy mogli pobrać wycinek +/- 10 próbek wokół przesunięcia 0 
    py.plot(tau[ind],f_corr[ind])
    py.title(fragment funkcji autokorelacji)
py.show()

Model AR

Rozważania na temat procesów AR są o tyle interesujące, że wiele sygnałów, które chcielibyśmy badać całkiem nieźle daje się opisać jako procesy AR . Wyobrażamy sobie wówczas, że rejestrowane sygnały są generowane przez pewien model AR (trochę tak jak funkcja realizacjaAR wytwarzała pojedyncze realizacje procesu). Pojawia się w tym momencie pytanie: jak możemy poznać wartości parametrów [math]a[/math] i [math]\sigma^2[/math], które pasują do badanych sygnałów?

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źć A — 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]\sigma_{\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]
sigma_eps = 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()* sigma_eps**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()* sigma_eps**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 \sigma_{\varepsilon} = 2[/math]

def generujAR(a, sigma_eps, 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] += sigma_eps*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

Ć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.