Ćwiczenia 5: Różnice pomiędzy wersjami

Z Brain-wiki
(Utworzono nową stronę "==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 k...")
 
 
(Nie pokazano 50 wersji utworzonych przez 4 użytkowników)
Linia 1: Linia 1:
==Model AR==
+
[[Analiza_sygnałów_-_ćwiczenia]]/AR_2
===Proces AR===
+
 
 +
= Notebook=
 +
Proszę pobrać i zapisać lkoalnie notebook:
 +
[https://drive.google.com/file/d/0BzwQ_Lscn8yDOG1icFRTeGMyRTQ/view?usp=sharing Ipython notebook z ćwiczeniami z modeli AR ]
 +
 
 +
 
 +
następnie w terminalu, w katalogu zawierającym notebook piszemy:
 +
> jupyter notebook
 +
 
 +
gdzie notebook to nazwa pobranego właśnie pliku
 +
 
 +
==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.
 +
<source lang = python>
 +
...# 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()
 +
</source>
 +
 
 +
===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.
 +
<source lang = python>
 +
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()
 +
</source>
  
 +
<!--
 
Proces autoregresyjny rzędu <math>p</math> jest zdefiniowany jako:
 
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>
 
::<math>x_t = \sum _{i=1}^p a_i x_{t-i} +\epsilon _t</math>
Linia 25: Linia 91:
 
py.show()
 
py.show()
 
</source>
 
</source>
 +
-->
 +
 +
==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===
 
===Estymacja parametrów===
 
 
Algorytmów służących do estymacji parametrów modelu AR jest kilka. Tu przedstawimy algorytm Yule-Walker'a:  
 
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>
 
* mnożymy stronami równania opisujące proces dla póbki <math>t</math> i <math>t-m</math>
Linia 57: Linia 127:
 
* można stąd wyliczyć <math>\sigma _\epsilon ^2 </math>
 
* można stąd wyliczyć <math>\sigma _\epsilon ^2 </math>
  
===Ćwiczenie estymacja parametrów procesu AR===
+
===Ćwiczenie: estymacja parametrów procesu AR===
 
* Wygeneruj 2000 próbek sygnału z modelu AR o parametrach a = {0.9, -0.6}, epsilon=2
 
* Wygeneruj 2000 próbek sygnału z modelu AR o parametrach a = {0.9, -0.6}, epsilon=2
 
* Oblicz funkcję autokorelacji tego sygnału:  <tt>ak = np.correlate(x,x,mode='full')</tt>
 
* Oblicz funkcję autokorelacji tego sygnału:  <tt>ak = np.correlate(x,x,mode='full')</tt>
* Znormalizuj tą funkcję dzieląc ją przez długość sygnału.
+
* Znormalizuj tą funkcję dzieląc ją przez liczbę pokrywających się próbek dla każdego przesunięcia <math>\tau</math>
* Oblicz parametry zgodnie ze wzorami z poprzedniego paragrafu dla modelu rzędu 2. (wypisz konkretną postać wzorów analitycznie a następnie zaimplementuje je)
+
* Oblicz parametry zgodnie ze wzorami z poprzedniego paragrafu dla modelu rzędu 2. (wypisz konkretną postać wzorów analitycznych a następnie zaimplementuj je)
wskazówka: <tt>R[0]=ak[N-1]</tt>
 
 
* Wypisz parametry prawdziwe i estymowane.
 
* 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)
 
* 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)
*
+
wskazówka: <tt>R[0]=ak[N-1]</tt>
 
<!--
 
<!--
 
<source lang = python>
 
<source lang = python>
Linia 114: Linia 183:
 
-->
 
-->
  
 +
<source lang = python>
 +
#wspolczynniki modelu AR
 +
a = np.array([0.9, -0.6])
 +
sigma = 2
 +
N = 200
 +
x = np.zeros(N);
 +
 +
#generujemy realizacje procesu
 +
for i in range(2,N):
 +
    x[i] = a[0]*x[i-1] + a[1]*x[i-...] + ...
 +
 +
py.subplot(2,1,1)
 +
py.plot(x)
 +
py.xlabel('numer próbki')
 +
py.title('wygenerowany sygnal')
 +
 +
py.subplot(2,1,2)
 +
ak = np.correlate(x,x,mode='full')
 +
# ak nieobciążona:
 +
norm_ak = ...
 +
ak /= norm_ak
 +
m = ... # przesunięcia
 +
py.plot(m, ak)
 +
py.xlabel('przesunięcie m')
 +
py.title('funkcja autokorelacij sygnalu x')
 +
 +
