Ćwiczenia 4

Analiza_sygnałów_-_ćwiczenia/AR_1

Funkcja autokorelacji

Funkcja korelacji wzajemnej

Dla sygnału [math]x(t)[/math] możemy badać jak jest on do siebie podobny dla przesunięciaczasowego [math]\tau [/math]:

[math]\sigma _{xx}(\tau ) = \int {x(t) x(t+\tau ) dt} [/math]

Zwróćmy uwagę na związek funkcji korelacji:

• z iloczynem skalarnym wektorów [math]x(t)[/math] i jego przesuniętej wersji [math]x(t +\tau)[/math]
• ze splotem.

Dla skończonych dyskretnych sygnałów mamy estymatory korelacji:

[math]R_{xx}(m) = E\lbrace x_{n+m}x^*_n \rbrace = E\lbrace x_{n}x^*_{n-m} \rbrace [/math]

oraz kowariancji:

[math]C_{xx}(m) = E\lbrace (x_{n+m}-\bar{x}) (x_n-\bar{y})^*\rbrace = R_{xx}(m) - \bar{x} \bar{x}^*[/math]

Funkcję [math]R[/math] można estymować z jednej realizacji procesu (zakładamy jego ergodyczność):

[math] \widehat{R}_{xx}(m) = \left\lbrace \begin{array}{ll} \sum _{n=0}^{N-m-1}{x_{n+m} x_n^*} & m \ge 0 \\ \widehat{R}_{xx}^*(-m) & m \lt 0 \end{array} \right.[/math]

Zadanie: Jak obliczyć funkcję autokorelacji?

W pythonie mamy wydajną implementację funkcji autokorelacji (tak naprawdę to funkcji korelacji wzajemnej, której szczególnym przypadkiem jest autokorelacja) w: numpy.correlate. Zanim jednak przejdziemy do rutynowego korzystania z tej funkcji warto sobie wyrobić pewną intuicję co do jej działania. Możemy sobie to działanie wyobrazić tak, że tworzymy dwie kopie naszego sygnału wydłużone zerami. jedna kopia ma oryginalny sygnał po środku, a w drugiej będzie on stopniowo przesuwany od pozycji takiej w której z nieruchomą kopią ma jeden punkt wspólny z lewej stron do takiej, kiedy ma jeden punkt wspólny z prawej strony. W każdym położeniu wzajemnym sygnałów musimy policzyć ich iloczyn skalarny. Poniższy kod pomoże nam to zilustrować:

import numpy as np
import pylab as py

# Niech nasz sygnał będzie:
x = np.array([1,2,3,2,1])

# Ma on długość
N = ...

# przygotujmy miejsce na kopie sygnału
# najpierw kopia nieruchoma na środku 
x_2 = np.zeros(3*N-2, dtype = 'int')
x_2[...] = x
print('liczymy f. autokorelacji sygnału x= ',x)
for k in range(-N,N-1):
    print(' ')
    x_1 = np.zeros(3*N-2, dtype = 'int')
    x_1[...]=x  # to jest kopia sygnału x cofnięta o k próbek w stosunku do nieruchomej kopii
    print('x1: \t'+str(x_1))
    print('x2: \t'+str(x_2))
    print(40*'-')
    print('x1*x2\t', x_1*x_2)
    f_corr[k+N] = ...
    print('f_corr(',k+1,') = \t',...)
# porównajmy wynik z wywołaniem funkcji bibliotecznej:
...

Zadanie: Funkcja auokorelacji funkcji okresowej

Funkcja okresowa jest do siebie podobna po przesunięciu o cały okres.

• Jak zatem wygląda funkcja autokorelacji fragmentu funkcji sinus?
• Jak zmienia się ta funkcja wraz z długością fragmentu sygnału?

import numpy as np
import pylab as py

Fs = 100
dt = 1/Fs
f = 10
T = 0.2
t = np.arange(0,T,dt)
s = np.sin(2*np.pi*f*t)

f_corr = ...# pełna wersja funkcji autokorelacji - za pomoca funkcji bibliotecznej
tau  =...# wektor przesunięć wzajemnych - powinien mieć 0 na środku

py.subplot(2,1,1)
py.stem(t,s)

py.subplot(2,1,2)
py.stem(tau,f_corr)

py.show()

Zadanie: Funkcja auokorelacji białego szumuj

Analogicznie do zadania powyżej proszę zapoznać się z funknją autokorelacji dla białego szumu. Dla ustalenia uwagi niech to będą niezależne próbki z rozkładu N(0,1).

Procesy AR

Dla przypomnienia: proces AR generowany jest tak, że koljna 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]. Dla współczynników [math]a[/math], dla których proces ten jest stacjonarny powinien się także charakteryzować pewną konkretną funkcją autokorelacji. Poiższe ćwiczenie powinno nam uświadomić:

• jak mogą wyglądać pojedyncze realizacje procesu AR
• jak może wyglądać jego estymowana funkcja autokorelacji
• jek estymata funkcji autokorelacji zależy od długości realizacji (czyli od ilości dostępnych informacji)

...# stosowne importy

def realizacjaAR(a,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])

N_realizacji = 5 # liczba realizacji
N=500 #liczba punktów w realizacji
#co dzieje się z f_autokorelacji dla poszczególnych realizacji gdy zwiększamy N od 50 do 5000?

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()

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()

Analiza_sygnałów_-_ćwiczenia/AR_1