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

Z Brain-wiki
Linia 560: Linia 560:
 
\epsilon_k(t)
 
\epsilon_k(t)
 
\end{array} \right)</math>
 
\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''&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.
+
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'':
 
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'':

Wersja z 08:49, 19 sty 2016

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 (x1, x2, ..., xk) w chwili czasu t i wartości k niezależnych procesów czysto losowych (ϵ1, ϵ2, ..., ϵ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 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×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.