R=ak[N-1:]
 +
r0=R[0]
 +
r1=R[1]
 +
r2=R[2]
 +
 +
# estymujemy wspolczynniki modelu na podstawie funkncji autokorelacji
 +
 +
a2 = ...
 +
a1 = ...
 +
s_2 = ...
 +
 +
print('prawdziwe wspolczynniki')
 +
print(  a[0], a[1], sigma)
 +
print('estymowane wspolczynniki')
 +
print( '%.3f, %.3f, %.3f'%(a1,  a2, s_2**0.5))
 +
 +
py.show()
 +
</source>
 +
===Estymacja parametrów dla modelu rzędu p ===
 
W przypadku modelu rzędu <math>p</math> estymację parametrów metodą Y-W można zaimplementować np. tak:
 
W przypadku modelu rzędu <math>p</math> estymację parametrów metodą Y-W można zaimplementować np. tak:
 
<source lang = python>
 
<source lang = python>
Linia 129: Linia 243:
 
     N = len(x)
 
     N = len(x)
 
     ak = np.correlate(x,x,mode='full')
 
     ak = np.correlate(x,x,mode='full')
     ak=ak/N
+
    norm_ak = np.hstack((np.arange(1,N+1,1),np.arange(N-1,0,-1)))
 +
     ak=ak/norm_ak
 
     R=ak[N-1:]
 
     R=ak[N-1:]
 
     RL  = R[1:1+p]
 
     RL  = R[1:1+p]
Linia 137: Linia 252:
 
         RP[i,:] = aa
 
         RP[i,:] = aa
 
     a=np.linalg.solve(RP,RL)
 
     a=np.linalg.solve(RP,RL)
     epsilon = (ak[N-1] - np.sum(a*ak[N:N+p]))**0.5
+
     sigma = (ak[N-1] - np.sum(a*ak[N:N+p]))**0.5
     return a,epsilon
+
     return a, sigma
 
</source>
 
</source>
  
Linia 147: Linia 262:
 
::<math>\mathrm{AIC}(p)= \frac{2p}{N} +\ln(V) </math>
 
::<math>\mathrm{AIC}(p)= \frac{2p}{N} +\ln(V) </math>
  
<math>p</math> - ilość parametrów modelu,
+
<math>p</math> - liczba parametrów modelu,
  
<math>N</math> - ilość próbek sygnału,
+
<math>N</math> - liczba próbek sygnału,
  
 
<math>V</math> - wariancja szumu.
 
<math>V</math> - wariancja szumu.
Linia 163: Linia 278:
 
     AIC = np.zeros(len(zakres_rzedow))
 
     AIC = np.zeros(len(zakres_rzedow))
 
     for p in zakres_rzedow:
 
     for p in zakres_rzedow:
         a,epsilon = parametryAR(x,p)
+
         a,sigma = parametryAR(x,p)
         AIC[p-1] = (2.0*p)/N + np.log(np.sqrt(epsilon))
+
         AIC[p-1] = (2.0*p)/N + np.log(np.sqrt(sigma))
         print 'p:', p, ' a:',a,' epsilon: ',epsilon
+
         print 'p:', p, ' a:',a,' sigma: ',sigma
 
     return AIC
 
     return AIC
  
 
</source>
 
</source>
  
Zobaczmy jak działa to na przykładowym syganle AR:
+
Zobaczmy jak działa to na przykładowym sygnale AR:
 
<source lang = python>
 
<source lang = python>
 
import numpy as np
 
import numpy as np
Linia 249: Linia 364:
 
::<math>S(f) = X(f)X^*(f)=H(f)VH^*(f)</math>
 
::<math>S(f) = X(f)X^*(f)=H(f)VH^*(f)</math>
  
=== Jak znaleźć <math>A</math> &mdash; transformata Z===
+
=== Jak znaleźć ''A'' &mdash; transformata Z===
  
 
Transformata Z jest dyskretnym odpowiednikiem transformaty Laplace'a:
 
Transformata Z jest dyskretnym odpowiednikiem transformaty Laplace'a:
Linia 282: Linia 397:
  
 
===Widmo procesu===
 
===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>:
+
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>:
 
<source lang = python>
 
<source lang = python>
 
a=[0.9, -0.6]
 
a=[0.9, -0.6]
epsilon = 2
+
sigma_eps = 2
 
w=np.arange(-np.pi,np.pi,0.1)
 
w=np.arange(-np.pi,np.pi,0.1)
 
z=np.exp(1j*w)
 
z=np.exp(1j*w)
Linia 297: Linia 412:
 
</source>
 
