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

Z Brain-wiki
 
(Nie pokazano 34 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 298: Linia 413:
 
i obliczamy widmo:
 
i obliczamy widmo:
 
<source lang = python>
 
<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===
 
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 &mdash; wartości wszystkich opisywanych modelem sygnałów (''x''<sub>1</sub>, ''x''<sub>2</sub>, ..., ''x''<sub>k</sub>) w chwili czasu ''t'' i wartości ''k'' niezależnych procesów czysto losowych (ϵ<sub>1</sub>, ϵ<sub>2</sub>, ..., ϵ<sub>k</sub>) w chwili czasu ''t'':
 
::<math>X(t)=\left(\begin{array}{l}
 
x_1(t)\\x_2(t)\\\vdots\\x_k(t)
 
\end{array} \right)\ \ \ E(t)=\left(\begin{array}{l}
 
\epsilon_1(t)\\
 
\epsilon_2(t)\\\vdots\\
 
\epsilon_k(t)
 
\end{array} \right)</math>
 
Proces opisujący ''jednocześnie'' więcej niż jeden sygnał nazywamy wielokanałowym lub wielozmiennym (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''&times;''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''&times;''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>
 
(symbol <sup>+</sup> oznacza transpozycję połączoną ze sprzężeniem zespolonym elementów macierzy).
 
 
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''&times;''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 ([http://www.fuw.edu.pl/~moira/other/dane_dtf_1.bin dane_dtf_1.bin], [http://www.fuw.edu.pl/~moira/other/dane_dtf_1.xml 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 ([http://www.fuw.edu.pl/~moira/other/dane_dtf_2.bin dane_dtf_2.bin], [http://www.fuw.edu.pl/~moira/other/dane_dtf_2.xml 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&rarr;2 i 1&rarr;3. Pytanie jest czy powinniśmy oczekiwać transmisji 2&rarr;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?
 
 
[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