</source>
 
i obliczamy widmo:
 
i obliczamy widmo:
<source lang = py >
+
<source lang = python>
Sp=H*H.conj()*epsilon**2
+
Sp=H*H.conj()* sigma_eps**2
 
Sp = Sp.real
 
Sp = Sp.real
 
</source>
 
</source>
Linia 315: Linia 430:
 
         A += parametry_a[i]*z**(-(i+1))
 
         A += parametry_a[i]*z**(-(i+1))
 
     H = 1./A
 
     H = 1./A
     Sp = H*H.conj()*epsilon**2 # widmo
+
     Sp = H*H.conj()* sigma_eps**2 # widmo
 
     Sp = Sp.real
 
     Sp = Sp.real
 
     return Sp, w
 
     return Sp, w
Linia 322: Linia 437:
 
==== Ćwiczenie ====
 
==== Ćwiczenie ====
 
Proszę:
 
Proszę:
* Wygenerować realizację modelu AR <math>a = \{0.6, -0.7, 0.3, -0.25\}, \quad \varepsilon = 2</math>
+
* Wygenerować realizację modelu AR <math>a = \{0.6, -0.7, 0.3, -0.25\}, \quad \sigma_{\varepsilon} = 2</math>
 
<source lang = python>
 
<source lang = python>
def generujAR(a, epsilon, N):
+
def generujAR(a, sigma_eps, N):
 
     x=np.zeros(N)
 
     x=np.zeros(N)
 
     rzad = len(a)
 
     rzad = len(a)
Linia 330: Linia 445:
 
         for p in range(len(a)):
 
         for p in range(len(a)):
 
             x[i] += a[p]*x[i-(p+1)]
 
             x[i] += a[p]*x[i-(p+1)]
         x[i] += epsilon*np.random.randn()
+
         x[i] += sigma_eps*np.random.randn()
 
     return x
 
     return x
 
</source>
 
</source>
Linia 345: Linia 460:
 
from numpy.fft import fft, fftshift
 
from numpy.fft import fft, fftshift
  
def generujAR(a, epsilon, N):
+
def generujAR(a, sigma_eps, N):
 
     x=np.zeros(N)
 
     x=np.zeros(N)
 
     rzad = len(a)
 
     rzad = len(a)
Linia 351: Linia 466:
 
         for p in range(len(a)):
 
         for p in range(len(a)):
 
             x[i] += a[p]*x[i-(p+1)]
 
             x[i] += a[p]*x[i-(p+1)]
         x[i] += epsilon*np.random.randn()
+
         x[i] += sigma_eps*np.random.randn()
 
     return x
 
     return x
 
      
 
      
Linia 379: Linia 494:
 
     return a,epsilon
 
     return a,epsilon
  
def widmoAR(parametry_a, epsilon, N_punktow):
+
def widmoAR(parametry_a, sigma_eps, N_punktow):
 
     w = np.linspace(-np.pi,np.pi,N_punktow)
 
     w = np.linspace(-np.pi,np.pi,N_punktow)
 
     z = np.exp(1j*w)
 
     z = np.exp(1j*w)
Linia 386: Linia 501:
 
         A += parametry_a[i]*z**(-(i+1))
 
         A += parametry_a[i]*z**(-(i+1))
 
     H = 1./A
 
     H = 1./A
     Sp = H*H.conj()*epsilon**2 # widmo
+
     Sp = H*H.conj()* sigma_eps**2 # widmo
 
     Sp = Sp.real
 
     Sp = Sp.real
 
     return Sp, w
 
     return Sp, w
Linia 394: Linia 509:
 
epsilon = 2
 
epsilon = 2
 
# obliczanie widma z modelu
 
# obliczanie widma z modelu
Sp, w = widmoAR(a,epsilon,200)
+
Sp, w = widmoAR(a, sigma_eps,200)
  
 
#generujemy realizacje procesu
 
#generujemy realizacje procesu
 
N=600
 
N=600
x = generujAR(a, epsilon, N)
+
x = generujAR(a, sigma_eps, N)
  
 
# estymujemy wspolczynniki modelu metodą Ypula-Walkera
 
# estymujemy wspolczynniki modelu metodą Ypula-Walkera
 
# obliczamy widmo dla estymowanego modelu
 
# obliczamy widmo dla estymowanego modelu
 
for p in range(3,7):
 
for p in range(3,7):
     a_est,epsilon_est = parametryAR(x,p)
+
     a_est, sigma_eps_est = parametryAR(x,p)
     Sp_est, w = widmoAR(a_est,epsilon_est,200)
+
     Sp_est, w = widmoAR(a_est, sigma_eps_est,200)
 
     py.plot(w,Sp_est )
 
     py.plot(w,Sp_est )
 
py.plot(w,Sp)
 
py.plot(w,Sp)
Linia 541: Linia 656:
  
 
</source>
 
</source>
 
+
== Co z tego musimy zapamiętać:==
====Ogólny schemat ====
+
* widmo sygnału stochastycznego estymujemy, a nie obliczamy
 +
* są dwie klasy technik:
 +
** nieparametryczne - '''widmo estymujemy''' bezpośrednio dla sygnału np. metodą periodogram, Welcha, wielookienkową
 +
**  parametryczne: najpierw '''estymujemy model''' opisujący dane, a nstępnie '''dla modelu obliczamy widmo'''
 +
===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:
 
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
 
* oszacować rząd modelu np. przy pomocy kryterium Akaikego
 
* wyestymować parametry modelu
 
* wyestymować parametry modelu
 
* obliczyć widmo dla estymowanego modelu
 
* obliczyć widmo dla estymowanego modelu
 
===Wielokanałowe modele AR===
 
[http://eeg.pl/DTF]
 
 
==Ć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ą
 
[http://brain.fuw.edu.pl/edu/TI:Programowanie_z_Pythonem/Wejście_i_wyjście#Przyk.C5.82ady tu].
 
Proszę zapisać sygnał c4spin.txt w swoim bieżącym katalogu, wczytać go i oszacować widmo metodą AR i metodą Welcha.
 

Aktualna wersja na dzień 08:51, 30 sty 2017

Analiza_sygnałów_-_ćwiczenia/AR_2

Notebook

Proszę pobrać i zapisać lkoalnie notebook: Ipython notebook z ćwiczeniami z modeli AR


następnie w terminalu, w katalogu zawierającym notebook piszemy: > jupyter notebook

gdzie notebook to nazwa pobranego właśnie pliku

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 liczbę pokrywających się próbek dla każdego przesunięcia [math]\tau[/math]
  • Oblicz parametry zgodnie ze wzorami z poprzedniego paragrafu dla modelu rzędu 2. (wypisz konkretną postać wzorów analitycznych a następnie zaimplementuj je)
  • 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)

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

#wspolczynniki modelu AR 
a = np.array([0.9, -0.6])
sigma = 2
N = 200
x = np.zeros(N);

#generujemy realizacje procesu
for i in range(2,N):
    x[i] = a[0]*x[i-1] + a[1]*x[i-...] + ...

py.subplot(2,1,1)
py.plot(x)
py.xlabel('numer próbki')
py.title('wygenerowany sygnal')

py.subplot(2,1,2)
ak = np.correlate(x,x,mode='full')
# ak nieobciążona:
norm_ak = ...
ak /= norm_ak
m = ... # przesunięcia
py.plot(m, ak)
py.xlabel('przesunięcie m')
py.title('funkcja autokorelacij sygnalu x')

R=ak[N-1:]
r0=R[0]
r1=R[1]
r2=R[2]

# estymujemy wspolczynniki modelu na podstawie funkncji autokorelacji

a2 = ...
a1 = ...
s_2 = ...

print('prawdziwe wspolczynniki')
print(  a[0], a[1], sigma)
print('estymowane wspolczynniki')
print( '%.3f, %.3f, %.3f'%(a1,  a2, s_2**0.5))

py.show()

Estymacja parametrów dla modelu rzędu p

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')
    norm_ak = np.hstack((np.arange(1,N+1,1),np.arange(N-1,0,-1)))
    ak=ak/norm_ak
    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)
    sigma = (ak[N-1] - np.sum(a*ak[N:N+p]))**0.5
    return a, sigma

Jak znaleźć rząd modelu?

Kryterium Akaike (AIC):

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

[math]p[/math] - liczba parametrów modelu,

[math]N[/math] - liczba 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,sigma = parametryAR(x,p)
        AIC[p-1] = (2.0*p)/N + np.log(np.sqrt(sigma))
        print 'p:', p, ' a:',a,' sigma: ',sigma
    return AIC

Zobaczmy jak działa to na przykładowym sygnale 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()

Co z tego musimy zapamiętać:

  • widmo sygnału stochastycznego estymujemy, a nie obliczamy
  • są dwie klasy technik:
    • nieparametryczne - widmo estymujemy bezpośrednio dla sygnału np. metodą periodogram, Welcha, wielookienkową
    • parametryczne: najpierw estymujemy model opisujący dane, a nstępnie dla modelu obliczamy widmo

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