http://brain.fuw.edu.pl/edu/api.php?action=feedcontributions&user=Prozanski&feedformat=atomBrain-wiki - Wkład użytkownika [pl]2024-03-29T14:50:16ZWkład użytkownikaMediaWiki 1.34.1http://brain.fuw.edu.pl/edu/index.php?title=%C4%86wiczenia_1&diff=5701Ćwiczenia 12016-10-11T08:58:44Z<p>Prozanski: zmiana linku</p>
<hr />
<div>[[Analiza_sygnałów_-_ćwiczenia]]/Sygnały<br />
<br />
==Narzędzia wykorzystywane na ćwiczeniach==<br />
===Python===<br />
<br />
Python jest językiem programowania wysokiego poziomu, który w połączeniu z bibliotekami NumPy i SciPy do obliczeń naukowych pozwala na szybkie i wygodne programowanie lub analizowanie danych w sposób interaktywny. Przykłady prezentowane w ramach zajęć powinny uruchamiać się zarówno w wersji 2 jak i 3 języka Python, jednak zachęcamy Państwa, aby od początku uczyć się i korzystać z wersji 3 języka.<br />
<br />
Szczególnie przydatne na analizie sygnałów będą moduły:<br />
*numpy<br />
*scipy<br />
*matplotlib<br />
<br />
Do reprezentowania sygnałów w programach będziemy stosować tablice numpy: <br />
* <tt>ndarray</tt> <br />
Jest to zarówno efektywne jeśli chodzi o pamięć jak i o szybkość wykonywania operacji matematycznych. <br />
<br />
==== Dokumentacja modułu scipy.signal ====<br />
Proszę zapoznać się z dokumentacją biblioteki scipy.signal:<br />
<br />
https://docs.scipy.org/doc/scipy/reference/<br />
<br />
===Svarog===<br />
Przydatnym narzędziem do analizy sygnałów, z którego będziemy korzystać na zajęciach, jest program SVAROG (pierwotnie skrót od Signal Viewer, Analyzer and Recorder On GPL). Program działa w środowisku Java, jest więc niezależny od systemu operacyjnego (Linux, Windows, OS X…). Svarog pozwala na wczytywanie i analizowanie sygnałów (nie tylko bioelektrycznych), zarówno przy użyciu prostych (FFT, spektrogram) jak i bardziej zaawansowanych (matching pursuit, ICA, DTF itd.) narzędzi. Dzięki współpracy z platformą OpenBCI, możliwa jest rejestracja sygnału (łącznie z metadanymi) bezpośrednio z poziomu graficznego interfejsu użytkownika.<br />
<br />
====Svarog: uruchamianie i konfiguracja====<br />
Aktualną wersję programu Svarog można pobrać [http://www.fuw.edu.pl/~durka/svarog/ stąd]. Program nie wymaga instalacji. Po rozpakowaniu paczki do dowolnego katalogu należy uruchomić skrypt „run-svarog.sh” lub uruchomić bezpośrednio plik *.jar.<br />
<br />
W przypadku pracy na własnych komputerach, do prawidłowego uruchomienia pluginu do analizy sygnałów, z którego będziemy korzystać w dalszej części ćwiczeń, konieczne jest zainstalowanie środowiska Oracle Java SE w wersji 8, które można pobrać [http://www.oracle.com/technetwork/java/javase/downloads/index.html ze strony wydawcy]. Alternatywnie, użytkownicy systemu Ubuntu lub pokrewnych dystrybucji mogą zainstalować środowisko Java według instrukcji [http://www.webupd8.org/2012/09/install-oracle-java-8-in-ubuntu-via-ppa.html dostępnych na tej stronie].<br />
<br />
==Sygnały ciągłe i dyskretne ==<br />
===Próbkowanie w czasie ===<br />
W tym ćwiczeniu zilustrujemy pojęcia:<br />
* częstość próbkowania<br />
* częstość Nyquista<br />
* aliasing<br />
<br />
====Próbkowanie====<br />
W poniższym ćwiczeniu chcemy zbadać efekt próbkowania sygnału w czasie. <br />
*W komputerach nie mamy dostępu do sygnału ciągłego. <br />
*Wartości sygnału znane są tylko w dyskretnych momentach czasu. <br />
*Najczęściej stosujemy równe odstępy czasu <math>dt</math><br />
*Odwrotność tych okresów to częstość próbkowania <math>Fs = \frac{1}{dt}</math><br />
*Bardzo wygodnie jest myśleć o sygnałach jako o wektorach lub macierzach (rys. na tablicy)<br />
*Do przechowywania sygnałów w pamięci używamy ndarray<br />
<br />
<br />
Czy próbkując sygnał z częstością <math>Fs = 100</math>[Hz] mogę odwzorować sygnał o dowolnej częstości?<br />
* wytwórz wektor t reprezentujący czas 1s próbkowany z częstością Fs<br />
<source lang =python><br />
# -*- coding: utf-8 -*-<br />
<br />
import pylab as py<br />
import numpy as np<br />
Fs =100<br />
dt = 1/Fs<br />
t = np.arange(0,1,dt) # czas 'prawie ciągły'<br />
</source><br />
<br />
* wytwórz sygnał s, sinus o częstości f =10Hz. Dla przypomnienia wyrażenie: <math> s(t) = \sin(2 \pi f t)</math> możemy w pythonie zapisać:<br />
<source lang =python><br />
s = np.sin(2*np.pi*f*t)<br />
</source><br />
<br />
* wykreśl ten sygnał za pomocą punktów i linii<br />
* wykreśl sygnały o częstościach 20, 40, 50, 90 Hz<br />
<br />
====Efekt aliasingu====<br />
Rzućmy okiem na przykład tego efektu: https://youtu.be/SFbINinFsxk<br />
* koła obracają się z pewną częstością zgodnie z kierunkiem ruchu wskazówek zegara;<br />
* na filmie widać, że w pewnych kierunek ten się zmienia, mimo, że samochód wciąż jedzie w tą samą stronę;<br />
* możemy przyjąć umownie, że obrotom zgodnie z kierunkiem wskazówek zegara przypiszemy wartości dodatnie, a obrotom w kierunku przeciwnym wartości ujemne.<br />
* Dlaczego tak się dzieje?<br />
Analogiczne zjawisko przeanalizujemy dla symulowanych sygnałów okresowych. <br />
Na nasze potrzeby wygenerujemy sygnały próbkowane z bardzo dużą częstością, które będą dla nas aproksymacją sygnałów ciągłych. Przy ich pomocy zaprezentujemy efekt utożsamiania (aliasingu). <br />
*Proszę wytworzyć wektor reprezentujący czas &bdquo;prawie&rdquo; ciągły. Będzie to u nas 1000 wartości z przedziału [0,1) wziętych z odstępem 0,001.<br />
*Teraz proszę wygenerować dwie sinusoidy: jedną o częstości <tt>-1</tt> a drugą o częstości <tt>9</tt>. <br />
*Proszę [[TI/Matplotlib#Kilka_wykres.C3.B3w_we_wsp.C3.B3lnych_osiach|wykreślić]] obie sinusoidy.<br />
*Teraz proszę spróbkować czas i nasze &bdquo;prawie&rdquo; ciągłe sinusoidy z okresem próbkowania 0,1. (Trzeba pobrać co 100 element, proszę posłużyć się [[Programowanie_z_Pythonem/Sekwencje#Wycinki|wycinkami]]) <!-- [[TI/Sekwencje|Struktury danych — sekwencje|wycinkami]]. --><br />
*Na tle &bdquo;prawie&rdquo; ciągłych sinusoid proszę dorysować punkty ze spróbkowanych sygnałów. Aby punkty były dobrze widoczne proponuję użyć markerów <tt>x</tt> oraz <tt>+</tt>.<br />
*Proszę zaobserwować wzajemne położenie punktów. <br />
*Czy można odróżnić sinusoidę o częstości &minus;1 od sinusoidy o częstości 9, jeśli obie są próbkowane z częstością 10? <br />
*Jak można uogólnić tą obserwację?<br />
<br />
<tt>*</tt><br />
<!--<br />
# -*- coding: utf-8 -*-<br />
<br />
import pylab as py<br />
import numpy as np<br />
<br />
t = np.arange(0,1,0.001) # czas 'prawie ciągły'<br />
f1 = -1 # częstość sygnału 1<br />
f2 = 9 # częstość sygnału 2<br />
s1 = np.sin(2*np.pi*f1*t) # prawie ciągły sygnał o częstości f1<br />
s2 = np.sin(2*np.pi*f2*t) # prawie ciągły sygnał o częstości f2<br />
<br />
T = 0.1 # okres próbkowania<br />
Fs = 1/T # częstość próbkowania<br />
FN = Fs/2 # częstość Nyquista<br />
T_samp = t[0::100] # czas pobierania próbek<br />
s1_samp = s1[0::100] # próbkowany sygnał o częstości f1<br />
s2_samp = s2[0::100] # próbkowany sygnał o częstości f2<br />
<br />
py.plot(t, s1, 'g')<br />
py.plot(t, s2, 'b')<br />
py.plot(T_samp, s1_samp, 'gx', markersize=10)<br />
py.plot(T_samp, s2_samp, 'r+', markersize=10)<br />
--><br />
<br />
<!--<br />
===Błąd kwantyzacji ===<br />
Kiedy mierzymy fizyczne wielkości w celu dalszej analizy najczęściej chcemy przypisać im pewne liczby. Liczby w systemach cyfrowych reprezentowane są ze skończoną dokładnością. Urządzenia dokonujące przypisania liczby do mierzonej wartości to przetworniki analogowo-cyfrowe (ang. ADC, Analog to Digital Converter).<br />
Charakteryzują się one określoną ilością bitów ''N'', za pomocą których reprezentują liczby. Pełen zakres wartości pomiarowych ''R'' jest dzielony na <math>2^N</math> poziomów. Błąd kwantyzacji szacujemy jako nie większy niż <math> \frac{R}{2^{N+1}} </math>. <br />
<br />
====Zadanie:====<br />
Proszę oszacować błąd kwantyzacji sygnału EEG mierzonego za pomocą 12-bitowego konwertera. Zakres pomiarowy tego urządzenia to &plusmn;200 &mu;V.<br />
====Zadanie:====<br />
Proszę zilustrować efekt kwantyzacji dla trzybitowego przetwornika o zakresie 2.<br />
W tym celu proszę wykonać następujące kroki:<br />
# wygenerować &bdquo;prawie&rdquo; ciągłą sinusoidę (częstość 1, czas trwania 1 próbkowanie co 0,001)<br />
# spróbkować tą sinusoidę co 0,1 (proszę zastosować wycinki)<br />
# proszę skwantować spróbkowane wartości <br />
#Proszę wykreślić na jednym rysunku <br />
#* oryginalny sygnał <br />
#* sygnał spróbkowany w czasie<br />
#* sygnał spróbkowany w czasie i o skwantowanej amplitudzie (skorzystać z funkcji <tt>py.step</tt>)<br />
<br />
{{ Wyjaśnienie|title= wskazówka: Kwantowanie | text = <br />
<source lang =python><br />
N_bits = 3<br />
zakres = 2.0<br />
dy = zakres/2**N_bits<br />
s1_kwantowany = np.floor(s1_samp/dy)*dy +0.5*dy<br />
</source><br />
}}<br />
<br />
<!--<br />
# -*- coding: utf-8 -*-<br />
<br />
import pylab as py<br />
import numpy as np<br />
<br />
t = np.arange(0,1,0.001) # czas 'prawie ciągły'<br />
f1 = 1 # częstość sygnału 1<br />
<br />
s1 = np.sin(2*np.pi*f1*t) # prawie ciągły sygnał o częstości f1<br />
T_samp = t[0::100]<br />
s1_samp = s1[0::100]<br />
<br />
N_bits = 3<br />
zakres = 2.0<br />
dy = zakres/2**N_bits<br />
s1_kwantowany = np.floor(s1_samp/dy)*dy +0.5*dy<br />
ax = py.subplot(111)<br />
py.plot(t,s1)<br />
py.plot(T_samp,s1_samp,'ko',markersize = 5)<br />
py.step(T_samp,s1_kwantowany,where = 'post')<br />
<br />
ax.set_xticks(T_samp)<br />
poziomy = np.arange(-1+0.5*dy,1,dy)<br />
ax.set_yticks(poziomy)<br />
py.grid('on')<br />
<br />
<br />
for i in xrange(0,8):<br />
py.text(1.02,poziomy[i],bin(i))<br />
--><br />
<br />
==Sygnały testowe==<br />
<br />
<br />
=== Generowanie sygnałów testowych ===<br />
<br />
Do badania różnych metod analizy sygnałów potrzebne nam będą sygnały o znanych własnościach. W szczególności dobrze jest umnieć nadać sygnałom występującym w postaci cyfrowej, oraz sztucznym sygnałom próbnym pewne własności fizyczne takie jak:<br />
<br />
<ul><br />
<br />
<li><br />
czętość próbkowania<br />
<br />
<br />
<li><br />
czas trwania<br />
<br />
<br />
<li><br />
amplituda<br />
<br />
</ul><br />
=== Przykład sinus===<br />
Sinus o zadanej częstości (w Hz), długości trwania, częstości próbkowania i fazie.<br />
Poniższy kod implementuje i testuje funkcję <br />
:<math> \sin(f,T,Fs,\phi) = \sin(2*\pi f t)</math> dla <math>t \in \{0,T\}</math> <br />
<br />
<source lang= python><br />
# -*- coding: utf-8 -*-<br />
<br />
import pylab as py<br />
import numpy as np<br />
<br />
def sin(f = 1, T = 1, Fs = 128, phi =0 ):<br />
'''sin o zadanej częstości (w Hz), długości, fazie i częstości próbkowania<br />
Domyślnie wytwarzany jest sygnał reprezentujący <br />
1 sekundę sinusa o częstości 1 Hz i zerowej fazie próbkowanego 128 Hz<br />
'''<br />
<br />
dt = 1.0/Fs<br />
t = np.arange(0,T,dt)<br />
s = np.sin(2*np.pi*f*t + phi)<br />
return (s,t)<br />
<br />
<br />
(s,t) = sin(f=10,Fs=1000)<br />
py.plot(t,s)<br />
py.show()<br />
</source><br />
=== Przykład: eksport sygnału do pliku binarnego ===<br />
*Przypominamy:<br />
** dtype http://docs.scipy.org/doc/numpy/reference/arrays.dtypes.html#arrays-dtypes-constructing<br />
** open/close<br />
** tofile<br />
** fromfile<br />
<br />
* Poniższy kod ilustruje sposób zapisu dwóch funkcji sinus o częstościach 10 Hz i 21 Hz do pliku binarnego:<br />
<br />
<source lang= python><br />
# -*- coding: utf-8 -*-<br />
<br />
import numpy as np<br />
<br />
T = 5<br />
Fs = 128<br />
<br />
(s1,t) = sin(f=10, T=T, Fs=Fs)<br />
(s2,t) = sin(f=21, T=T, Fs=Fs) <br />
<br />
signal = np.zeros((T*Fs, 2), dtype='<f')<br />
signal[:, 0] = s1<br />
signal[:, 1] = s2<br />
<br />
# zapis sygnału binarnego:<br />
f_out = open('test_signal1.bin', 'wb') # otwieramy plik do pisania binarnego: 'wb'<br />
signal.tofile(f_out) # zrzucamy zawartość tablicy ''signal'' do pliku identyfikowanego przez ''f_out''<br />
f_out.close() # zamykamy plik<br />
<br />
# czynność odwrotna - wczytanie sygnału binarnego<br />
ch = 2 # liczba kanałów<br />
fin = open('test_signal1.bin', 'rb') # otwieramy plik do czytania binarnego: 'rb'<br />
s = np.fromfile(fin, dtype='<f') # tworzymy tablicę sig o typie określonym przez ''dtype'' wkładając do niej bity z pliku ''fin'' interpretowane zgodnie z ''dtype''<br />
fin.close() # zamykamy plik<br />
s = np.reshape(s,(len(s)/ch,ch)) # zmieniamy tablicę z jednowymiarowej na dwuwymiarową<br />
</source><br />
<br />
===Przykład: wczytanie sygnału do Svaroga ===<br />
W celu wczytania zapisanego binarnie sygnału do programu Svarog, po wybraniu File -> Open signal, należy wprowadzić częstość próbkowania sygnału oraz liczbę kanałów. <br />
<br />
[[Plik:svarog_open_signal.png|center|800px|thumb|<figure id="uid3" />Wczytywanie sygnału w programie Svarog.]]<br />
<br />
=== Zadanie: Nieznany format danych===<br />
Posiadamy 3 kanałowy plik binarny z sygnałem EKG i EEG. Niestety zapomnieliśmy jaki jest format danych. Proszę wczytać plik i przy pomocy biblioteki matplotlib oraz grafik z fragmentami sygnałów odgadnąć, która z pośród wymienionych zmiennych dtype jest prawdziwa:<br />
<ul><br />
<li> 'float32'<br />
<br />
<li> '>H'<br />
<br />
<li> 'unint64'<br />
<br />
<li> 'float64'<br />
</ul><br />
Plik binarny oraz grafiki można pobrać z [https://drive.google.com/drive/folders/0ByBe6Y9KQQXYRnVVQVhtdnJkWHM?usp=sharing stąd]<br />
<br />
=== Zadanie: Nieznana liczba kanałów===<br />
Posiadamy n kanałowy plik binarny z sygnałem EKG, EEG i EMG. Niestety zapomnieliśmy jaka jest liczba kanałów. Proszę wczytać plik i przy pomocy biblioteki matplotlib oraz grafik z fragmentami sygnałów odgadnąć liczbę kanałów. Plik binarny oraz grafiki można pobrać z [https://drive.google.com/drive/folders/0ByBe6Y9KQQXYVFpOV3A4ZllYUGc?usp=sharing stąd]<br />
<br />
=== Zadanie:===<br />
Proszę stworzyć nową macierz z n+3 kanałami, zawierającą sygnały z dwóch poprzednich zadań. Następnie proszę zapisać macierz sygnałów do pliku binarnego i otworzyć w programie Svarog.<br />
<br />
=== Delta ===<br />
Podobnie można zdefiniować funkcję delta o zadanym czasie trwania, częstości próbkowania i momencie wystąpienia impulsu:<br />
<br />
:<math> \delta(t_0) = \left\{^{1 \quad t=t_0} _{0 \quad t \ne t_0} \right.</math><br />
<br />
<source lang= python><br />
def delta(t0=0.5, T=1 ,Fs = 128):<br />
dt = 1.0/Fs<br />
t = np.arange(0,T,dt)<br />
d = np.zeros(len(t))<br />
d[np.ceil(t0*Fs)]=1<br />
return (d,t)<br />
</source><br />
<br />
=== Zadanie:===<br />
Analogicznie do powyższych przykładów proszę zaimplementować i przetestować funkcje generujące:<br />
* funkcję Gabora (funkcja Gaussa modulowana cosinusem) o zadanej częstości i standardowym odchyleniu w czasie, momencie wystąpienia, długości, częstości próbkowania i fazie. <br />
:<math> g = \exp\left(-\frac{1}{2}\left(\frac{t-t_0}{\sigma}\right)^2 \right) \cdot \cos(2 \pi f t + \phi); </math><br />
<br />
* szum gaussowski o zadanej średniej, odchyleniu standardowym, długości i częstości próbkowania.<br />
<br />
<!-- * pochodną funkcji Gaussa<br />
<br />
* połówkę funkcji Gaussa<br />
--><br />
<br />
==Co musimy z tego zapamiętać? ==<br />
* sygnały dyskretne - dyskretne chwile czasu<br />
* tworzenie sygnałów o konkretnej interpretacji fizycznej (częstość sygnału i częstość próbkowania)<br />
* sygnały w pythonie przechowujemy w tablicach numpy<br />
* zapis sygnałów do pliku i wczytywanie sygnałów z plików binarnych<br />
* wczytywanie sygnałów do SVAROGA<br />
<br />
[[Analiza_sygnałów_-_ćwiczenia]]/Sygnały</div>Prozanskihttp://brain.fuw.edu.pl/edu/index.php?title=%C4%86wiczenia_5&diff=4832Ćwiczenia 52016-03-02T21:23:10Z<p>Prozanski: cofnięcie testu</p>
<hr />
<div>==Model AR==<br />
===Proces AR===<br />
<br />
Proces autoregresyjny rzędu <math>p</math> jest zdefiniowany jako:<br />
::<math>x_t = \sum _{i=1}^p a_i x_{t-i} +\epsilon _t</math><br />
Poniższy kod demonstruje generację pięciu realizacji procesu AR drugiego rzędu, każda o długości 500 punktów:<br />
<source lang = python><br />
import numpy as np<br />
import pylab as py<br />
#wspolczynniki modelu AR<br />
<br />
a = np.array([0.9, -0.6])<br />
<br />
# chcemy badac N punktow procesu<br />
<br />
for j in range(1,6):<br />
N=500<br />
x=np.zeros(N);<br />
#generujemy realizacje procesu<br />
for i in range(2,N):<br />
x[i]=a[0]*x[i-1]+a[1]*x[i-2]+np.random.randn()<br />
py.subplot(5,1,j)<br />
py.plot(x)<br />
py.title('realizacja'+str(j))<br />
py.show()<br />
</source><br />
<br />
===Estymacja parametrów===<br />
<br />
Algorytmów służących do estymacji parametrów modelu AR jest kilka. Tu przedstawimy algorytm Yule-Walker'a: <br />
* mnożymy stronami równania opisujące proces dla póbki <math>t</math> i <math>t-m</math><br />
::<math>x_t x_{t-m} = \sum _{i=1}^p a_i x_{t-i} x_{t-m} +\epsilon _t x_{t-m} </math><br />
* 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:<br />
::<math>R(m) = \sum _{i=1}^p a_i R(m-i)+ \sigma _\epsilon ^2 \delta (m)</math><br />
gdzie <math>m=0,\dots ,p.</math><br />
* Dla <math>m>0</math> możemy zapisać stąd układ równań:<br />
::<math>\left[\begin{array}{c}<br />
R(1)\\<br />
R(2)\\<br />
\vdots \\<br />
R(p)<br />
\end{array}\right]=<br />
\left[\begin{array}{cccc}<br />
R(0)& R(-1) &\dots &\\<br />
R(1)& R(0) &R(-1) \dots &\\<br />
\vdots & & &\\<br />
R(p-1) & &\dots &R(0)<br />
\end{array}\right] \left[\begin{array}{c}<br />
a_1\\<br />
a_2\\<br />
\vdots \\<br />
a_p<br />
\end{array} \right] </math><br />
* stąd wyliczamy współczynniki <math>a</math>, <br />
dla <math>m=0</math> mamy<br />
::<math>R(0) = \sum _{k=1}^p a_k R(-k) + \sigma _\epsilon ^2</math><br />
* można stąd wyliczyć <math>\sigma _\epsilon ^2 </math><br />
<br />
===Ćwiczenie: estymacja parametrów procesu AR===<br />
* Wygeneruj 2000 próbek sygnału z modelu AR o parametrach a = {0.9, -0.6}, epsilon=2<br />
* Oblicz funkcję autokorelacji tego sygnału: <tt>ak = np.correlate(x,x,mode='full')</tt><br />
* Znormalizuj tą funkcję dzieląc ją przez długość sygnału.<br />
* Oblicz parametry zgodnie ze wzorami z poprzedniego paragrafu dla modelu rzędu 2. (wypisz konkretną postać wzorów analitycznie a następnie zaimplementuje je)<br />
wskazówka: <tt>R[0]=ak[N-1]</tt><br />
* Wypisz parametry prawdziwe i estymowane.<br />
* 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)<br />
*<br />
<!--<br />
<source lang = python><br />
import numpy as np<br />
import pylab as py<br />
#wspolczynniki modelu AR<br />
<br />
a = np.array([0.9, -0.6])<br />
epsilon=2<br />
N=2000<br />
x=np.zeros(N);<br />
#generujemy realizacje procesu<br />
for i in range(2,N):<br />
x[i]=a[0]*x[i-1]+a[1]*x[i-2]+epsilon*np.random.randn()<br />
<br />
py.subplot(2,1,1)<br />
py.plot(x)<br />
py.title('wygenerowany sygnal')<br />
py.subplot(2,1,2)<br />
ak = np.correlate(x,x,mode='full')<br />
# ak unbiased<br />
ak=ak/N<br />
#for i in range(1,int(np.floor(len(ak)/2.0))):<br />
# ak[N-1-i]=ak[N-1-i]/(N-i)<br />
# ak[N-1+i]=ak[N-1+i]/(N-i)<br />
py.plot(ak)<br />
py.title('funkcja autokorelacij sygnalu x')<br />
<br />
R=ak[N-1:]<br />
r0=R[0]<br />
r1=R[1]<br />
r2=R[2]<br />
<br />
# estymujemy wspolczynniki modelu na podstawie funkncji autokorelacji<br />
<br />
a2=(r1**2-r0*r2)/(r1**2-r0**2)<br />
a1=r1/r0-r1/r0*a2<br />
s_2=(r0-a1*r1-a2*r2)<br />
<br />
print 'prawdziwe wspolczynniki'<br />
print a[0], a[1], epsilon<br />
print 'estymowane wspolczynniki'<br />
print a1, a2, s_2**0.5<br />
py.show()<br />
<br />
</source><br />
<br />
--><br />
<br />
W przypadku modelu rzędu <math>p</math> estymację parametrów metodą Y-W można zaimplementować np. tak:<br />
<source lang = python><br />
def parametryAR(x,p):<br />
'''funkcja estymująca parametry modelu AR <br />
argumenty:<br />
x- sygnał<br />
p - rząd modelu<br />
f. zwraca:<br />
a - wektor współczynników modelu<br />
epsilon - estymowana wariancja szumu<br />
<br />
funkcja wymaga zaimportowania modułu numpy as np<br />
'''<br />
N = len(x)<br />
ak = np.correlate(x,x,mode='full')<br />
ak=ak/N<br />
R=ak[N-1:]<br />
RL = R[1:1+p]<br />
RP = np.zeros((p,p))<br />
for i in range(p):<br />
aa = ak[N-1-i:N-1-i+p]<br />
RP[i,:] = aa<br />
a=np.linalg.solve(RP,RL)<br />
epsilon = (ak[N-1] - np.sum(a*ak[N:N+p]))**0.5<br />
return a,epsilon<br />
</source><br />
<br />
===Jak znaleźć rząd modelu?===<br />
<br />
Kryterium Akaike (AIC):<br />
<br />
::<math>\mathrm{AIC}(p)= \frac{2p}{N} +\ln(V) </math><br />
<br />
<math>p</math> - ilość parametrów modelu,<br />
<br />
<math>N</math> - ilość próbek sygnału,<br />
<br />
<math>V</math> - wariancja szumu.<br />
<br />
Kryterium to karze za zwiększanie ilości parametrów i nagradza za zmniejszanie niewytłumaczonej wariancji.<br />
<br />
Poniższy kod jest przykładową implementacją kryterium AIC:<br />
<source lang = python><br />
<br />
def kryterium_AIC(x,maxymalnyRzad):<br />
zakres_rzedow = range(1,maxymalnyRzad)<br />
N = len(x)<br />
AIC = np.zeros(len(zakres_rzedow))<br />
for p in zakres_rzedow:<br />
a,epsilon = parametryAR(x,p)<br />
AIC[p-1] = (2.0*p)/N + np.log(np.sqrt(epsilon))<br />
print 'p:', p, ' a:',a,' epsilon: ',epsilon<br />
return AIC<br />
<br />
</source><br />
<br />
Zobaczmy jak działa to na przykładowym syganle AR:<br />
<source lang = python><br />
import numpy as np<br />
import pylab as py<br />
from numpy.fft import fft, fftshift<br />
def parametryAR(x,p):<br />
'''funkcja estymująca parametry modelu AR <br />
argumenty:<br />
x- sygnał<br />
p - rząd modelu<br />
f. zwraca:<br />
a - wektor współczynników modelu<br />
epsilon - estymowana wariancja szumu<br />
<br />
funkcja wymaga zaimportowania modułu numpy as np<br />
'''<br />
N = len(x)<br />
ak = np.correlate(x,x,mode='full')<br />
ak=ak/N<br />
R=ak[N-1:]<br />
RL = R[1:1+p]<br />
RP = np.zeros((p,p))<br />
for i in range(p):<br />
aa = ak[N-1-i:N-1-i+p]<br />
RP[i,:] = aa<br />
a=np.linalg.solve(RP,RL)<br />
epsilon = (ak[N-1] - np.sum(a*ak[N:N+p]))**0.5<br />
return a,epsilon<br />
<br />
<br />
def kryterium_AIC(x,maxymalnyRzad):<br />
zakres_rzedow = range(1,maxymalnyRzad)<br />
AIC = np.zeros(len(zakres_rzedow))<br />
for p in zakres_rzedow:<br />
a,epsilon = parametryAR(x,p)<br />
AIC[p-1] = (2.0*p)/N + np.log(np.sqrt(epsilon))<br />
print 'p:', p, ' a:',a,' epsilon: ',epsilon<br />
return AIC<br />
<br />
#wspolczynniki modelu AR <br />
a = np.array([0.9, -0.7])<br />
epsilon=2<br />
N=600<br />
x=np.zeros(N);<br />
#generujemy realizacje procesu<br />
for i in range(2,N):<br />
x[i]=a[0]*x[i-1]+a[1]*x[i-2] +epsilon*np.random.randn()<br />
<br />
py.subplot(2,1,1)<br />
py.plot(x)<br />
py.title('wygenerowany sygnal')<br />
py.subplot(2,1,2)<br />
<br />
AIC = kryterium_AIC(x,6)<br />
py.plot(range(1,len(AIC)+1),AIC)<br />
py.title('Kryterium AIC')<br />
py.show()<br />
<br />
</source><br />
<br />
== Widmo modelu AR==<br />
<br />
Widmo modelu można wyliczyć analitycznie znając jego współczynniki:<br />
Przepisujemy równanie modelu:<br />
<br />
::<math>x_t = \sum _{i=1}^p a_i x_{t-i} +\epsilon _t</math><br />
<br />
::<math>\sum _{i=0}^p a_i x_{t-i} =\epsilon _t</math><br />
<br />
biorąc transformaty <math>Z</math> obu stron mamy równanie algebraiczne:<br />
<br />
::<math>A(f)X(f) =E(f)</math><br />
<br />
::<math>X(f)=A^{-1}(f) E(f)=H(f) E(f)</math><br />
<br />
Stąd widmo:<br />
<br />
::<math>S(f) = X(f)X^*(f)=H(f)VH^*(f)</math><br />
<br />
=== Jak znaleźć ''A'' &mdash; transformata Z===<br />
<br />
Transformata Z jest dyskretnym odpowiednikiem transformaty Laplace'a:<br />
<br />
::<math>X(z) = Z\lbrace x[n]\rbrace = \sum _{n=0}^\infty {x[n]z^{-n}}</math><br />
<br />
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> [[Ćwiczenia_2#Dyskretna_Transformata_Fouriera_.28DFT.29|porównaj]].<br />
<br />
=== Własności transformaty Z===<br />
<br />
Transformata ta jest liniowa tzn.<br />
<br />
::<math>\mathrm{Z}\lbrace a_1x_1[n] +a_2x_2[n]\rbrace =a_1X_1(z)+a_2X_2(z)</math><br />
<br />
jak ją policzyć od sygnału przesuniętego w czasie to:<br />
<br />
::<math>\mathrm{Z}\lbrace x[n-k]\rbrace = z^{-k}X(z)</math><br />
<br />
dla impulsu:<br />
<br />
::<math>\mathrm{Z}\lbrace \delta [n]\rbrace =1</math><br />
<br />
więc<br />
<br />
::<math>\mathrm{Z}\lbrace \delta [n-n0]\rbrace = z^{-n0} </math><br />
<br />
Stosując tą transfomatę do procesu AR dostajemy:<br />
<br />
::<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><br />
<br />
::<math>=A(z)X(z)</math><br />
<br />
===Widmo procesu===<br />
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>:<br />
<source lang = python><br />
a=[0.9, -0.6]<br />
epsilon = 2<br />
w=np.arange(-np.pi,np.pi,0.1)<br />
z=np.exp(1j*w)<br />
# dla zadanego modelu<br />
A=-1+a[0]*z**(-1)+a[1]*z**(-2);<br />
</source><br />
<br />
Następnie obliczamy odwrotność ''A'' :<br />
<source lang =python><br />
H=1./A<br />
</source><br />
i obliczamy widmo:<br />
<source lang = python><br />
Sp=H*H.conj()*epsilon**2<br />
Sp = Sp.real<br />
</source><br />
Możemy je wykreślić w funkcji częstości <math>\omega</math>.<br />
<source lang = python><br />
py.plot(w,Sp )<br />
</source><br />
<br />
Operacje te możemy uogólnić i zaimplementować jako funkcję do obliczania widma modelu zadanego przez parametry:<br />
<source lang = python><br />
def widmoAR(parametry_a, epsilon, N_punktow):<br />
w = np.linspace(-np.pi,np.pi,N_punktow)<br />
z = np.exp(1j*w)<br />
A = -1 * np.ones(N_punktow) + 1j*np.zeros(N_punktow)<br />
for i in range(len(parametry_a)):<br />
A += parametry_a[i]*z**(-(i+1))<br />
H = 1./A<br />
Sp = H*H.conj()*epsilon**2 # widmo<br />
Sp = Sp.real<br />
return Sp, w<br />
</source><br />
<br />
==== Ćwiczenie ====<br />
Proszę:<br />
* Wygenerować realizację modelu AR <math>a = \{0.6, -0.7, 0.3, -0.25\}, \quad \varepsilon = 2</math><br />
<source lang = python><br />
def generujAR(a, epsilon, N):<br />
x=np.zeros(N)<br />
rzad = len(a)<br />
for i in range(rzad,N):<br />
for p in range(len(a)):<br />
x[i] += a[p]*x[i-(p+1)]<br />
x[i] += epsilon*np.random.randn()<br />
return x<br />
</source><br />
* Obliczyć widmo dla tego modelu<br />
* Wyestymować parametry modelu na podstawie sygału, zakładając, że rząd jest p = 3,4,5,6<br />
* Obliczyć widmo dla wyestymowanego modelu<br />
* Wykreślić widma prawdziwego modelu i modeli estymowanych<br />
<br />
*<br />
<!--<br />
<source lang = python><br />
import numpy as np<br />
import pylab as py<br />
from numpy.fft import fft, fftshift<br />
<br />
def generujAR(a, epsilon, N):<br />
x=np.zeros(N)<br />
rzad = len(a)<br />
for i in range(rzad,N):<br />
for p in range(len(a)):<br />
x[i] += a[p]*x[i-(p+1)]<br />
x[i] += epsilon*np.random.randn()<br />
return x<br />
<br />
def parametryAR(x,p):<br />
'''funkcja estymująca parametry modelu AR <br />
argumenty:<br />
x- sygnał<br />
p - rząd modelu<br />
f. zwraca:<br />
a - wektor współczynników modelu<br />
epsilon - estymowana wariancja szumu<br />
<br />
funkcja wymaga zaimportowania modułu numpy as np<br />
'''<br />
<br />
ak = np.correlate(x,x,mode='full')<br />
N = len(x)<br />
ak = ak/N<br />
R = ak[N-1:]<br />
RL = R[1:1+p]<br />
RP = np.zeros((p,p))<br />
for i in range(p):<br />
aa = ak[N-1-i:N-1-i+p]<br />
RP[i,:] = aa<br />
a = np.linalg.solve(RP,RL)<br />
epsilon = (ak[N-1] - np.sum(a*ak[N:N+p]))**0.5<br />
return a,epsilon<br />
<br />
def widmoAR(parametry_a, epsilon, N_punktow):<br />
w = np.linspace(-np.pi,np.pi,N_punktow)<br />
z = np.exp(1j*w)<br />
A = -1 * np.ones(N_punktow) + 1j*np.zeros(N_punktow)<br />
for i in range(len(parametry_a)):<br />
A += parametry_a[i]*z**(-(i+1))<br />
H = 1./A<br />
Sp = H*H.conj()*epsilon**2 # widmo<br />
Sp = Sp.real<br />
return Sp, w<br />
<br />
#wspolczynniki modelu AR <br />
a = np.array([0.6, -0.7, 0.3, -0.25])<br />
epsilon = 2<br />
# obliczanie widma z modelu<br />
Sp, w = widmoAR(a,epsilon,200)<br />
<br />
#generujemy realizacje procesu<br />
N=600<br />
x = generujAR(a, epsilon, N)<br />
<br />
# estymujemy wspolczynniki modelu metodą Ypula-Walkera<br />
# obliczamy widmo dla estymowanego modelu<br />
for p in range(3,7):<br />
a_est,epsilon_est = parametryAR(x,p)<br />
Sp_est, w = widmoAR(a_est,epsilon_est,200)<br />
py.plot(w,Sp_est )<br />
py.plot(w,Sp)<br />
py.legend(('p = 3','p = 4','p = 5','p = 6','prawdziwy'))<br />
py.title('widmo z modelu')<br />
<br />
py.show()<br />
</source><br />
--><br />
<br />
==== Ćwiczenie ====<br />
<br />
Dla modelu z poprzedniego ćwiczenia proszę wygenerować realizację sygnału długości 1000 punktów. Proszę porównać widma:<br />
* prawdziwe, obliczone z prawdziwych parametrów modelu<br />
* obliczone z estymowanego modelu<br />
* obliczone przez periodogram<br />
* obliczone metodą Welcha <br />
* obliczone metodą wielookienkową Thomsona<br />
<br />
<source lang = python><br />
import numpy as np<br />
import pylab as py<br />
from numpy.fft import fft, fftshift,fftfreq<br />
import gendpss as dpss<br />
<br />
def generujAR(a, epsilon, N):<br />
x=np.zeros(N)<br />
r = len(a)<br />
for i in range(r,N):<br />
for p in range(len(a)):<br />
x[i] += a[p]*x[i-(p+1)]<br />
x[i] += epsilon*np.random.randn()<br />
return x<br />
<br />
def parametryAR(x,p):<br />
N = len(x)<br />
ak = np.correlate(x,x,mode='full')<br />
ak = ak/N<br />
R = ak[N-1:]<br />
RL = R[1:1+p]<br />
RP = np.zeros((p,p))<br />
for i in range(p):<br />
aa = ak[N-1-i:N-1-i+p]<br />
RP[i,:] = aa<br />
a = np.linalg.solve(RP,RL)<br />
epsilon = (ak[N-1] - np.sum(a*ak[N:N+p]))**0.5<br />
return a,epsilon<br />
<br />
def widmoAR(parametry_a, epsilon, N_punktow):<br />
w = np.linspace(-np.pi,np.pi,N_punktow)<br />
z = np.exp(1j*w)<br />
A = -1 * np.ones(N_punktow) + 1j*np.zeros(N_punktow)<br />
for i in range(len(parametry_a)):<br />
A += parametry_a[i]*z**(-(i+1))<br />
H = 1./A<br />
Sp = H*H.conj()*epsilon**2 # widmo<br />
Sp = Sp.real<br />
return Sp, w<br />
<br />
def periodogram(s, okno , Fs):<br />
s = s*okno<br />
N_fft = len(s)<br />
S = fft(s,N_fft)#/np.sqrt(N_fft)<br />
P = S*S.conj()/np.sum(okno**2)<br />
P = P.real # P i tak ma zerowe wartośći urojone, ale trzeba ykonać konwersję typów<br />
F = fftfreq(N_fft, 1/Fs)<br />
return (fftshift(P),fftshift(F))<br />
<br />
def pwelch(s,okno, przesuniencie, Fs):<br />
N = len(s)<br />
N_s = len(okno)<br />
<br />
start_fragmentow = np.arange(0,N-N_s+1,przesuniencie)<br />
ile_fragmentow = len(start_fragmentow)<br />
ile_przekrycia = N_s*ile_fragmentow/float(N)<br />
print ile_przekrycia, ile_fragmentow<br />
P_sredni = np.zeros(N_s)<br />
for i in range(ile_fragmentow):<br />
s_fragment = s[start_fragmentow[i]:start_fragmentow[i]+N_s]<br />
(P, F) = periodogram(s_fragment,okno,Fs)<br />
P_sredni += P<br />
return (P_sredni/ile_przekrycia,F)#(P_sredni/ile_przekrycia,F)<br />
<br />
def mtm(s, NW = 3, Fs = 128):<br />
K = int(2*NW-1)<br />
N = len(s)<br />
w = dpss.gendpss(N,NW,K)<br />
S=np.zeros(N)<br />
for i in range(K):<br />
Si = np.abs(fft(s*w.dpssarray[i]))**2<br />
S[:] += Si.real<br />
S = S/K<br />
F = fftfreq(N,1.0/Fs)<br />
return (fftshift(S),fftshift(F))<br />
<br />
<br />
#wspolczynniki modelu AR <br />
a = np.array([0.3, 0.2, 0.5, -0.25 ,-0.3])<br />
epsilon = 2<br />
N=256<br />
# obliczanie widma z modelu<br />
Sp, w = widmoAR(a,epsilon,N)<br />
<br />
#generujemy realizacje procesu<br />
<br />
x = generujAR(a, epsilon, N)<br />
<br />
# estymujemy wspolczynniki modelu metodą Yula-Walkera<br />
# obliczamy widmo dla estymowanego modelu<br />
a_est,epsilon_est = parametryAR(x,5)<br />
Sp_est, w = widmoAR(a_est,epsilon_est,N)<br />
<br />
okno = np.blackman(N)<br />
Fs = 2*np.pi <br />
P_periodogram,F_periodogram = periodogram(x, okno , Fs=Fs)<br />
okno = np.blackman(N/4)<br />
P_welch, F_welch = pwelch(x,okno, len(okno)/2, Fs=Fs) <br />
P_mtm, F_mtm = mtm(x, NW = 4.5, Fs =Fs)<br />
<br />
py.plot(w,Sp)<br />
py.plot(w,Sp_est)<br />
py.plot(F_periodogram,P_periodogram)<br />
py.plot(F_welch, P_welch/(len(P_periodogram)/len(P_welch)))#uwaga Welch ma inne df<br />
py.plot(F_mtm, P_mtm)<br />
<br />
#py.legend(('prawdziwy','estymowany z AR','periodogram','Welch','mtm'))<br />
<br />
print 'enenrgia sygnału: ', np.sum(x**2)<br />
print 'enenrgia spektrum AR',np.sum(Sp)<br />
print 'enenrgia est',np.sum(Sp_est)<br />
print 'enenrgia mtm',np.sum(P_mtm)<br />
print 'enenrgia welch',np.sum(P_welch)<br />
print 'enenrgia period',np.sum(P_periodogram)<br />
py.show()<br />
<br />
<br />
</source><br />
<br />
====Ogólny schemat ====<br />
W praktyce musimy najczęściej estymować widmo mając dany tylko sygnał. Wówczas powinniśmy pwstąpić według następującego algorytmu:<br />
* oszacować rząd modelu np. przy pomocy kryterium Akaikego<br />
* wyestymować parametry modelu<br />
* obliczyć widmo dla estymowanego modelu<br />
<br />
====Ćwiczenie====<br />
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ą<br />
[http://brain.fuw.edu.pl/edu/TI:Programowanie_z_Pythonem/Wejście_i_wyjście#Przyk.C5.82ady tu].<br />
Proszę zapisać sygnał c4spin.txt w swoim bieżącym katalogu, wczytać go i oszacować widmo metodą AR i metodą Welcha.<br />
<br />
==Wielokanałowe modele AR==<br />
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ć<br />
wzór modelu ''k''-kanałowego jako<br />
::<math>X(t) = \sum _{i=1}^p A(i) X(t-i) +E(t)</math><br />
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'':<br />
::<math>X(t)=\left(\begin{array}{l}<br />
x_1(t)\\x_2(t)\\\vdots\\x_k(t)<br />
\end{array} \right)\ \ \ E(t)=\left(\begin{array}{l}<br />
\epsilon_1(t)\\<br />
\epsilon_2(t)\\\vdots\\<br />
\epsilon_k(t)<br />
\end{array} \right)</math><br />
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.<br />
<br />
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'', natomiast ''X''(''f'') i ''E''(''f'') są wektorami o rozmiarze ''k'':<br />
::<math>A(f)X(f) =E(f)</math><br />
<br />
::<math>X(f)=A^{-1}(f) E(f)=H(f) E(f)</math><br />
<br />
::<math>S(f) = X(f)X^+(f)=H(f)VH^+(f)</math><br />
(symbol <sup>+</sup> oznacza transpozycję połączoną ze sprzężeniem zespolonym elementów macierzy).<br />
<br />
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''.<br />
<br />
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. <br />
<br />
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.<br />
::<math>\gamma_{ij}(f)=\frac{\left|H_{ij}(f)\right|^2}{\sum_{m=1}^{k}\left|H_{im}(f)\right|^2}</math><br />
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 <br />
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.<br />
<br />
===Ćwiczenia===<br />
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.<br />
<br />
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 układ 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 zaś szumem 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?<br />
*<br />
<!--<br />
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).<br />
--><br />
<br />
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ż układ 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ę. <br />
<br />
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?<br />
<br />
[http://eeg.pl/DTF]</div>Prozanskihttp://brain.fuw.edu.pl/edu/index.php?title=%C4%86wiczenia_5&diff=4831Ćwiczenia 52016-03-02T21:21:57Z<p>Prozanski: test</p>
<hr />
<div>==Model AR==<br />
===Proces AR===<br />
<br />
Proces autoregresyjny rzędu <math>p</math> jest zdefiniowany jako:<br />
::<math>x_t = \sum _{i=1}^p a_i x_{t-i} +\epsilon _t</math><br />
Poniższy kod demonstruje generację pięciu realizacji procesu AR drugiego rzędu, każda o długości 500 punktów:<br />
<source lang = python><br />
import numpy as np<br />
import pylab as py<br />
#wspolczynniki modelu AR<br />
<br />
a = np.array([0.9, -0.6])<br />
<br />
# chcemy badac N punktow procesu<br />
<br />
for j in range(1,6):<br />
N=500<br />
x=np.zeros(N);<br />
#generujemy realizacje procesu<br />
for i in range(2,N):<br />
x[i]=a[0]*x[i-1]+a[1]*x[i-2]+np.random.randn()<br />
py.subplot(5,1,j)<br />
py.plot(x)<br />
py.title('realizacja'+str(j))<br />
py.show()<br />
</source><br />
<br />
===Estymacja parametrów===<br />
<br />
Algorytmów służących do estymacji parametrów modelu AR jest kilka. Tu przedstawimy algorytm Yule-Walker'a: <br />
* mnożymy stronami równania opisujące proces dla póbki <math>t</math> i <math>t-m</math><br />
::<math>x_t x_{t-m} = \sum _{i=1}^p a_i x_{t-i} x_{t-m} +\epsilon _t x_{t-m} </math><br />
* 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:<br />
::<math>R(m) = \sum _{i=1}^p a_i R(m-i)+ \sigma _\epsilon ^2 \delta (m)</math><br />
gdzie <math>m=0,\dots ,p.</math><br />
* Dla <math>m>0</math> możemy zapisać stąd układ równań:<br />
::<math>\left[\begin{array}{c}<br />
R(1)\\<br />
R(2)\\<br />
\vdots \\<br />
R(p)<br />
\end{array}\right]=<br />
\left[\begin{array}{cccc}<br />
R(0)& R(-1) &\dots &\\<br />
R(1)& R(0) &R(-1) \dots &\\<br />
\vdots & & &\\<br />
R(p-1) & &\dots &R(0)<br />
\end{array}\right] \left[\begin{array}{c}<br />
a_1\\<br />
a_2\\<br />
\vdots \\<br />
a_p<br />
\end{array} \right] </math><br />
* stąd wyliczamy współczynniki <math>a</math>, <br />
dla <math>m=0</math> mamy<br />
::<math>R(0) = \sum _{k=1}^p a_k R(-k) + \sigma _\epsilon ^2</math><br />
* można stąd wyliczyć <math>\sigma _\epsilon ^2 </math><br />
<br />
===Ćwiczenie: estymacja parametrów procesu AR===<br />
* Wygeneruj 2000 próbek sygnału z modelu AR o parametrach a = {0.9, -0.6}, epsilon=2<br />
* Oblicz funkcję autokorelacji tego sygnału: <tt>ak = np.correlate(x,x,mode='full')</tt><br />
* Znormalizuj tą funkcję dzieląc ją przez długość sygnału.<br />
* Oblicz parametry zgodnie ze wzorami z poprzedniego paragrafu dla modelu rzędu 2. (wypisz konkretną postać wzorów analitycznie a następnie zaimplementuje je)<br />
wskazówka: <tt>R[0]=ak[N-1]</tt><br />
* Wypisz parametry prawdziwe i estymowane.<br />
* 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)<br />
*<br />
<!--<br />
<source lang = python><br />
import numpy as np<br />
import pylab as py<br />
#wspolczynniki modelu AR<br />
<br />
a = np.array([0.9, -0.6])<br />
epsilon=2<br />
N=2000<br />
x=np.zeros(N);<br />
#generujemy realizacje procesu<br />
for i in range(2,N):<br />
x[i]=a[0]*x[i-1]+a[1]*x[i-2]+epsilon*np.random.randn()<br />
<br />
py.subplot(2,1,1)<br />
py.plot(x)<br />
py.title('wygenerowany sygnal')<br />
py.subplot(2,1,2)<br />
ak = np.correlate(x,x,mode='full')<br />
# ak unbiased<br />
ak=ak/N<br />
#for i in range(1,int(np.floor(len(ak)/2.0))):<br />
# ak[N-1-i]=ak[N-1-i]/(N-i)<br />
# ak[N-1+i]=ak[N-1+i]/(N-i)<br />
py.plot(ak)<br />
py.title('funkcja autokorelacij sygnalu x')<br />
<br />
R=ak[N-1:]<br />
r0=R[0]<br />
r1=R[1]<br />
r2=R[2]<br />
<br />
# estymujemy wspolczynniki modelu na podstawie funkncji autokorelacji<br />
<br />
a2=(r1**2-r0*r2)/(r1**2-r0**2)<br />
a1=r1/r0-r1/r0*a2<br />
s_2=(r0-a1*r1-a2*r2)<br />
<br />
print 'prawdziwe wspolczynniki'<br />
print a[0], a[1], epsilon<br />
print 'estymowane wspolczynniki'<br />
print a1, a2, s_2**0.5<br />
py.show()<br />
<br />
</source><br />
<br />
--><br />
<br />
W przypadku modelu rzędu <math>p</math> estymację parametrów metodą Y-W można zaimplementować np. tak:<br />
<source lang = python><br />
def parametryAR(x,p):<br />
'''funkcja estymująca parametry modelu AR <br />
argumenty:<br />
x- sygnał<br />
p - rząd modelu<br />
f. zwraca:<br />
a - wektor współczynników modelu<br />
epsilon - estymowana wariancja szumu<br />
<br />
funkcja wymaga zaimportowania modułu numpy as np<br />
'''<br />
N = len(x)<br />
ak = np.correlate(x,x,mode='full')<br />
ak=ak/N<br />
R=ak[N-1:]<br />
RL = R[1:1+p]<br />
RP = np.zeros((p,p))<br />
for i in range(p):<br />
aa = ak[N-1-i:N-1-i+p]<br />
RP[i,:] = aa<br />
a=np.linalg.solve(RP,RL)<br />
epsilon = (ak[N-1] - np.sum(a*ak[N:N+p]))**0.5<br />
return a,epsilon<br />
</source><br />
<br />
===Jak znaleźć rząd modelu?===<br />
<br />
Kryterium Akaike (AIC):<br />
<br />
::<math>\mathrm{AIC}(p)= \frac{2p}{N} +\ln(V) </math><br />
<br />
<math>p</math> - ilość parametrów modelu,<br />
<br />
<math>N</math> - ilość próbek sygnału,<br />
<br />
<math>V</math> - wariancja szumu.<br />
<br />
Kryterium to karze za zwiększanie ilości parametrów i nagradza za zmniejszanie niewytłumaczonej wariancji.<br />
<br />
Poniższy kod jest przykładową implementacją kryterium AIC:<br />
<source lang = python><br />
<br />
def kryterium_AIC(x,maxymalnyRzad):<br />
zakres_rzedow = range(1,maxymalnyRzad)<br />
N = len(x)<br />
AIC = np.zeros(len(zakres_rzedow))<br />
for p in zakres_rzedow:<br />
a,epsilon = parametryAR(x,p)<br />
AIC[p-1] = (2.0*p)/N + np.log(np.sqrt(epsilon))<br />
print 'p:', p, ' a:',a,' epsilon: ',epsilon<br />
return AIC<br />
<br />
</source><br />
<br />
Zobaczmy jak działa to na przykładowym syganle AR:<br />
<source lang = python><br />
import numpy as np<br />
import pylab as py<br />
from numpy.fft import fft, fftshift<br />
def parametryAR(x,p):<br />
'''funkcja estymująca parametry modelu AR <br />
argumenty:<br />
x- sygnał<br />
p - rząd modelu<br />
f. zwraca:<br />
a - wektor współczynników modelu<br />
epsilon - estymowana wariancja szumu<br />
<br />
funkcja wymaga zaimportowania modułu numpy as np<br />
'''<br />
N = len(x)<br />
ak = np.correlate(x,x,mode='full')<br />
ak=ak/N<br />
R=ak[N-1:]<br />
RL = R[1:1+p]<br />
RP = np.zeros((p,p))<br />
for i in range(p):<br />
aa = ak[N-1-i:N-1-i+p]<br />
RP[i,:] = aa<br />
a=np.linalg.solve(RP,RL)<br />
epsilon = (ak[N-1] - np.sum(a*ak[N:N+p]))**0.5<br />
return a,epsilon<br />
<br />
<br />
def kryterium_AIC(x,maxymalnyRzad):<br />
zakres_rzedow = range(1,maxymalnyRzad)<br />
AIC = np.zeros(len(zakres_rzedow))<br />
for p in zakres_rzedow:<br />
a,epsilon = parametryAR(x,p)<br />
AIC[p-1] = (2.0*p)/N + np.log(np.sqrt(epsilon))<br />
print 'p:', p, ' a:',a,' epsilon: ',epsilon<br />
return AIC<br />
<br />
#wspolczynniki modelu AR <br />
a = np.array([0.9, -0.7])<br />
epsilon=2<br />
N=600<br />
x=np.zeros(N);<br />
#generujemy realizacje procesu<br />
for i in range(2,N):<br />
x[i]=a[0]*x[i-1]+a[1]*x[i-2] +epsilon*np.random.randn()<br />
<br />
py.subplot(2,1,1)<br />
py.plot(x)<br />
py.title('wygenerowany sygnal')<br />
py.subplot(2,1,2)<br />
<br />
AIC = kryterium_AIC(x,6)<br />
py.plot(range(1,len(AIC)+1),AIC)<br />
py.title('Kryterium AIC')<br />
py.show()<br />
<br />
</source><br />
<br />
== Widmo modelu AR==<br />
<br />
Widmo modelu można wyliczyć analitycznie znając jego współczynniki:<br />
Przepisujemy równanie modelu:<br />
<br />
::<math>x_t = \sum _{i=1}^p a_i x_{t-i} +\epsilon _t</math><br />
<br />
::<math>\sum _{i=0}^p a_i x_{t-i} =\epsilon _t</math><br />
<br />
biorąc transformaty <math>Z</math> obu stron mamy równanie algebraiczne:<br />
<br />
::<math>A(f)X(f) =E(f)</math><br />
<br />
::<math>X(f)=A^{-1}(f) E(f)=H(f) E(f)</math><br />
<br />
Stąd widmo:<br />
<br />
::<math>S(f) = X(f)X^*(f)=H(f)VH^*(f)</math><br />
<br />
=== Jak znaleźć ''A'' &mdash; transformata Z===<br />
<br />
Transformata Z jest dyskretnym odpowiednikiem transformaty Laplace'a:<br />
<br />
::<math>X(z) = Z\lbrace x[n]\rbrace = \sum _{n=0}^\infty {x[n]z^{-n}}</math><br />
<br />
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> [[Ćwiczenia_2#Dyskretna_Transformata_Fouriera_.28DFT.29|porównaj]].<br />
<br />
=== Własności transformaty Z===<br />
<br />
Transformata ta jest liniowa tzn.<br />
<br />
::<math>\mathrm{Z}\lbrace a_1x_1[n] +a_2x_2[n]\rbrace =a_1X_1(z)+a_2X_2(z)</math><br />
<br />
jak ją policzyć od sygnału przesuniętego w czasie to:<br />
<br />
::<math>\mathrm{Z}\lbrace x[n-k]\rbrace = z^{-k}X(z)</math><br />
<br />
dla impulsu:<br />
<br />
::<math>\mathrm{Z}\lbrace \delta [n]\rbrace =1</math><br />
<br />
więc<br />
<br />
::<math>\mathrm{Z}\lbrace \delta [n-n0]\rbrace = z^{-n0} </math><br />
<br />
Stosując tą transfomatę do procesu AR dostajemy:<br />
<br />
::<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><br />
<br />
::<math>=A(z)X(z)</math><br />
<br />
===Widmo procesu===<br />
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>:<br />
<source lang = python><br />
a=[0.9, -0.6]<br />
epsilon = 2<br />
w=np.arange(-np.pi,np.pi,0.1)<br />
z=np.exp(1j*w)<br />
# dla zadanego modelu<br />
A=-1+a[0]*z**(-1)+a[1]*z**(-2);<br />
</source><br />
<br />
Następnie obliczamy odwrotność ''A'' :<br />
<source lang =python><br />
H=1./A<br />
</source><br />
i obliczamy widmo:<br />
<source lang = python><br />
Sp=H*H.conj()*epsilon**2<br />
Sp = Sp.real<br />
</source><br />
Możemy je wykreślić w funkcji częstości <math>\omega</math>.<br />
<source lang = python><br />
py.plot(w,Sp )<br />
</source><br />
<br />
Operacje te możemy uogólnić i zaimplementować jako funkcję do obliczania widma modelu zadanego przez parametry:<br />
<source lang = python><br />
def widmoAR(parametry_a, epsilon, N_punktow):<br />
w = np.linspace(-np.pi,np.pi,N_punktow)<br />
z = np.exp(1j*w)<br />
A = -1 * np.ones(N_punktow) + 1j*np.zeros(N_punktow)<br />
for i in range(len(parametry_a)):<br />
A += parametry_a[i]*z**(-(i+1))<br />
H = 1./A<br />
Sp = H*H.conj()*epsilon**2 # widmo<br />
Sp = Sp.real<br />
return Sp, w<br />
</source><br />
<br />
==== Ćwiczenie ====<br />
Proszę:<br />
* Wygenerować realizację modelu AR <math>a = \{0.6, -0.7, 0.3, -0.25\}, \quad \varepsilon = 2</math><br />
<source lang = python><br />
def generujAR(a, epsilon, N):<br />
x=np.zeros(N)<br />
rzad = len(a)<br />
for i in range(rzad,N):<br />
for p in range(len(a)):<br />
x[i] += a[p]*x[i-(p+1)]<br />
x[i] += epsilon*np.random.randn()<br />
return x<br />
</source><br />
* Obliczyć widmo dla tego modelu<br />
* Wyestymować parametry modelu na podstawie sygału, zakładając, że rząd jest p = 3,4,5,6<br />
* Obliczyć widmo dla wyestymowanego modelu<br />
* Wykreślić widma prawdziwego modelu i modeli estymowanych<br />
<br />
*<br />
<!--<br />
<source lang = python><br />
import numpy as np<br />
import pylab as py<br />
from numpy.fft import fft, fftshift<br />
<br />
def generujAR(a, epsilon, N):<br />
x=np.zeros(N)<br />
rzad = len(a)<br />
for i in range(rzad,N):<br />
for p in range(len(a)):<br />
x[i] += a[p]*x[i-(p+1)]<br />
x[i] += epsilon*np.random.randn()<br />
return x<br />
<br />
def parametryAR(x,p):<br />
'''funkcja estymująca parametry modelu AR <br />
argumenty:<br />
x- sygnał<br />
p - rząd modelu<br />
f. zwraca:<br />
a - wektor współczynników modelu<br />
epsilon - estymowana wariancja szumu<br />
<br />
funkcja wymaga zaimportowania modułu numpy as np<br />
'''<br />
<br />
ak = np.correlate(x,x,mode='full')<br />
N = len(x)<br />
ak = ak/N<br />
R = ak[N-1:]<br />
RL = R[1:1+p]<br />
RP = np.zeros((p,p))<br />
for i in range(p):<br />
aa = ak[N-1-i:N-1-i+p]<br />
RP[i,:] = aa<br />
a = np.linalg.solve(RP,RL)<br />
epsilon = (ak[N-1] - np.sum(a*ak[N:N+p]))**0.5<br />
return a,epsilon<br />
<br />
def widmoAR(parametry_a, epsilon, N_punktow):<br />
w = np.linspace(-np.pi,np.pi,N_punktow)<br />
z = np.exp(1j*w)<br />
A = -1 * np.ones(N_punktow) + 1j*np.zeros(N_punktow)<br />
for i in range(len(parametry_a)):<br />
A += parametry_a[i]*z**(-(i+1))<br />
H = 1./A<br />
Sp = H*H.conj()*epsilon**2 # widmo<br />
Sp = Sp.real<br />
return Sp, w<br />
<br />
#wspolczynniki modelu AR <br />
a = np.array([0.6, -0.7, 0.3, -0.25])<br />
epsilon = 2<br />
# obliczanie widma z modelu<br />
Sp, w = widmoAR(a,epsilon,200)<br />
<br />
#generujemy realizacje procesu<br />
N=600<br />
x = generujAR(a, epsilon, N)<br />
<br />
# estymujemy wspolczynniki modelu metodą Ypula-Walkera<br />
# obliczamy widmo dla estymowanego modelu<br />
for p in range(3,7):<br />
a_est,epsilon_est = parametryAR(x,p)<br />
Sp_est, w = widmoAR(a_est,epsilon_est,200)<br />
py.plot(w,Sp_est )<br />
py.plot(w,Sp)<br />
py.legend(('p = 3','p = 4','p = 5','p = 6','prawdziwy'))<br />
py.title('widmo z modelu')<br />
<br />
py.show()<br />
</source><br />
--><br />
<br />
==== Ćwiczenie ====<br />
<br />
Dla modelu z poprzedniego ćwiczenia proszę wygenerować realizację sygnału długości 1000 punktów. Proszę porównać widma:<br />
* prawdziwe, obliczone z prawdziwych parametrów modelu<br />
* obliczone z estymowanego modelu<br />
* obliczone przez periodogram<br />
* obliczone metodą Welcha <br />
* obliczone metodą wielookienkową Thomsona<br />
<br />
<source lang = python><br />
import numpy as np<br />
import pylab as py<br />
from numpy.fft import fft, fftshift,fftfreq<br />
import gendpss as dpss<br />
<br />
def generujAR(a, epsilon, N):<br />
x=np.zeros(N)<br />
r = len(a)<br />
for i in range(r,N):<br />
for p in range(len(a)):<br />
x[i] += a[p]*x[i-(p+1)]<br />
x[i] += epsilon*np.random.randn()<br />
return x<br />
<br />
def parametryAR(x,p):<br />
N = len(x)<br />
ak = np.correlate(x,x,mode='full')<br />
ak = ak/N<br />
R = ak[N-1:]<br />
RL = R[1:1+p]<br />
RP = np.zeros((p,p))<br />
for i in range(p):<br />
aa = ak[N-1-i:N-1-i+p]<br />
RP[i,:] = aa<br />
a = np.linalg.solve(RP,RL)<br />
epsilon = (ak[N-1] - np.sum(a*ak[N:N+p]))**0.5<br />
return a,epsilon<br />
<br />
def widmoAR(parametry_a, epsilon, N_punktow):<br />
w = np.linspace(-np.pi,np.pi,N_punktow)<br />
z = np.exp(1j*w)<br />
A = -1 * np.ones(N_punktow) + 1j*np.zeros(N_punktow)<br />
for i in range(len(parametry_a)):<br />
A += parametry_a[i]*z**(-(i+1))<br />
H = 1./A<br />
Sp = H*H.conj()*epsilon**2 # widmo<br />
Sp = Sp.real<br />
return Sp, w<br />
<br />
def periodogram(s, okno , Fs):<br />
s = s*okno<br />
N_fft = len(s)<br />
S = fft(s,N_fft)#/np.sqrt(N_fft)<br />
P = S*S.conj()/np.sum(okno**2)<br />
P = P.real # P i tak ma zerowe wartośći urojone, ale trzeba ykonać konwersję typów<br />
F = fftfreq(N_fft, 1/Fs)<br />
return (fftshift(P),fftshift(F))<br />
<br />
def pwelch(s,okno, przesuniencie, Fs):<br />
N = len(s)<br />
N_s = len(okno)<br />
<br />
start_fragmentow = np.arange(0,N-N_s+1,przesuniencie)<br />
ile_fragmentow = len(start_fragmentow)<br />
ile_przekrycia = N_s*ile_fragmentow/float(N)<br />
print ile_przekrycia, ile_fragmentow<br />
P_sredni = np.zeros(N_s)<br />
for i in range(ile_fragmentow):<br />
s_fragment = s[start_fragmentow[i]:start_fragmentow[i]+N_s]<br />
(P, F) = periodogram(s_fragment,okno,Fs)<br />
P_sredni += P<br />
return (P_sredni/ile_przekrycia,F)#(P_sredni/ile_przekrycia,F)<br />
<br />
def mtm(s, NW = 3, Fs = 128):<br />
K = int(2*NW-1)<br />
N = len(s)<br />
w = dpss.gendpss(N,NW,K)<br />
S=np.zeros(N)<br />
for i in range(K):<br />
Si = np.abs(fft(s*w.dpssarray[i]))**2<br />
S[:] += Si.real<br />
S = S/K<br />
F = fftfreq(N,1.0/Fs)<br />
return (fftshift(S),fftshift(F))<br />
<br />
<br />
#wspolczynniki modelu AR <br />
a = np.array([0.3, 0.2, 0.5, -0.25 ,-0.3])<br />
epsilon = 2<br />
N=256<br />
# obliczanie widma z modelu<br />
Sp, w = widmoAR(a,epsilon,N)<br />
<br />
#generujemy realizacje procesu<br />
<br />
x = generujAR(a, epsilon, N)<br />
<br />
# estymujemy wspolczynniki modelu metodą Yula-Walkera<br />
# obliczamy widmo dla estymowanego modelu<br />
a_est,epsilon_est = parametryAR(x,5)<br />
Sp_est, w = widmoAR(a_est,epsilon_est,N)<br />
<br />
okno = np.blackman(N)<br />
Fs = 2*np.pi <br />
P_periodogram,F_periodogram = periodogram(x, okno , Fs=Fs)<br />
okno = np.blackman(N/4)<br />
P_welch, F_welch = pwelch(x,okno, len(okno)/2, Fs=Fs) <br />
P_mtm, F_mtm = mtm(x, NW = 4.5, Fs =Fs)<br />
<br />
py.plot(w,Sp)<br />
py.plot(w,Sp_est)<br />
py.plot(F_periodogram,P_periodogram)<br />
py.plot(F_welch, P_welch/(len(P_periodogram)/len(P_welch)))#uwaga Welch ma inne df<br />
py.plot(F_mtm, P_mtm)<br />
<br />
#py.legend(('prawdziwy','estymowany z AR','periodogram','Welch','mtm'))<br />
<br />
print 'enenrgia sygnału: ', np.sum(x**2)<br />
print 'enenrgia spektrum AR',np.sum(Sp)<br />
print 'enenrgia est',np.sum(Sp_est)<br />
print 'enenrgia mtm',np.sum(P_mtm)<br />
print 'enenrgia welch',np.sum(P_welch)<br />
print 'enenrgia period',np.sum(P_periodogram)<br />
py.show()<br />
<br />
<br />
</source><br />
<br />
====Ogólny schemat ====<br />
W praktyce musimy najczęściej estymować widmo mając dany tylko sygnał. Wówczas powinniśmy pwstąpić według następującego algorytmu:<br />
* oszacować rząd modelu np. przy pomocy kryterium Akaikego<br />
* wyestymować parametry modelu<br />
* obliczyć widmo dla estymowanego modelu<br />
<br />
====Ćwiczenie====<br />
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ą<br />
[http://brain.fuw.edu.pl/edu/TI:Programowanie_z_Pythonem/Wejście_i_wyjście#Przyk.C5.82ady tu].<br />
Proszę zapisać sygnał c4spin.txt w swoim bieżącym katalogu, wczytać go i oszacować widmo metodą AR i metodą Welcha.<br />
<br />
==Wielokanałowe modele AR==<br />
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ć<br />
wzór modelu ''k''-kanałowego jako<br />
::<math>X(t) = \sum _{i=1}^p A(i) X(t-i) +E(t)</math><br />
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'':<br />
::<math>X(t)=\left(\begin{array}{l}<br />
x_1(t)\\x_2(t)\\\vdots\\x_k(t)<br />
\end{array} \right)\ \ \ E(t)=\left(\begin{array}{l}<br />
\epsilon_1(t)\\<br />
\epsilon_2(t)\\\vdots\\<br />
\epsilon_k(t)<br />
\end{array} \right)</math><br />
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.<br />
<br />
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'', natomiast ''X''(''f'') i ''E''(''f'') są wektorami o rozmiarze ''k'':<br />
::<math>A(f)X(f) =E(f)</math><br />
<br />
::<math>X(f)=A^{-1}(f) E(f)=H(f) E(f)</math><br />
<br />
::<math>S(f) = X(f)X^+(f)=H(f)VH^+(f)</math><br />
(symbol <sup>+</sup> oznacza transpozycję połączoną ze sprzężeniem zespolonym elementów macierzy).<br />
<br />
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''.<br />
<br />
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. <br />
<br />
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.<br />
::<math>\gamma_{ij}(f)=\frac{\left|H_{ij}(f)\right|^2}{\sum_{m=1}^{k}\left|H_{im}(f)\right|^2}</math><br />
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 <br />
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.<br />
<br />
===Ćwiczenia===<br />
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.<br />
<br />
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 układ 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 zaś szumem 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?<br />
*<br />
<!--<br />
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).<br />
--><br />
<br />
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ż układ 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ę. <br />
<br />
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?<br />
<br />
[http://eeg.pl/DTF]<br />
[[Plik:dane_dtf_2.bin]]</div>Prozanskihttp://brain.fuw.edu.pl/edu/index.php?title=%C4%86wiczenia_1&diff=4743Ćwiczenia 12016-02-18T01:48:29Z<p>Prozanski: /* Svarog */</p>
<hr />
<div>==Narzędzia wykorzystywane na ćwiczeniach==<br />
===Python===<br />
<br />
Python jest językiem programowania wysokiego poziomu, który w połączeniu z bibliotekami NumPy i SciPy do obliczeń naukowych pozwala na szybkie i wygodne programowanie lub analizowanie danych w sposób interaktywny. Przykłady prezentowane w ramach zajęć powinny uruchamiać się zarówno w wersji 2 jak i 3 języka Python, jednak zachęcamy Państwa, aby od początku uczyć się i korzystać z wersji 3 języka.<br />
<br />
==== Dokumentacja modułu scipy.signal ====<br />
Proszę zapoznać się z dokumentacją biblioteki scipy.signal:<br />
<br />
http://docs.scipy.org/doc/scipy-0.14.0/reference/signal.html<br />
<br />
===Svarog===<br />
Przydatnym narzędziem do analizy sygnałów, z którego będziemy korzystać na zajęciach, jest program SVAROG (pierwotnie skrót od Signal Viewer, Analyzer and Recorder On GPL). Program działa w środowisku Java, jest więc niezależny od systemu operacyjnego (Linux, Windows, OS X…). Svarog pozwala na wczytywanie i analizowanie sygnałów (nie tylko bioelektrycznych), zarówno przy użyciu prostych (FFT, spektrogram) jak i bardziej zaawansowanych (matching pursuit, ICA, DFT itd.) narzędzi. Dzięki współpracy z platformą OpenBCI, możliwa jest rejestracja sygnału (łącznie z metadanymi) bezpośrednio z poziomu graficznego interfejsu użytkownika.<br />
<br />
====Svarog: uruchamianie i konfiguracja====<br />
Aktualną wersję programu Svarog można pobrać [http://www.mimuw.edu.pl/~ptr/fid stąd]. Program nie wymaga instalacji. Po rozpakowaniu paczki do dowolnego katalogu należy uruchomić skrypt „run-svarog.sh” lub uruchomić bezpośrednio plik *.jar.<br />
<br />
W przypadku pracy na własnych komputerach, do prawidłowego uruchomienia pluginu do analizy sygnałów, z którego będziemy korzystać w dalszej części ćwiczeń, konieczne jest zainstalowanie środowiska Oracle Java SE w wersji 8, które można pobrać [http://www.oracle.com/technetwork/java/javase/downloads/index.html ze strony wydawcy]. Alternatywnie, użytkownicy systemu Ubuntu lub pokrewnych dystrybucji mogą zainstalować środowisko Java według instrukcji [http://www.webupd8.org/2012/09/install-oracle-java-8-in-ubuntu-via-ppa.html dostępnych na tej stronie].<br />
<br />
==Sygnały ciągłe i dyskretne ==<br />
===Próbkowanie w czasie ===<br />
Proszę powtórzyć sobie pojęcia:<br />
* częstość próbkowania<br />
* częstość Nyquista<br />
* aliasing<br />
<br />
W poniższym ćwiczeniu chcemy zbadać efekt próbkowania sygnału w czasie. W komputerach nie mamy dostępu do sygnału ciągłego. Na nasze potrzeby wygenerujemy sygnały próbkowane z bardzo dużą częstością, które będą dla nas aproksymacją sygnałów ciągłych. Przy ich pomocy zaprezentujemy efekt utożsamiania (aliasingu).<br />
<br />
Proszę wytworzyć wektor reprezentujący czas &bdquo;prawie&rdquo; ciągły. Będzie to u nas 1000 wartości z przedziału [0,1) wziętych z odstępem 0,001.<br />
<br />
<source lang =python><br />
# -*- coding: utf-8 -*-<br />
<br />
import pylab as py<br />
import numpy as np<br />
<br />
t = np.arange(0,1,0.001) # czas 'prawie ciągły'<br />
</source><br />
<br />
Teraz proszę wygenerować dwie sinusoidy: jedną o częstości <tt>-1</tt> a drugą o częstości <tt>9</tt>. Dla przypomnienia wyrażenie:<br />
:<math> s(t) = \sin(2 \pi f t)</math> możemy w pythonie zapisać:<br />
<source lang =python><br />
s = np.sin(2*np.pi*f*t)<br />
</source><br />
<br />
Proszę [[TI/Matplotlib#Kilka_wykres.C3.B3w_we_wsp.C3.B3lnych_osiach|wykreślić]] obie sinusoidy.<br />
<br />
Teraz proszę spróbkować czas i nasze &bdquo;prawie&rdquo; ciągłe sinusoidy z okresem próbkowania 0,1. (Trzeba pobrać co 100 element, proszę posłużyć się [[Programowanie_z_Pythonem/Sekwencje#Wycinki|wycinkami]]) <!-- [[TI/Sekwencje|Struktury danych — sekwencje|wycinkami]]. --><br />
Na tle &bdquo;prawie&rdquo; ciągłych sinusoid proszę dorysować punkty ze spróbkowanych sygnałów. Aby punkty były dobrze widoczne proponuję użyć markerów <tt>x</tt> oraz <tt>+</tt>.<br />
<br />
Proszę zaobserwować wzajemne położenie punktów. Czy można odróżnić sinusoidę o częstości &minus;1 od sinusoidy o częstości 9, jeśli obie są próbkowane z częstością 10? Jak można uogólnić tą obserwację?<br />
<tt>*</tt><br />
<!--<br />
# -*- coding: utf-8 -*-<br />
<br />
import pylab as py<br />
import numpy as np<br />
<br />
t = np.arange(0,1,0.001) # czas 'prawie ciągły'<br />
f1 = -1 # częstość sygnału 1<br />
f2 = 9 # częstość sygnału 2<br />
s1 = np.sin(2*np.pi*f1*t) # prawie ciągły sygnał o częstości f1<br />
s2 = np.sin(2*np.pi*f2*t) # prawie ciągły sygnał o częstości f2<br />
<br />
T = 0.1 # okres próbkowania<br />
Fs = 1/T # częstość próbkowania<br />
FN = Fs/2 # częstość Nyquista<br />
T_samp = t[0::100] # czas pobierania próbek<br />
s1_samp = s1[0::100] # próbkowany sygnał o częstości f1<br />
s2_samp = s2[0::100] # próbkowany sygnał o częstości f2<br />
<br />
py.plot(t, s1, 'g')<br />
py.plot(t, s2, 'b')<br />
py.plot(T_samp, s1_samp, 'gx', markersize=10)<br />
py.plot(T_samp, s2_samp, 'r+', markersize=10)<br />
--><br />
<br />
<!--<br />
===Błąd kwantyzacji ===<br />
Kiedy mierzymy fizyczne wielkości w celu dalszej analizy najczęściej chcemy przypisać im pewne liczby. Liczby w systemach cyfrowych reprezentowane są ze skończoną dokładnością. Urządzenia dokonujące przypisania liczby do mierzonej wartości to przetworniki analogowo-cyfrowe (ang. ADC, Analog to Digital Converter).<br />
Charakteryzują się one określoną ilością bitów ''N'', za pomocą których reprezentują liczby. Pełen zakres wartości pomiarowych ''R'' jest dzielony na <math>2^N</math> poziomów. Błąd kwantyzacji szacujemy jako nie większy niż <math> \frac{R}{2^{N+1}} </math>. <br />
<br />
====Zadanie:====<br />
Proszę oszacować błąd kwantyzacji sygnału EEG mierzonego za pomocą 12-bitowego konwertera. Zakres pomiarowy tego urządzenia to &plusmn;200 &mu;V.<br />
====Zadanie:====<br />
Proszę zilustrować efekt kwantyzacji dla trzybitowego przetwornika o zakresie 2.<br />
W tym celu proszę wykonać następujące kroki:<br />
# wygenerować &bdquo;prawie&rdquo; ciągłą sinusoidę (częstość 1, czas trwania 1 próbkowanie co 0,001)<br />
# spróbkować tą sinusoidę co 0,1 (proszę zastosować wycinki)<br />
# proszę skwantować spróbkowane wartości <br />
#Proszę wykreślić na jednym rysunku <br />
#* oryginalny sygnał <br />
#* sygnał spróbkowany w czasie<br />
#* sygnał spróbkowany w czasie i o skwantowanej amplitudzie (skorzystać z funkcji <tt>py.step</tt>)<br />
<br />
{{ Wyjaśnienie|title= wskazówka: Kwantowanie | text = <br />
<source lang =python><br />
N_bits = 3<br />
zakres = 2.0<br />
dy = zakres/2**N_bits<br />
s1_kwantowany = np.floor(s1_samp/dy)*dy +0.5*dy<br />
</source><br />
}}<br />
<br />
<!--<br />
# -*- coding: utf-8 -*-<br />
<br />
import pylab as py<br />
import numpy as np<br />
<br />
t = np.arange(0,1,0.001) # czas 'prawie ciągły'<br />
f1 = 1 # częstość sygnału 1<br />
<br />
s1 = np.sin(2*np.pi*f1*t) # prawie ciągły sygnał o częstości f1<br />
T_samp = t[0::100]<br />
s1_samp = s1[0::100]<br />
<br />
N_bits = 3<br />
zakres = 2.0<br />
dy = zakres/2**N_bits<br />
s1_kwantowany = np.floor(s1_samp/dy)*dy +0.5*dy<br />
ax = py.subplot(111)<br />
py.plot(t,s1)<br />
py.plot(T_samp,s1_samp,'ko',markersize = 5)<br />
py.step(T_samp,s1_kwantowany,where = 'post')<br />
<br />
ax.set_xticks(T_samp)<br />
poziomy = np.arange(-1+0.5*dy,1,dy)<br />
ax.set_yticks(poziomy)<br />
py.grid('on')<br />
<br />
<br />
for i in xrange(0,8):<br />
py.text(1.02,poziomy[i],bin(i))<br />
--><br />
<br />
==Sygnały testowe==<br />
<br />
<br />
=== Generowanie sygnałów testowych ===<br />
<br />
Do badania różnych metod analizy sygnałów potrzebne nam będą sygnały o znanych własnościach. W szczególności dobrze jest umnieć nadać sygnałom występującym w postaci cyfrowej, oraz sztucznym sygnałom próbnym pewne własności fizyczne takie jak:<br />
<br />
<ul><br />
<br />
<li><br />
czętość próbkowania<br />
<br />
<br />
<li><br />
czas trwania<br />
<br />
<br />
<li><br />
amplituda<br />
<br />
</ul><br />
=== Przykład sinus===<br />
Sinus o zadanej częstości (w Hz), długości trwania, częstości próbkowania i fazie.<br />
Poniższy kod implementuje i testuje funkcję <br />
:<math> \sin(f,T,Fs,\phi) = \sin(2*\pi f t)</math> dla <math>t \in \{0,T\}</math> <br />
<br />
<source lang= python><br />
# -*- coding: utf-8 -*-<br />
<br />
import pylab as py<br />
import numpy as np<br />
<br />
def sin(f = 1, T = 1, Fs = 128, phi =0 ):<br />
'''sin o zadanej częstości (w Hz), długości, fazie i częstości próbkowania<br />
Domyślnie wytwarzany jest sygnał reprezentujący <br />
1 sekundę sinusa o częstości 1 Hz i zerowej fazie próbkowanego 128 Hz<br />
'''<br />
<br />
dt = 1.0/Fs<br />
t = np.arange(0,T,dt)<br />
s = np.sin(2*np.pi*f*t + phi)<br />
return (s,t)<br />
<br />
<br />
(s,t) = sin(f=10,Fs=1000)<br />
py.plot(t,s)<br />
py.show()<br />
</source><br />
=== Przykład: eksport sygnału do pliku binarnego ===<br />
* Poniższy kod ilustruje sposób zapisu dwóch funkcji sinus o częstościach 10 Hz i 21 Hz do pliku binarnego:<br />
<br />
<source lang= python><br />
# -*- coding: utf-8 -*-<br />
<br />
import numpy as np<br />
<br />
T = 5<br />
Fs = 128.0<br />
<br />
s1 = sin(f=10, T=T, Fs=Fs)<br />
s2 = sin(f=21, T=T, Fs=Fs) <br />
<br />
signal = np.zeros((T*Fs, 2), dtype='<f')<br />
signal[:, 0] = s1<br />
signal[:, 1] = s2<br />
<br />
with open('test_signal.bin', 'wb') as f:<br />
signal.tofile(f)<br />
<br />
</source><br />
<br />
===Przykład: wczytanie sygnału do Svaroga ===<br />
W celu wczytania zapisanego binarnie sygnału do programu Svarog, po wybraniu File -> Open signal, należy wprowadzić częstość próbkowania sygnału oraz liczbę kanałów. <br />
<br />
[[Plik:svarog_open_signal.png|center|800px|thumb|<figure id="uid3" />Wczytywanie sygnału w programie Svarog.]]<br />
<br />
=== Delta ===<br />
Podobnie można zdefiniować funkcję delta o zadanym czasie trwania, częstości próbkowania i momencie wystąpienia impulsu:<br />
<br />
:<math> \delta(t_0) = \left\{^{1 \quad t=t_0} _{0 \quad t \ne t_0} \right.</math><br />
<br />
<source lang= python><br />
def delta(t0=0.5, T=1 ,Fs = 128):<br />
dt = 1.0/Fs<br />
t = np.arange(0,T,dt)<br />
d = np.zeros(len(t))<br />
d[np.ceil(t0*Fs)]=1<br />
return (d,t)<br />
</source><br />
<br />
=== Zadanie:===<br />
Analogicznie do powyższych przykładów proszę zaimplementować i przetestować funkcje generujące:<br />
* funkcję Gabora (funkcja Gaussa modulowana cosinusem) o zadanej częstości i standardowym odchyleniu w czasie, momencie wystąpienia, długości, częstości próbkowania i fazie. <br />
:<math> g = \exp\left(-\frac{1}{2}\left(\frac{t-t_0}{\sigma}\right)^2 \right) \cdot \cos(2 \pi f t + \phi); </math><br />
<br />
* szum gaussowski o zadanej średniej, odchyleniu standardowym, długości i częstości próbkowania.<br />
<br />
<!-- * pochodną funkcji Gaussa<br />
<br />
* połówkę funkcji Gaussa<br />
--></div>Prozanskihttp://brain.fuw.edu.pl/edu/index.php?title=%C4%86wiczenia_1&diff=4742Ćwiczenia 12016-02-18T01:47:56Z<p>Prozanski: /* Python */</p>
<hr />
<div>==Narzędzia wykorzystywane na ćwiczeniach==<br />
===Python===<br />
<br />
Python jest językiem programowania wysokiego poziomu, który w połączeniu z bibliotekami NumPy i SciPy do obliczeń naukowych pozwala na szybkie i wygodne programowanie lub analizowanie danych w sposób interaktywny. Przykłady prezentowane w ramach zajęć powinny uruchamiać się zarówno w wersji 2 jak i 3 języka Python, jednak zachęcamy Państwa, aby od początku uczyć się i korzystać z wersji 3 języka.<br />
<br />
==== Dokumentacja modułu scipy.signal ====<br />
Proszę zapoznać się z dokumentacją biblioteki scipy.signal:<br />
<br />
http://docs.scipy.org/doc/scipy-0.14.0/reference/signal.html<br />
<br />
===Svarog===<br />
Przydatnym narzędziem do analizy sygnałów, z którego będziemy korzystać na zajęciach, jest program SVAROG (pierwotnie skrót od Signal Viewer, Analyzer and Recorder On GPL). Program działa w środowisku Java, jest więc niezależny od systemu operacyjnego (Linux, Windows, OS X…). Svarog pozwala na wczytywanie i analizowanie sygnałów (nie tylko bioelektrycznych), zarówno przy użyciu prostych (FFT, spektrogram) jak i bardziej zaawansowanych (matching pursuit, ICA, DFT itd.) narzędzi. Dzięki współpracy z platformą OpenBCI, możliwa jest rejestracja sygnału (łącznie z metadanymi) bezpośrednio z poziomu graficznego interfejsu użytkownika.<br />
<br />
Signal Viewer, ...<br />
====Svarog: uruchamianie i konfiguracja====<br />
Aktualną wersję programu Svarog można pobrać [http://www.mimuw.edu.pl/~ptr/fid stąd]. Program nie wymaga instalacji. Po rozpakowaniu paczki do dowolnego katalogu należy uruchomić skrypt „run-svarog.sh” lub uruchomić bezpośrednio plik *.jar.<br />
<br />
W przypadku pracy na własnych komputerach, do prawidłowego uruchomienia pluginu do analizy sygnałów, z którego będziemy korzystać w dalszej części ćwiczeń, konieczne jest zainstalowanie środowiska Oracle Java SE w wersji 8, które można pobrać [http://www.oracle.com/technetwork/java/javase/downloads/index.html ze strony wydawcy]. Alternatywnie, użytkownicy systemu Ubuntu lub pokrewnych dystrybucji mogą zainstalować środowisko Java według instrukcji [http://www.webupd8.org/2012/09/install-oracle-java-8-in-ubuntu-via-ppa.html dostępnych na tej stronie].<br />
<br />
==Sygnały ciągłe i dyskretne ==<br />
===Próbkowanie w czasie ===<br />
Proszę powtórzyć sobie pojęcia:<br />
* częstość próbkowania<br />
* częstość Nyquista<br />
* aliasing<br />
<br />
W poniższym ćwiczeniu chcemy zbadać efekt próbkowania sygnału w czasie. W komputerach nie mamy dostępu do sygnału ciągłego. Na nasze potrzeby wygenerujemy sygnały próbkowane z bardzo dużą częstością, które będą dla nas aproksymacją sygnałów ciągłych. Przy ich pomocy zaprezentujemy efekt utożsamiania (aliasingu).<br />
<br />
Proszę wytworzyć wektor reprezentujący czas &bdquo;prawie&rdquo; ciągły. Będzie to u nas 1000 wartości z przedziału [0,1) wziętych z odstępem 0,001.<br />
<br />
<source lang =python><br />
# -*- coding: utf-8 -*-<br />
<br />
import pylab as py<br />
import numpy as np<br />
<br />
t = np.arange(0,1,0.001) # czas 'prawie ciągły'<br />
</source><br />
<br />
Teraz proszę wygenerować dwie sinusoidy: jedną o częstości <tt>-1</tt> a drugą o częstości <tt>9</tt>. Dla przypomnienia wyrażenie:<br />
:<math> s(t) = \sin(2 \pi f t)</math> możemy w pythonie zapisać:<br />
<source lang =python><br />
s = np.sin(2*np.pi*f*t)<br />
</source><br />
<br />
Proszę [[TI/Matplotlib#Kilka_wykres.C3.B3w_we_wsp.C3.B3lnych_osiach|wykreślić]] obie sinusoidy.<br />
<br />
Teraz proszę spróbkować czas i nasze &bdquo;prawie&rdquo; ciągłe sinusoidy z okresem próbkowania 0,1. (Trzeba pobrać co 100 element, proszę posłużyć się [[Programowanie_z_Pythonem/Sekwencje#Wycinki|wycinkami]]) <!-- [[TI/Sekwencje|Struktury danych — sekwencje|wycinkami]]. --><br />
Na tle &bdquo;prawie&rdquo; ciągłych sinusoid proszę dorysować punkty ze spróbkowanych sygnałów. Aby punkty były dobrze widoczne proponuję użyć markerów <tt>x</tt> oraz <tt>+</tt>.<br />
<br />
Proszę zaobserwować wzajemne położenie punktów. Czy można odróżnić sinusoidę o częstości &minus;1 od sinusoidy o częstości 9, jeśli obie są próbkowane z częstością 10? Jak można uogólnić tą obserwację?<br />
<tt>*</tt><br />
<!--<br />
# -*- coding: utf-8 -*-<br />
<br />
import pylab as py<br />
import numpy as np<br />
<br />
t = np.arange(0,1,0.001) # czas 'prawie ciągły'<br />
f1 = -1 # częstość sygnału 1<br />
f2 = 9 # częstość sygnału 2<br />
s1 = np.sin(2*np.pi*f1*t) # prawie ciągły sygnał o częstości f1<br />
s2 = np.sin(2*np.pi*f2*t) # prawie ciągły sygnał o częstości f2<br />
<br />
T = 0.1 # okres próbkowania<br />
Fs = 1/T # częstość próbkowania<br />
FN = Fs/2 # częstość Nyquista<br />
T_samp = t[0::100] # czas pobierania próbek<br />
s1_samp = s1[0::100] # próbkowany sygnał o częstości f1<br />
s2_samp = s2[0::100] # próbkowany sygnał o częstości f2<br />
<br />
py.plot(t, s1, 'g')<br />
py.plot(t, s2, 'b')<br />
py.plot(T_samp, s1_samp, 'gx', markersize=10)<br />
py.plot(T_samp, s2_samp, 'r+', markersize=10)<br />
--><br />
<br />
<!--<br />
===Błąd kwantyzacji ===<br />
Kiedy mierzymy fizyczne wielkości w celu dalszej analizy najczęściej chcemy przypisać im pewne liczby. Liczby w systemach cyfrowych reprezentowane są ze skończoną dokładnością. Urządzenia dokonujące przypisania liczby do mierzonej wartości to przetworniki analogowo-cyfrowe (ang. ADC, Analog to Digital Converter).<br />
Charakteryzują się one określoną ilością bitów ''N'', za pomocą których reprezentują liczby. Pełen zakres wartości pomiarowych ''R'' jest dzielony na <math>2^N</math> poziomów. Błąd kwantyzacji szacujemy jako nie większy niż <math> \frac{R}{2^{N+1}} </math>. <br />
<br />
====Zadanie:====<br />
Proszę oszacować błąd kwantyzacji sygnału EEG mierzonego za pomocą 12-bitowego konwertera. Zakres pomiarowy tego urządzenia to &plusmn;200 &mu;V.<br />
====Zadanie:====<br />
Proszę zilustrować efekt kwantyzacji dla trzybitowego przetwornika o zakresie 2.<br />
W tym celu proszę wykonać następujące kroki:<br />
# wygenerować &bdquo;prawie&rdquo; ciągłą sinusoidę (częstość 1, czas trwania 1 próbkowanie co 0,001)<br />
# spróbkować tą sinusoidę co 0,1 (proszę zastosować wycinki)<br />
# proszę skwantować spróbkowane wartości <br />
#Proszę wykreślić na jednym rysunku <br />
#* oryginalny sygnał <br />
#* sygnał spróbkowany w czasie<br />
#* sygnał spróbkowany w czasie i o skwantowanej amplitudzie (skorzystać z funkcji <tt>py.step</tt>)<br />
<br />
{{ Wyjaśnienie|title= wskazówka: Kwantowanie | text = <br />
<source lang =python><br />
N_bits = 3<br />
zakres = 2.0<br />
dy = zakres/2**N_bits<br />
s1_kwantowany = np.floor(s1_samp/dy)*dy +0.5*dy<br />
</source><br />
}}<br />
<br />
<!--<br />
# -*- coding: utf-8 -*-<br />
<br />
import pylab as py<br />
import numpy as np<br />
<br />
t = np.arange(0,1,0.001) # czas 'prawie ciągły'<br />
f1 = 1 # częstość sygnału 1<br />
<br />
s1 = np.sin(2*np.pi*f1*t) # prawie ciągły sygnał o częstości f1<br />
T_samp = t[0::100]<br />
s1_samp = s1[0::100]<br />
<br />
N_bits = 3<br />
zakres = 2.0<br />
dy = zakres/2**N_bits<br />
s1_kwantowany = np.floor(s1_samp/dy)*dy +0.5*dy<br />
ax = py.subplot(111)<br />
py.plot(t,s1)<br />
py.plot(T_samp,s1_samp,'ko',markersize = 5)<br />
py.step(T_samp,s1_kwantowany,where = 'post')<br />
<br />
ax.set_xticks(T_samp)<br />
poziomy = np.arange(-1+0.5*dy,1,dy)<br />
ax.set_yticks(poziomy)<br />
py.grid('on')<br />
<br />
<br />
for i in xrange(0,8):<br />
py.text(1.02,poziomy[i],bin(i))<br />
--><br />
<br />
==Sygnały testowe==<br />
<br />
<br />
=== Generowanie sygnałów testowych ===<br />
<br />
Do badania różnych metod analizy sygnałów potrzebne nam będą sygnały o znanych własnościach. W szczególności dobrze jest umnieć nadać sygnałom występującym w postaci cyfrowej, oraz sztucznym sygnałom próbnym pewne własności fizyczne takie jak:<br />
<br />
<ul><br />
<br />
<li><br />
czętość próbkowania<br />
<br />
<br />
<li><br />
czas trwania<br />
<br />
<br />
<li><br />
amplituda<br />
<br />
</ul><br />
=== Przykład sinus===<br />
Sinus o zadanej częstości (w Hz), długości trwania, częstości próbkowania i fazie.<br />
Poniższy kod implementuje i testuje funkcję <br />
:<math> \sin(f,T,Fs,\phi) = \sin(2*\pi f t)</math> dla <math>t \in \{0,T\}</math> <br />
<br />
<source lang= python><br />
# -*- coding: utf-8 -*-<br />
<br />
import pylab as py<br />
import numpy as np<br />
<br />
def sin(f = 1, T = 1, Fs = 128, phi =0 ):<br />
'''sin o zadanej częstości (w Hz), długości, fazie i częstości próbkowania<br />
Domyślnie wytwarzany jest sygnał reprezentujący <br />
1 sekundę sinusa o częstości 1 Hz i zerowej fazie próbkowanego 128 Hz<br />
'''<br />
<br />
dt = 1.0/Fs<br />
t = np.arange(0,T,dt)<br />
s = np.sin(2*np.pi*f*t + phi)<br />
return (s,t)<br />
<br />
<br />
(s,t) = sin(f=10,Fs=1000)<br />
py.plot(t,s)<br />
py.show()<br />
</source><br />
=== Przykład: eksport sygnału do pliku binarnego ===<br />
* Poniższy kod ilustruje sposób zapisu dwóch funkcji sinus o częstościach 10 Hz i 21 Hz do pliku binarnego:<br />
<br />
<source lang= python><br />
# -*- coding: utf-8 -*-<br />
<br />
import numpy as np<br />
<br />
T = 5<br />
Fs = 128.0<br />
<br />
s1 = sin(f=10, T=T, Fs=Fs)<br />
s2 = sin(f=21, T=T, Fs=Fs) <br />
<br />
signal = np.zeros((T*Fs, 2), dtype='<f')<br />
signal[:, 0] = s1<br />
signal[:, 1] = s2<br />
<br />
with open('test_signal.bin', 'wb') as f:<br />
signal.tofile(f)<br />
<br />
</source><br />
<br />
===Przykład: wczytanie sygnału do Svaroga ===<br />
W celu wczytania zapisanego binarnie sygnału do programu Svarog, po wybraniu File -> Open signal, należy wprowadzić częstość próbkowania sygnału oraz liczbę kanałów. <br />
<br />
[[Plik:svarog_open_signal.png|center|800px|thumb|<figure id="uid3" />Wczytywanie sygnału w programie Svarog.]]<br />
<br />
=== Delta ===<br />
Podobnie można zdefiniować funkcję delta o zadanym czasie trwania, częstości próbkowania i momencie wystąpienia impulsu:<br />
<br />
:<math> \delta(t_0) = \left\{^{1 \quad t=t_0} _{0 \quad t \ne t_0} \right.</math><br />
<br />
<source lang= python><br />
def delta(t0=0.5, T=1 ,Fs = 128):<br />
dt = 1.0/Fs<br />
t = np.arange(0,T,dt)<br />
d = np.zeros(len(t))<br />
d[np.ceil(t0*Fs)]=1<br />
return (d,t)<br />
</source><br />
<br />
=== Zadanie:===<br />
Analogicznie do powyższych przykładów proszę zaimplementować i przetestować funkcje generujące:<br />
* funkcję Gabora (funkcja Gaussa modulowana cosinusem) o zadanej częstości i standardowym odchyleniu w czasie, momencie wystąpienia, długości, częstości próbkowania i fazie. <br />
:<math> g = \exp\left(-\frac{1}{2}\left(\frac{t-t_0}{\sigma}\right)^2 \right) \cdot \cos(2 \pi f t + \phi); </math><br />
<br />
* szum gaussowski o zadanej średniej, odchyleniu standardowym, długości i częstości próbkowania.<br />
<br />
<!-- * pochodną funkcji Gaussa<br />
<br />
* połówkę funkcji Gaussa<br />
--></div>Prozanskihttp://brain.fuw.edu.pl/edu/index.php?title=%C4%86wiczenia_1&diff=4741Ćwiczenia 12016-02-18T01:40:41Z<p>Prozanski: /* Python */</p>
<hr />
<div>==Narzędzia wykorzystywane na ćwiczeniach==<br />
===Python===<br />
<br />
Python jest językiem wysokiego poziomu, który w połączeniu z bibliotekami NumPy i SciPy do obliczeń naukowych pozwala na szybkie i wygodne programowanie lub analizowanie danych w sposób interaktywny. Przykłady prezentowane w ramach zajęć powinny uruchamiać się zarówno w wersji 2 jak i 3 języka Python, jednak zachęcamy Państwa, aby od początku uczyć się i korzystać z wersji 3 języka.<br />
<br />
==== Dokumentacja modułu scipy.signal ====<br />
Proszę zapoznać się z dokumentacją biblioteki scipy.signal:<br />
<br />
http://docs.scipy.org/doc/scipy-0.14.0/reference/signal.html<br />
<br />
===Svarog===<br />
Przydatnym narzędziem do analizy sygnałów, z którego będziemy korzystać na zajęciach, jest program SVAROG (pierwotnie skrót od Signal Viewer, Analyzer and Recorder On GPL). Program działa w środowisku Java, jest więc niezależny od systemu operacyjnego (Linux, Windows, OS X…). Svarog pozwala na wczytywanie i analizowanie sygnałów (nie tylko bioelektrycznych), zarówno przy użyciu prostych (FFT, spektrogram) jak i bardziej zaawansowanych (matching pursuit, ICA, DFT itd.) narzędzi. Dzięki współpracy z platformą OpenBCI, możliwa jest rejestracja sygnału (łącznie z metadanymi) bezpośrednio z poziomu graficznego interfejsu użytkownika.<br />
<br />
Signal Viewer, ...<br />
====Svarog: uruchamianie i konfiguracja====<br />
Aktualną wersję programu Svarog można pobrać [http://www.mimuw.edu.pl/~ptr/fid stąd]. Program nie wymaga instalacji. Po rozpakowaniu paczki do dowolnego katalogu należy uruchomić skrypt „run-svarog.sh” lub uruchomić bezpośrednio plik *.jar.<br />
<br />
W przypadku pracy na własnych komputerach, do prawidłowego uruchomienia pluginu do analizy sygnałów, z którego będziemy korzystać w dalszej części ćwiczeń, konieczne jest zainstalowanie środowiska Oracle Java SE w wersji 8, które można pobrać [http://www.oracle.com/technetwork/java/javase/downloads/index.html ze strony wydawcy]. Alternatywnie, użytkownicy systemu Ubuntu lub pokrewnych dystrybucji mogą zainstalować środowisko Java według instrukcji [http://www.webupd8.org/2012/09/install-oracle-java-8-in-ubuntu-via-ppa.html dostępnych na tej stronie].<br />
<br />
==Sygnały ciągłe i dyskretne ==<br />
===Próbkowanie w czasie ===<br />
Proszę powtórzyć sobie pojęcia:<br />
* częstość próbkowania<br />
* częstość Nyquista<br />
* aliasing<br />
<br />
W poniższym ćwiczeniu chcemy zbadać efekt próbkowania sygnału w czasie. W komputerach nie mamy dostępu do sygnału ciągłego. Na nasze potrzeby wygenerujemy sygnały próbkowane z bardzo dużą częstością, które będą dla nas aproksymacją sygnałów ciągłych. Przy ich pomocy zaprezentujemy efekt utożsamiania (aliasingu).<br />
<br />
Proszę wytworzyć wektor reprezentujący czas &bdquo;prawie&rdquo; ciągły. Będzie to u nas 1000 wartości z przedziału [0,1) wziętych z odstępem 0,001.<br />
<br />
<source lang =python><br />
# -*- coding: utf-8 -*-<br />
<br />
import pylab as py<br />
import numpy as np<br />
<br />
t = np.arange(0,1,0.001) # czas 'prawie ciągły'<br />
</source><br />
<br />
Teraz proszę wygenerować dwie sinusoidy: jedną o częstości <tt>-1</tt> a drugą o częstości <tt>9</tt>. Dla przypomnienia wyrażenie:<br />
:<math> s(t) = \sin(2 \pi f t)</math> możemy w pythonie zapisać:<br />
<source lang =python><br />
s = np.sin(2*np.pi*f*t)<br />
</source><br />
<br />
Proszę [[TI/Matplotlib#Kilka_wykres.C3.B3w_we_wsp.C3.B3lnych_osiach|wykreślić]] obie sinusoidy.<br />
<br />
Teraz proszę spróbkować czas i nasze &bdquo;prawie&rdquo; ciągłe sinusoidy z okresem próbkowania 0,1. (Trzeba pobrać co 100 element, proszę posłużyć się [[Programowanie_z_Pythonem/Sekwencje#Wycinki|wycinkami]]) <!-- [[TI/Sekwencje|Struktury danych — sekwencje|wycinkami]]. --><br />
Na tle &bdquo;prawie&rdquo; ciągłych sinusoid proszę dorysować punkty ze spróbkowanych sygnałów. Aby punkty były dobrze widoczne proponuję użyć markerów <tt>x</tt> oraz <tt>+</tt>.<br />
<br />
Proszę zaobserwować wzajemne położenie punktów. Czy można odróżnić sinusoidę o częstości &minus;1 od sinusoidy o częstości 9, jeśli obie są próbkowane z częstością 10? Jak można uogólnić tą obserwację?<br />
<tt>*</tt><br />
<!--<br />
# -*- coding: utf-8 -*-<br />
<br />
import pylab as py<br />
import numpy as np<br />
<br />
t = np.arange(0,1,0.001) # czas 'prawie ciągły'<br />
f1 = -1 # częstość sygnału 1<br />
f2 = 9 # częstość sygnału 2<br />
s1 = np.sin(2*np.pi*f1*t) # prawie ciągły sygnał o częstości f1<br />
s2 = np.sin(2*np.pi*f2*t) # prawie ciągły sygnał o częstości f2<br />
<br />
T = 0.1 # okres próbkowania<br />
Fs = 1/T # częstość próbkowania<br />
FN = Fs/2 # częstość Nyquista<br />
T_samp = t[0::100] # czas pobierania próbek<br />
s1_samp = s1[0::100] # próbkowany sygnał o częstości f1<br />
s2_samp = s2[0::100] # próbkowany sygnał o częstości f2<br />
<br />
py.plot(t, s1, 'g')<br />
py.plot(t, s2, 'b')<br />
py.plot(T_samp, s1_samp, 'gx', markersize=10)<br />
py.plot(T_samp, s2_samp, 'r+', markersize=10)<br />
--><br />
<br />
<!--<br />
===Błąd kwantyzacji ===<br />
Kiedy mierzymy fizyczne wielkości w celu dalszej analizy najczęściej chcemy przypisać im pewne liczby. Liczby w systemach cyfrowych reprezentowane są ze skończoną dokładnością. Urządzenia dokonujące przypisania liczby do mierzonej wartości to przetworniki analogowo-cyfrowe (ang. ADC, Analog to Digital Converter).<br />
Charakteryzują się one określoną ilością bitów ''N'', za pomocą których reprezentują liczby. Pełen zakres wartości pomiarowych ''R'' jest dzielony na <math>2^N</math> poziomów. Błąd kwantyzacji szacujemy jako nie większy niż <math> \frac{R}{2^{N+1}} </math>. <br />
<br />
====Zadanie:====<br />
Proszę oszacować błąd kwantyzacji sygnału EEG mierzonego za pomocą 12-bitowego konwertera. Zakres pomiarowy tego urządzenia to &plusmn;200 &mu;V.<br />
====Zadanie:====<br />
Proszę zilustrować efekt kwantyzacji dla trzybitowego przetwornika o zakresie 2.<br />
W tym celu proszę wykonać następujące kroki:<br />
# wygenerować &bdquo;prawie&rdquo; ciągłą sinusoidę (częstość 1, czas trwania 1 próbkowanie co 0,001)<br />
# spróbkować tą sinusoidę co 0,1 (proszę zastosować wycinki)<br />
# proszę skwantować spróbkowane wartości <br />
#Proszę wykreślić na jednym rysunku <br />
#* oryginalny sygnał <br />
#* sygnał spróbkowany w czasie<br />
#* sygnał spróbkowany w czasie i o skwantowanej amplitudzie (skorzystać z funkcji <tt>py.step</tt>)<br />
<br />
{{ Wyjaśnienie|title= wskazówka: Kwantowanie | text = <br />
<source lang =python><br />
N_bits = 3<br />
zakres = 2.0<br />
dy = zakres/2**N_bits<br />
s1_kwantowany = np.floor(s1_samp/dy)*dy +0.5*dy<br />
</source><br />
}}<br />
<br />
<!--<br />
# -*- coding: utf-8 -*-<br />
<br />
import pylab as py<br />
import numpy as np<br />
<br />
t = np.arange(0,1,0.001) # czas 'prawie ciągły'<br />
f1 = 1 # częstość sygnału 1<br />
<br />
s1 = np.sin(2*np.pi*f1*t) # prawie ciągły sygnał o częstości f1<br />
T_samp = t[0::100]<br />
s1_samp = s1[0::100]<br />
<br />
N_bits = 3<br />
zakres = 2.0<br />
dy = zakres/2**N_bits<br />
s1_kwantowany = np.floor(s1_samp/dy)*dy +0.5*dy<br />
ax = py.subplot(111)<br />
py.plot(t,s1)<br />
py.plot(T_samp,s1_samp,'ko',markersize = 5)<br />
py.step(T_samp,s1_kwantowany,where = 'post')<br />
<br />
ax.set_xticks(T_samp)<br />
poziomy = np.arange(-1+0.5*dy,1,dy)<br />
ax.set_yticks(poziomy)<br />
py.grid('on')<br />
<br />
<br />
for i in xrange(0,8):<br />
py.text(1.02,poziomy[i],bin(i))<br />
--><br />
<br />
==Sygnały testowe==<br />
<br />
<br />
=== Generowanie sygnałów testowych ===<br />
<br />
Do badania różnych metod analizy sygnałów potrzebne nam będą sygnały o znanych własnościach. W szczególności dobrze jest umnieć nadać sygnałom występującym w postaci cyfrowej, oraz sztucznym sygnałom próbnym pewne własności fizyczne takie jak:<br />
<br />
<ul><br />
<br />
<li><br />
czętość próbkowania<br />
<br />
<br />
<li><br />
czas trwania<br />
<br />
<br />
<li><br />
amplituda<br />
<br />
</ul><br />
=== Przykład sinus===<br />
Sinus o zadanej częstości (w Hz), długości trwania, częstości próbkowania i fazie.<br />
Poniższy kod implementuje i testuje funkcję <br />
:<math> \sin(f,T,Fs,\phi) = \sin(2*\pi f t)</math> dla <math>t \in \{0,T\}</math> <br />
<br />
<source lang= python><br />
# -*- coding: utf-8 -*-<br />
<br />
import pylab as py<br />
import numpy as np<br />
<br />
def sin(f = 1, T = 1, Fs = 128, phi =0 ):<br />
'''sin o zadanej częstości (w Hz), długości, fazie i częstości próbkowania<br />
Domyślnie wytwarzany jest sygnał reprezentujący <br />
1 sekundę sinusa o częstości 1 Hz i zerowej fazie próbkowanego 128 Hz<br />
'''<br />
<br />
dt = 1.0/Fs<br />
t = np.arange(0,T,dt)<br />
s = np.sin(2*np.pi*f*t + phi)<br />
return (s,t)<br />
<br />
<br />
(s,t) = sin(f=10,Fs=1000)<br />
py.plot(t,s)<br />
py.show()<br />
</source><br />
=== Przykład: eksport sygnału do pliku binarnego ===<br />
* Poniższy kod ilustruje sposób zapisu dwóch funkcji sinus o częstościach 10 Hz i 21 Hz do pliku binarnego:<br />
<br />
<source lang= python><br />
# -*- coding: utf-8 -*-<br />
<br />
import numpy as np<br />
<br />
T = 5<br />
Fs = 128.0<br />
<br />
s1 = sin(f=10, T=T, Fs=Fs)<br />
s2 = sin(f=21, T=T, Fs=Fs) <br />
<br />
signal = np.zeros((T*Fs, 2), dtype='<f')<br />
signal[:, 0] = s1<br />
signal[:, 1] = s2<br />
<br />
with open('test_signal.bin', 'wb') as f:<br />
signal.tofile(f)<br />
<br />
</source><br />
<br />
===Przykład: wczytanie sygnału do Svaroga ===<br />
W celu wczytania zapisanego binarnie sygnału do programu Svarog, po wybraniu File -> Open signal, należy wprowadzić częstość próbkowania sygnału oraz liczbę kanałów. <br />
<br />
[[Plik:svarog_open_signal.png|center|800px|thumb|<figure id="uid3" />Wczytywanie sygnału w programie Svarog.]]<br />
<br />
=== Delta ===<br />
Podobnie można zdefiniować funkcję delta o zadanym czasie trwania, częstości próbkowania i momencie wystąpienia impulsu:<br />
<br />
:<math> \delta(t_0) = \left\{^{1 \quad t=t_0} _{0 \quad t \ne t_0} \right.</math><br />
<br />
<source lang= python><br />
def delta(t0=0.5, T=1 ,Fs = 128):<br />
dt = 1.0/Fs<br />
t = np.arange(0,T,dt)<br />
d = np.zeros(len(t))<br />
d[np.ceil(t0*Fs)]=1<br />
return (d,t)<br />
</source><br />
<br />
=== Zadanie:===<br />
Analogicznie do powyższych przykładów proszę zaimplementować i przetestować funkcje generujące:<br />
* funkcję Gabora (funkcja Gaussa modulowana cosinusem) o zadanej częstości i standardowym odchyleniu w czasie, momencie wystąpienia, długości, częstości próbkowania i fazie. <br />
:<math> g = \exp\left(-\frac{1}{2}\left(\frac{t-t_0}{\sigma}\right)^2 \right) \cdot \cos(2 \pi f t + \phi); </math><br />
<br />
* szum gaussowski o zadanej średniej, odchyleniu standardowym, długości i częstości próbkowania.<br />
<br />
<!-- * pochodną funkcji Gaussa<br />
<br />
* połówkę funkcji Gaussa<br />
--></div>Prozanskihttp://brain.fuw.edu.pl/edu/index.php?title=%C4%86wiczenia_1&diff=4740Ćwiczenia 12016-02-18T01:36:47Z<p>Prozanski: /* Svarog */</p>
<hr />
<div>==Narzędzia wykorzystywane na ćwiczeniach==<br />
===Python===<br />
==== Dokumentacja modułu scipy.signal ====<br />
Proszę zapoznać się z dokumentacją biblioteki scipy.signal:<br />
<br />
http://docs.scipy.org/doc/scipy-0.14.0/reference/signal.html<br />
<br />
===Svarog===<br />
Przydatnym narzędziem do analizy sygnałów, z którego będziemy korzystać na zajęciach, jest program SVAROG (pierwotnie skrót od Signal Viewer, Analyzer and Recorder On GPL). Program działa w środowisku Java, jest więc niezależny od systemu operacyjnego (Linux, Windows, OS X…). Svarog pozwala na wczytywanie i analizowanie sygnałów (nie tylko bioelektrycznych), zarówno przy użyciu prostych (FFT, spektrogram) jak i bardziej zaawansowanych (matching pursuit, ICA, DFT itd.) narzędzi. Dzięki współpracy z platformą OpenBCI, możliwa jest rejestracja sygnału (łącznie z metadanymi) bezpośrednio z poziomu graficznego interfejsu użytkownika.<br />
<br />
Signal Viewer, ...<br />
====Svarog: uruchamianie i konfiguracja====<br />
Aktualną wersję programu Svarog można pobrać [http://www.mimuw.edu.pl/~ptr/fid stąd]. Program nie wymaga instalacji. Po rozpakowaniu paczki do dowolnego katalogu należy uruchomić skrypt „run-svarog.sh” lub uruchomić bezpośrednio plik *.jar.<br />
<br />
W przypadku pracy na własnych komputerach, do prawidłowego uruchomienia pluginu do analizy sygnałów, z którego będziemy korzystać w dalszej części ćwiczeń, konieczne jest zainstalowanie środowiska Oracle Java SE w wersji 8, które można pobrać [http://www.oracle.com/technetwork/java/javase/downloads/index.html ze strony wydawcy]. Alternatywnie, użytkownicy systemu Ubuntu lub pokrewnych dystrybucji mogą zainstalować środowisko Java według instrukcji [http://www.webupd8.org/2012/09/install-oracle-java-8-in-ubuntu-via-ppa.html dostępnych na tej stronie].<br />
<br />
==Sygnały ciągłe i dyskretne ==<br />
===Próbkowanie w czasie ===<br />
Proszę powtórzyć sobie pojęcia:<br />
* częstość próbkowania<br />
* częstość Nyquista<br />
* aliasing<br />
<br />
W poniższym ćwiczeniu chcemy zbadać efekt próbkowania sygnału w czasie. W komputerach nie mamy dostępu do sygnału ciągłego. Na nasze potrzeby wygenerujemy sygnały próbkowane z bardzo dużą częstością, które będą dla nas aproksymacją sygnałów ciągłych. Przy ich pomocy zaprezentujemy efekt utożsamiania (aliasingu).<br />
<br />
Proszę wytworzyć wektor reprezentujący czas &bdquo;prawie&rdquo; ciągły. Będzie to u nas 1000 wartości z przedziału [0,1) wziętych z odstępem 0,001.<br />
<br />
<source lang =python><br />
# -*- coding: utf-8 -*-<br />
<br />
import pylab as py<br />
import numpy as np<br />
<br />
t = np.arange(0,1,0.001) # czas 'prawie ciągły'<br />
</source><br />
<br />
Teraz proszę wygenerować dwie sinusoidy: jedną o częstości <tt>-1</tt> a drugą o częstości <tt>9</tt>. Dla przypomnienia wyrażenie:<br />
:<math> s(t) = \sin(2 \pi f t)</math> możemy w pythonie zapisać:<br />
<source lang =python><br />
s = np.sin(2*np.pi*f*t)<br />
</source><br />
<br />
Proszę [[TI/Matplotlib#Kilka_wykres.C3.B3w_we_wsp.C3.B3lnych_osiach|wykreślić]] obie sinusoidy.<br />
<br />
Teraz proszę spróbkować czas i nasze &bdquo;prawie&rdquo; ciągłe sinusoidy z okresem próbkowania 0,1. (Trzeba pobrać co 100 element, proszę posłużyć się [[Programowanie_z_Pythonem/Sekwencje#Wycinki|wycinkami]]) <!-- [[TI/Sekwencje|Struktury danych — sekwencje|wycinkami]]. --><br />
Na tle &bdquo;prawie&rdquo; ciągłych sinusoid proszę dorysować punkty ze spróbkowanych sygnałów. Aby punkty były dobrze widoczne proponuję użyć markerów <tt>x</tt> oraz <tt>+</tt>.<br />
<br />
Proszę zaobserwować wzajemne położenie punktów. Czy można odróżnić sinusoidę o częstości &minus;1 od sinusoidy o częstości 9, jeśli obie są próbkowane z częstością 10? Jak można uogólnić tą obserwację?<br />
<tt>*</tt><br />
<!--<br />
# -*- coding: utf-8 -*-<br />
<br />
import pylab as py<br />
import numpy as np<br />
<br />
t = np.arange(0,1,0.001) # czas 'prawie ciągły'<br />
f1 = -1 # częstość sygnału 1<br />
f2 = 9 # częstość sygnału 2<br />
s1 = np.sin(2*np.pi*f1*t) # prawie ciągły sygnał o częstości f1<br />
s2 = np.sin(2*np.pi*f2*t) # prawie ciągły sygnał o częstości f2<br />
<br />
T = 0.1 # okres próbkowania<br />
Fs = 1/T # częstość próbkowania<br />
FN = Fs/2 # częstość Nyquista<br />
T_samp = t[0::100] # czas pobierania próbek<br />
s1_samp = s1[0::100] # próbkowany sygnał o częstości f1<br />
s2_samp = s2[0::100] # próbkowany sygnał o częstości f2<br />
<br />
py.plot(t, s1, 'g')<br />
py.plot(t, s2, 'b')<br />
py.plot(T_samp, s1_samp, 'gx', markersize=10)<br />
py.plot(T_samp, s2_samp, 'r+', markersize=10)<br />
--><br />
<br />
<!--<br />
===Błąd kwantyzacji ===<br />
Kiedy mierzymy fizyczne wielkości w celu dalszej analizy najczęściej chcemy przypisać im pewne liczby. Liczby w systemach cyfrowych reprezentowane są ze skończoną dokładnością. Urządzenia dokonujące przypisania liczby do mierzonej wartości to przetworniki analogowo-cyfrowe (ang. ADC, Analog to Digital Converter).<br />
Charakteryzują się one określoną ilością bitów ''N'', za pomocą których reprezentują liczby. Pełen zakres wartości pomiarowych ''R'' jest dzielony na <math>2^N</math> poziomów. Błąd kwantyzacji szacujemy jako nie większy niż <math> \frac{R}{2^{N+1}} </math>. <br />
<br />
====Zadanie:====<br />
Proszę oszacować błąd kwantyzacji sygnału EEG mierzonego za pomocą 12-bitowego konwertera. Zakres pomiarowy tego urządzenia to &plusmn;200 &mu;V.<br />
====Zadanie:====<br />
Proszę zilustrować efekt kwantyzacji dla trzybitowego przetwornika o zakresie 2.<br />
W tym celu proszę wykonać następujące kroki:<br />
# wygenerować &bdquo;prawie&rdquo; ciągłą sinusoidę (częstość 1, czas trwania 1 próbkowanie co 0,001)<br />
# spróbkować tą sinusoidę co 0,1 (proszę zastosować wycinki)<br />
# proszę skwantować spróbkowane wartości <br />
#Proszę wykreślić na jednym rysunku <br />
#* oryginalny sygnał <br />
#* sygnał spróbkowany w czasie<br />
#* sygnał spróbkowany w czasie i o skwantowanej amplitudzie (skorzystać z funkcji <tt>py.step</tt>)<br />
<br />
{{ Wyjaśnienie|title= wskazówka: Kwantowanie | text = <br />
<source lang =python><br />
N_bits = 3<br />
zakres = 2.0<br />
dy = zakres/2**N_bits<br />
s1_kwantowany = np.floor(s1_samp/dy)*dy +0.5*dy<br />
</source><br />
}}<br />
<br />
<!--<br />
# -*- coding: utf-8 -*-<br />
<br />
import pylab as py<br />
import numpy as np<br />
<br />
t = np.arange(0,1,0.001) # czas 'prawie ciągły'<br />
f1 = 1 # częstość sygnału 1<br />
<br />
s1 = np.sin(2*np.pi*f1*t) # prawie ciągły sygnał o częstości f1<br />
T_samp = t[0::100]<br />
s1_samp = s1[0::100]<br />
<br />
N_bits = 3<br />
zakres = 2.0<br />
dy = zakres/2**N_bits<br />
s1_kwantowany = np.floor(s1_samp/dy)*dy +0.5*dy<br />
ax = py.subplot(111)<br />
py.plot(t,s1)<br />
py.plot(T_samp,s1_samp,'ko',markersize = 5)<br />
py.step(T_samp,s1_kwantowany,where = 'post')<br />
<br />
ax.set_xticks(T_samp)<br />
poziomy = np.arange(-1+0.5*dy,1,dy)<br />
ax.set_yticks(poziomy)<br />
py.grid('on')<br />
<br />
<br />
for i in xrange(0,8):<br />
py.text(1.02,poziomy[i],bin(i))<br />
--><br />
<br />
==Sygnały testowe==<br />
<br />
<br />
=== Generowanie sygnałów testowych ===<br />
<br />
Do badania różnych metod analizy sygnałów potrzebne nam będą sygnały o znanych własnościach. W szczególności dobrze jest umnieć nadać sygnałom występującym w postaci cyfrowej, oraz sztucznym sygnałom próbnym pewne własności fizyczne takie jak:<br />
<br />
<ul><br />
<br />
<li><br />
czętość próbkowania<br />
<br />
<br />
<li><br />
czas trwania<br />
<br />
<br />
<li><br />
amplituda<br />
<br />
</ul><br />
=== Przykład sinus===<br />
Sinus o zadanej częstości (w Hz), długości trwania, częstości próbkowania i fazie.<br />
Poniższy kod implementuje i testuje funkcję <br />
:<math> \sin(f,T,Fs,\phi) = \sin(2*\pi f t)</math> dla <math>t \in \{0,T\}</math> <br />
<br />
<source lang= python><br />
# -*- coding: utf-8 -*-<br />
<br />
import pylab as py<br />
import numpy as np<br />
<br />
def sin(f = 1, T = 1, Fs = 128, phi =0 ):<br />
'''sin o zadanej częstości (w Hz), długości, fazie i częstości próbkowania<br />
Domyślnie wytwarzany jest sygnał reprezentujący <br />
1 sekundę sinusa o częstości 1 Hz i zerowej fazie próbkowanego 128 Hz<br />
'''<br />
<br />
dt = 1.0/Fs<br />
t = np.arange(0,T,dt)<br />
s = np.sin(2*np.pi*f*t + phi)<br />
return (s,t)<br />
<br />
<br />
(s,t) = sin(f=10,Fs=1000)<br />
py.plot(t,s)<br />
py.show()<br />
</source><br />
=== Przykład: eksport sygnału do pliku binarnego ===<br />
* Poniższy kod ilustruje sposób zapisu dwóch funkcji sinus o częstościach 10 Hz i 21 Hz do pliku binarnego:<br />
<br />
<source lang= python><br />
# -*- coding: utf-8 -*-<br />
<br />
import numpy as np<br />
<br />
T = 5<br />
Fs = 128.0<br />
<br />
s1 = sin(f=10, T=T, Fs=Fs)<br />
s2 = sin(f=21, T=T, Fs=Fs) <br />
<br />
signal = np.zeros((T*Fs, 2), dtype='<f')<br />
signal[:, 0] = s1<br />
signal[:, 1] = s2<br />
<br />
with open('test_signal.bin', 'wb') as f:<br />
signal.tofile(f)<br />
<br />
</source><br />
<br />
===Przykład: wczytanie sygnału do Svaroga ===<br />
W celu wczytania zapisanego binarnie sygnału do programu Svarog, po wybraniu File -> Open signal, należy wprowadzić częstość próbkowania sygnału oraz liczbę kanałów. <br />
<br />
[[Plik:svarog_open_signal.png|center|800px|thumb|<figure id="uid3" />Wczytywanie sygnału w programie Svarog.]]<br />
<br />
=== Delta ===<br />
Podobnie można zdefiniować funkcję delta o zadanym czasie trwania, częstości próbkowania i momencie wystąpienia impulsu:<br />
<br />
:<math> \delta(t_0) = \left\{^{1 \quad t=t_0} _{0 \quad t \ne t_0} \right.</math><br />
<br />
<source lang= python><br />
def delta(t0=0.5, T=1 ,Fs = 128):<br />
dt = 1.0/Fs<br />
t = np.arange(0,T,dt)<br />
d = np.zeros(len(t))<br />
d[np.ceil(t0*Fs)]=1<br />
return (d,t)<br />
</source><br />
<br />
=== Zadanie:===<br />
Analogicznie do powyższych przykładów proszę zaimplementować i przetestować funkcje generujące:<br />
* funkcję Gabora (funkcja Gaussa modulowana cosinusem) o zadanej częstości i standardowym odchyleniu w czasie, momencie wystąpienia, długości, częstości próbkowania i fazie. <br />
:<math> g = \exp\left(-\frac{1}{2}\left(\frac{t-t_0}{\sigma}\right)^2 \right) \cdot \cos(2 \pi f t + \phi); </math><br />
<br />
* szum gaussowski o zadanej średniej, odchyleniu standardowym, długości i częstości próbkowania.<br />
<br />
<!-- * pochodną funkcji Gaussa<br />
<br />
* połówkę funkcji Gaussa<br />
--></div>Prozanskihttp://brain.fuw.edu.pl/edu/index.php?title=STATLAB/Zadanie_zaliczeniowe5&diff=4479STATLAB/Zadanie zaliczeniowe52015-12-22T01:08:31Z<p>Prozanski: </p>
<hr />
<div><br />
'''Przepróbkowanie sygnału w dół'''<br />
<br />
===Przygotowanie danych===<br />
Wygneruj sygnał o długości 2 sekund i przebiegu zadanym wzorem:<br />
<br />
<math><br />
x(t) = \cos(2 \pi f_1 \, t) + 5 \cos(2 \pi f_2 \, t)<br />
</math><br />
<br />
gdzie: <math>f_1</math>=5 Hz, <math>f_2</math>=450 Hz. Sygnał wygeneruj z częstością próbkowania 192 000 Hz.<br />
<br />
===Przepróbkowanie sygnału===<br />
Dokonaj przebrókowania sygnału do częstosci 750 Hz (repróbkowanie 256 razy), przeprowadzając je na dwa sposoby:<br />
<br />
====Metoda A: jednokrotne====<br />
Repróbkowania dokonaj w jednym kroku. W tym celu zaprojektuj filtr typu butter rzędu 7. Częstość odcięcia ustaw w częstosci Nyquista docelowej częstości próbkowania. Przefiltruj oryginalny sygnał, a następnie wybierz z przefiltrowanego sygnału co 256 próbkę.<br />
<br />
====Metoda B: Wielokrotne====<br />
Repróbkowania dokonaj wielokrotnie (najlepiej w pętli) w następujących krokach:<br />
# Zaprojektuj filtr typu butter i częstości odcięcia '''k'''-krotnie niższej niż bieżąca częstości próbkowania i takim samym rzędzie jak w metodzie jednokrotnej. Rozważ '''k''' = 2 oraz '''k''' = 4. <br />
# Przefiltruj sygnał zaprojektowanym filtrem.<br />
# Wybierz z przefiltrowanego na końcu sygnału co '''k'''-tą próbkę.<br />
<br />
Powtórz powtórz powyższą procedurę z odpowiednimi częstosciami próbkowania i częstościami odcięcia aż uzyskasz częstość próbkowania 750 Hz.<br />
<br />
''Przykład: ''<br />
* ''dla '''k'''=2 w pierwszym wywołaniu tej procedury częstość próbkowania wynosi 192 000 Hz, częstość odcięcia 96 000 Hz''<br />
* ''dla '''k'''=4 w pierwszym wywołaniu tej procedury częstość próbkowania wynosi 192 000 Hz, częstość odcięcia 48 000 Hz''<br />
<br />
===Prezentacja wyniku i analiza ===<br />
Narysuj:<br />
# Charakterystyki amplitudowe i fazowe zaprojektowanych w metodzie A i B filtrów.<br />
# Narysuj sygnał oryginalny oraz sygnały repróbkowane metodami A i B.<br />
# Narysuj widma amplitudowe sygnałów repróbkowanych metodami A i B.<br />
# Która metoda repróbkowania dała lepszy wynik ? Uzasadnij wybór.<br />
# Czy inny rząd filtru zmienia wynik?</div>Prozanskihttp://brain.fuw.edu.pl/edu/index.php?title=STATLAB/Zadanie_zaliczeniowe5&diff=4478STATLAB/Zadanie zaliczeniowe52015-12-22T01:07:04Z<p>Prozanski: typografia we wzorach</p>
<hr />
<div><br />
'''Przepróbkowanie sygnału w dół'''<br />
<br />
===Przygotowanie danych===<br />
Wygneruj sygnał o długości 2 sekund i przebiegu zadanym wzorem:<br />
<br />
<math><br />
x(t) = cos(2 \pi f_1 t) + 5*cos(2 \pi f_2 t)<br />
</math><br />
<br />
gdzie: <math>f_1</math>=5 Hz, <math>f_2</math>=450 Hz. Sygnał wygeneruj z częstością próbkowania 192 000 Hz. <br />
<br />
===Przepróbkowanie sygnału===<br />
Dokonaj przebrókowania sygnału do częstosci 750 Hz (repróbkowanie 256 razy), przeprowadzając je na dwa sposoby:<br />
<br />
====Metoda A: jednokrotne====<br />
Repróbkowania dokonaj w jednym kroku. W tym celu zaprojektuj filtr typu butter rzędu 7. Częstość odcięcia ustaw w częstosci Nyquista docelowej częstości próbkowania. Przefiltruj oryginalny sygnał, a następnie wybierz z przefiltrowanego sygnału co 256 próbkę.<br />
<br />
====Metoda B: Wielokrotne====<br />
Repróbkowania dokonaj wielokrotnie (najlepiej w pętli) w następujących krokach:<br />
# Zaprojektuj filtr typu butter i częstości odcięcia '''k'''-krotnie niższej niż bieżąca częstości próbkowania i takim samym rzędzie jak w metodzie jednokrotnej. Rozważ '''k''' = 2 oraz '''k''' = 4. <br />
# Przefiltruj sygnał zaprojektowanym filtrem.<br />
# Wybierz z przefiltrowanego na końcu sygnału co '''k'''-tą próbkę.<br />
<br />
Powtórz powtórz powyższą procedurę z odpowiednimi częstosciami próbkowania i częstościami odcięcia aż uzyskasz częstość próbkowania 750 Hz.<br />
<br />
''Przykład: ''<br />
* ''dla '''k'''=2 w pierwszym wywołaniu tej procedury częstość próbkowania wynosi 192 000 Hz, częstość odcięcia 96 000 Hz''<br />
* ''dla '''k'''=4 w pierwszym wywołaniu tej procedury częstość próbkowania wynosi 192 000 Hz, częstość odcięcia 48 000 Hz''<br />
<br />
===Prezentacja wyniku i analiza ===<br />
Narysuj:<br />
# Charakterystyki amplitudowe i fazowe zaprojektowanych w metodzie A i B filtrów.<br />
# Narysuj sygnał oryginalny oraz sygnały repróbkowane metodami A i B.<br />
# Narysuj widma amplitudowe sygnałów repróbkowanych metodami A i B.<br />
# Która metoda repróbkowania dała lepszy wynik ? Uzasadnij wybór.<br />
# Czy inny rząd filtru zmienia wynik?</div>Prozanskihttp://brain.fuw.edu.pl/edu/index.php?title=Analiza_sygna%C5%82%C3%B3w_-_wyk%C5%82ad&diff=4371Analiza sygnałów - wykład2015-10-25T22:19:16Z<p>Prozanski: podmiana linka</p>
<hr />
<div>[[Category:Przedmioty specjalizacyjne]]<br />
=Analiza Sygnałów=<br />
<br />
# [[wstep | Wstęp]]<br />
<br />
Klasyczna analiza sygnałów -- zob. też http://durka.info/ksiazki/as/as_klasyczna.pdf<br />
<br />
#[[Szereg Fouriera|Szereg Fouriera]]<br />
#[[Przekształcenie Fouriera|Przekształcenie Fouriera]]<br />
#[[Systemy liniowe niezmiennicze w czasie (LTI)|Systemy liniowe niezmiennicze w czasie]]<br />
#[[Twierdzenia o splocie i o próbkowaniu (aliasing)|Twierdzenie o splocie]]<br />
#[[Twierdzenie Wienera-Chinczyna|Twierdzenie Wienera-Chinczyna]]<br />
#[[Aliasing|Aliasing]]<br />
#[[Twierdzenie o próbkowaniu|Twierdzenie o próbkowaniu]]<br />
#[[Funkcja_systemu|Funkcja systemu]]<br />
#[[Model autoregresyjny (AR)|Model autoregresyjny]]<br />
<!-- ##[[Procesy_stochastyczne|Procesy stochastyczne]] --><br />
Pomiędzy czasem a częstością<br />
# [[Zasada nieoznaczoności|Zasada nieoznaczoności]]<br />
# [[Transformata Wignera|Transformata Wignera]]<br />
# [[Spektrogram|Spektrogram — oknowana transformata Fouriera]]<br />
# [[Falki (wavelets)|Falki]]<br />
# [[Reprezentacje czas-częstość|Reprezentacje czas-częstość]]<br />
# [[Reprezentacje przybliżone | Przybliżenia adaptacyjne i algorytm matching pursuit]]<br />
Inne<br />
#[[Analiza sygnałów wielowymiarowych|Analiza sygnałów wielozmiennych]]<br />
<br />
Ciekawostki<br />
#[[Sztuczne sieci neuronowe (ANN )|Sztuczne sieci neuronowe]]<br />
#[[Algorytmy Genetyczne|Algorytmy Genetyczne]]<br />
<!-- # [[appendix | Dodatek ]]--><br />
<br />
[http://www.dropbox.com/sh/7b0jzdx7hbv0z13/AAD8KKySea36vGcpZYe1uSeya?dl=0 Slajdy z wykładów 2014/15]<br />
<br />
----<br />
<!--<br />
[[ORGANIZACJA_FUW_AS|Sprawy organizacyjne: materiały do powtórzenia przed egzaminem, organizacja egzaminu, wyniki...]]<br />
--><br />
<br />
Zapraszamy do pobrania aktualnej wersji omawianego na wykładzie i ćwiczeniach narzędzia do eksperymentowania z metodami analizy sygnałów, [http://www.mimuw.edu.pl/~ptr/fid/ czyli programu SVAROG]<br />
----<br />
<br />
{{color|green|'''Całość podręcznika jest udostępniona na licencji [http://creativecommons.org/licenses/by-sa/3.0/pl Creative Commons Uznanie autorstwa-Na tych samych zasadach 3.0 Polska].'''}} [[Grafika:CC-88x31.png]]<br />
Autor: [http://durka.name Piotr Durka]. Podręcznik powstał częściowo w oparciu o skrypty udostępniane wcześniej na stronach http://www.fuw.edu.pl/~durka/ksiazki/as i książkę [http://www.fuw.edu.pl/~durka/ksiazki/UNI.html MP an unification in EEG analysis]</div>Prozanskihttp://brain.fuw.edu.pl/edu/index.php?title=ZasadyZaliczenia&diff=4301ZasadyZaliczenia2015-10-24T12:26:08Z<p>Prozanski: dopisanie informacji o pierwszym projekcie</p>
<hr />
<div>== Zasady zaliczenia ćwiczeń z [[Analiza_sygnałów - exercises|Analizy Sygnałów]] ==<br />
<br />
W trakcie semestru odbędą się dwa kolokwia. Z każdego z nich będzie można uzyskać 10 punktów, z czego do zaliczenia przedmiotu konieczne będzie zdobycie łącznie 10 punktów z kolokwiów (czyli 50%). Kolokwia będą wspólne dla wszystkich grup, będą się odbywać w poniedziałki przed południem.<br />
<br />
Wymagane będzie wykonanie dwóch projektów indywidualnych. Jeden z nich będzie wymagany w połowie semestru, a drugi do końca semestru. Nie będzie możliwości oddawania projektów po zakończeniu zajęć semestru zimowego. Z każdego projektu będzie można uzyskać 10 punktów, z czego do zaliczenia przedmiotu konieczne jest zdobycie łącznie 10 punktów z projektów (czyli 50%).<br />
<br />
# Pierwszy projekt [[STATLAB/Zadanie zaliczeniowe3|Analiza czasowa sygnału audio]] należy wykonać i oddać '''do końca listopada'''. Dokładniejsze terminy zaliczenia przekażą poszczególni prowadzący grup.<br />
# Drugi projekt zostanie ogłoszony niebawem.<br />
<br />
Na większości (10) zajęć będą miały miejsce kartkówki na podstawie materiału z zajęć poprzednich. Pytania nie będą wymagające, będą raczej sprawdzać uwagę i systematyczną pracę. Łącznie za kartkówki można uzyskać 10 punktów (maksymalnie 1 punkt za każdą), natomiast przy zaliczeniu przedmiotu dozwolone jest niezaliczenie maksymalnie trzech kartkówek (nieusprawiedliwiona nieobecność również pociąga za sobą niezaliczenie kartkówki). Punktów z kartkówek nie można poprawiać.<br />
<br />
Na podstawie łącznej liczby punktów (maksymalnie do zdobycia jest 50) zostanie obliczona ocena z ćwiczeń:<br />
# [25–30) pkt: 3.0 (dst)<br />
# [30–35) pkt: 3.5 (dst+)<br />
# [35–40) pkt: 4.0 (db)<br />
# [40–45) pkt: 4.5 (db+)<br />
# [45–50] pkt: 5.0 (bdb)<br />
<br />
Dozwolone są maksymalnie dwie nieusprawiedliwione nieobecności.</div>Prozanskihttp://brain.fuw.edu.pl/edu/index.php?title=%C4%86wiczenia_1&diff=4255Ćwiczenia 12015-10-05T22:16:03Z<p>Prozanski: info o Svarogu i Javie</p>
<hr />
<div>== Dokumentacja modułu scipy.signal ==<br />
Proszę zapoznać się z dokumentacją biblioteki scipy.signal:<br />
<br />
http://docs.scipy.org/doc/scipy-0.14.0/reference/signal.html<br />
<br />
== Svarog: uruchamianie i konfiguracja ==<br />
Zajęcia będziemy ilustrować przykładami realizowanymi w programie Svarog, który można pobrać [http://www.fuw.edu.pl/~durka/svarog-1.0.10.zip stąd]. Po rozpakowaniu paczki do dowolnego katalogu należy uruchomić skrypt „run-svarog.sh”.<br />
<br />
W przypadku pracy na własnych komputerach, do prawidłowego uruchomienia pluginu do analizy sygnałów, z którego będziemy korzystać w dalszej części ćwiczeń, konieczne jest zainstalowanie środowiska Oracle Java SE w wersji 8, które można pobrać [http://www.oracle.com/technetwork/java/javase/downloads/index.html ze strony wydawcy]. Alternatywnie, użytkownicy systemu Ubuntu lub pokrewnych dystrybucji mogą zainstalować środowisko Java według instrukcji [http://www.webupd8.org/2012/09/install-oracle-java-8-in-ubuntu-via-ppa.html dostępnych na tej stronie].<br />
<br />
==Sygnały ciągłe i dyskretne ==<br />
===Próbkowanie w czasie ===<br />
Proszę powtórzyć sobie pojęcia:<br />
* częstość próbkowania<br />
* częstość Nyquista<br />
* aliasing<br />
<br />
W poniższym ćwiczeniu chcemy zbadać efekt próbkowania sygnału w czasie. W komputerach nie mamy dostępu do sygnału ciągłego. Na nasze potrzeby wygenerujemy sygnały próbkowane z bardzo dużą częstością, które będą dla nas aproksymacją sygnałów ciągłych. Przy ich pomocy zaprezentujemy efekt utożsamiania (aliasingu).<br />
<br />
Proszę wytworzyć wektor reprezentujący czas &bdquo;prawie&rdquo; ciągły. Będzie to u nas 1000 wartości z przedziału [0,1) wziętych z odstępem 0,001.<br />
<br />
<source lang =python><br />
# -*- coding: utf-8 -*-<br />
<br />
import pylab as py<br />
import numpy as np<br />
<br />
t = np.arange(0,1,0.001) # czas 'prawie ciągły'<br />
</source><br />
<br />
Teraz proszę wygenerować dwie sinusoidy: jedną o częstości <tt>-1</tt> a drugą o częstości <tt>9</tt>. Dla przypomnienia wyrażenie:<br />
:<math> s(t) = \sin(2 \pi f t)</math> możemy w pythonie zapisać:<br />
<source lang =python><br />
s = np.sin(2*np.pi*f*t)<br />
</source><br />
<br />
Proszę [[TI/Matplotlib#Kilka_wykres.C3.B3w_we_wsp.C3.B3lnych_osiach|wykreślić]] obie sinusoidy.<br />
<br />
Teraz proszę spróbkować czas i nasze &bdquo;prawie&rdquo; ciągłe sinusoidy z okresem próbkowania 0,1. (Trzeba pobrać co 100 element, proszę posłużyć się [[Programowanie_z_Pythonem/Sekwencje#Wycinki|wycinkami]]) <!-- [[TI/Sekwencje|Struktury danych — sekwencje|wycinkami]]. --><br />
Na tle &bdquo;prawie&rdquo; ciągłych sinusoid proszę dorysować punkty ze spróbkowanych sygnałów. Aby punkty były dobrze widoczne proponuję użyć markerów <tt>x</tt> oraz <tt>+</tt>.<br />
<br />
Proszę zaobserwować wzajemne położenie punktów. Czy można odróżnić sinusoidę o częstości &minus;1 od sinusoidy o częstości 9, jeśli obie są próbkowane z częstością 10? Jak można uogólnić tą obserwację?<br />
<tt>*</tt><br />
<!--<br />
# -*- coding: utf-8 -*-<br />
<br />
import pylab as py<br />
import numpy as np<br />
<br />
t = np.arange(0,1,0.001) # czas 'prawie ciągły'<br />
f1 = -1 # częstość sygnału 1<br />
f2 = 9 # częstość sygnału 2<br />
s1 = np.sin(2*np.pi*f1*t) # prawie ciągły sygnał o częstości f1<br />
s2 = np.sin(2*np.pi*f2*t) # prawie ciągły sygnał o częstości f2<br />
<br />
T = 0.1 # okres próbkowania<br />
Fs = 1/T # częstość próbkowania<br />
FN = Fs/2 # częstość Nyquista<br />
T_samp = t[0::100] # czas pobierania próbek<br />
s1_samp = s1[0::100] # próbkowany sygnał o częstości f1<br />
s2_samp = s2[0::100] # próbkowany sygnał o częstości f2<br />
<br />
py.plot(t, s1, 'g')<br />
py.plot(t, s2, 'b')<br />
py.plot(T_samp, s1_samp, 'gx', markersize=10)<br />
py.plot(T_samp, s2_samp, 'r+', markersize=10)<br />
--><br />
<br />
<!--<br />
===Błąd kwantyzacji ===<br />
Kiedy mierzymy fizyczne wielkości w celu dalszej analizy najczęściej chcemy przypisać im pewne liczby. Liczby w systemach cyfrowych reprezentowane są ze skończoną dokładnością. Urządzenia dokonujące przypisania liczby do mierzonej wartości to przetworniki analogowo-cyfrowe (ang. ADC, Analog to Digital Converter).<br />
Charakteryzują się one określoną ilością bitów ''N'', za pomocą których reprezentują liczby. Pełen zakres wartości pomiarowych ''R'' jest dzielony na <math>2^N</math> poziomów. Błąd kwantyzacji szacujemy jako nie większy niż <math> \frac{R}{2^{N+1}} </math>. <br />
<br />
====Zadanie:====<br />
Proszę oszacować błąd kwantyzacji sygnału EEG mierzonego za pomocą 12-bitowego konwertera. Zakres pomiarowy tego urządzenia to &plusmn;200 &mu;V.<br />
====Zadanie:====<br />
Proszę zilustrować efekt kwantyzacji dla trzybitowego przetwornika o zakresie 2.<br />
W tym celu proszę wykonać następujące kroki:<br />
# wygenerować &bdquo;prawie&rdquo; ciągłą sinusoidę (częstość 1, czas trwania 1 próbkowanie co 0,001)<br />
# spróbkować tą sinusoidę co 0,1 (proszę zastosować wycinki)<br />
# proszę skwantować spróbkowane wartości <br />
#Proszę wykreślić na jednym rysunku <br />
#* oryginalny sygnał <br />
#* sygnał spróbkowany w czasie<br />
#* sygnał spróbkowany w czasie i o skwantowanej amplitudzie (skorzystać z funkcji <tt>py.step</tt>)<br />
<br />
{{ Wyjaśnienie|title= wskazówka: Kwantowanie | text = <br />
<source lang =python><br />
N_bits = 3<br />
zakres = 2.0<br />
dy = zakres/2**N_bits<br />
s1_kwantowany = np.floor(s1_samp/dy)*dy +0.5*dy<br />
</source><br />
}}<br />
<br />
<!--<br />
# -*- coding: utf-8 -*-<br />
<br />
import pylab as py<br />
import numpy as np<br />
<br />
t = np.arange(0,1,0.001) # czas 'prawie ciągły'<br />
f1 = 1 # częstość sygnału 1<br />
<br />
s1 = np.sin(2*np.pi*f1*t) # prawie ciągły sygnał o częstości f1<br />
T_samp = t[0::100]<br />
s1_samp = s1[0::100]<br />
<br />
N_bits = 3<br />
zakres = 2.0<br />
dy = zakres/2**N_bits<br />
s1_kwantowany = np.floor(s1_samp/dy)*dy +0.5*dy<br />
ax = py.subplot(111)<br />
py.plot(t,s1)<br />
py.plot(T_samp,s1_samp,'ko',markersize = 5)<br />
py.step(T_samp,s1_kwantowany,where = 'post')<br />
<br />
ax.set_xticks(T_samp)<br />
poziomy = np.arange(-1+0.5*dy,1,dy)<br />
ax.set_yticks(poziomy)<br />
py.grid('on')<br />
<br />
<br />
for i in xrange(0,8):<br />
py.text(1.02,poziomy[i],bin(i))<br />
--><br />
<br />
==Sygnały testowe==<br />
<br />
<br />
=== Generowanie sygnałów testowych ===<br />
<br />
Do badania różnych metod analizy sygnałów potrzebne nam będą sygnały o znanych własnościach. W szczególności dobrze jest umnieć nadać sygnałom występującym w postaci cyfrowej, oraz sztucznym sygnałom próbnym pewne własności fizyczne takie jak:<br />
<br />
<ul><br />
<br />
<li><br />
czętość próbkowania<br />
<br />
<br />
<li><br />
czas trwania<br />
<br />
<br />
<li><br />
amplituda<br />
<br />
</ul><br />
=== Przykład sinus===<br />
Sinus o zadanej częstości (w Hz), długości trwania, częstości próbkowania i fazie.<br />
Poniższy kod implementuje i testuje funkcję <br />
:<math> \sin(f,T,Fs,\phi) = \sin(2*\pi f t)</math> dla <math>t \in \{0,T\}</math> <br />
<br />
<source lang= python><br />
# -*- coding: utf-8 -*-<br />
<br />
import pylab as py<br />
import numpy as np<br />
<br />
def sin(f = 1, T = 1, Fs = 128, phi =0 ):<br />
'''sin o zadanej częstości (w Hz), długości, fazie i częstości próbkowania<br />
Domyślnie wytwarzany jest sygnał reprezentujący <br />
1 sekundę sinusa o częstości 1 Hz i zerowej fazie próbkowanego 128 Hz<br />
'''<br />
<br />
dt = 1.0/Fs<br />
t = np.arange(0,T,dt)<br />
s = np.sin(2*np.pi*f*t + phi)<br />
return (s,t)<br />
<br />
<br />
(s,t) = sin(f=10,Fs=1000)<br />
py.plot(t,s)<br />
py.show()<br />
</source><br />
=== Przykład: eksport sygnału do pliku binarnego ===<br />
* Poniższy kod ilustruje sposób zapisu dwóch funkcji sinus o częstościach 10 Hz i 21 Hz do pliku binarnego:<br />
<br />
<source lang= python><br />
# -*- coding: utf-8 -*-<br />
<br />
import numpy as np<br />
<br />
if __name__ == "__main__":<br />
T = 5<br />
Fs = 128.0<br />
<br />
s1 = sin(f=10, T=T, Fs=Fs)<br />
s2 = sin(f=21, T=T, Fs=Fs) <br />
<br />
signal = np.zeros((T*Fs, 2), dtype='<f')<br />
signal[:, 0] = s1<br />
signal[:, 1] = s2<br />
<br />
with open('test_signal.bin', 'wb') as f:<br />
signal.tofile(f)<br />
<br />
</source><br />
<br />
===Przykład: wczytanie sygnału do Svaroga ===<br />
W celu wczytania zapisanego binarnie sygnału do programu Svarog, po wybraniu File -> Open signal, należy wprowadzić częstość próbkowania sygnału oraz liczbę kanałów. <br />
<br />
[[Plik:svarog_open_signal.png|center|800px|thumb|<figure id="uid3" />Wczytywanie sygnału w programie Svarog.]]<br />
<br />
=== Delta ===<br />
Podobnie można zdefiniować funkcję delta o zadanym czasie trwania, częstości próbkowania i momencie wystąpienia impulsu:<br />
<br />
:<math> \delta(t_0) = \left\{^{1 \quad t=t_0} _{0 \quad t \ne t_0} \right.</math><br />
<br />
<source lang= python><br />
def delta(t0=0.5, T=1 ,Fs = 128):<br />
dt = 1.0/Fs<br />
t = np.arange(0,T,dt)<br />
d = np.zeros(len(t))<br />
d[np.ceil(t0*Fs)]=1<br />
return (d,t)<br />
</source><br />
<br />
=== Zadanie:===<br />
Analogicznie do powyższych przykładów proszę zaimplementować i przetestować funkcje generujące:<br />
* funkcję Gabora (funkcja Gaussa modulowana cosinusem) o zadanej częstości i standardowym odchyleniu w czasie, momencie wystąpienia, długości, częstości próbkowania i fazie. <br />
:<math> g = \exp\left(-\frac{1}{2}\left(\frac{t-t_0}{\sigma}\right)^2 \right) \cdot \cos(2 \pi f t + \phi); </math><br />
<br />
* szum gaussowski o zadanej średniej, odchyleniu standardowym, długości i częstości próbkowania.<br />
<br />
<!-- * pochodną funkcji Gaussa<br />
<br />
* połówkę funkcji Gaussa<br />
--></div>Prozanskihttp://brain.fuw.edu.pl/edu/index.php?title=ZasadyZaliczenia&diff=4164ZasadyZaliczenia2015-10-01T17:06:35Z<p>Prozanski: tytuł</p>
<hr />
<div>== Zasady zaliczenia ćwiczeń z Analizy Sygnałów ==<br />
<br />
W trakcie semestru odbędą się dwa kolokwia. Z każdego z nich będzie można uzyskać 10 punktów, z czego do zaliczenia przedmiotu konieczne będzie zdobycie łącznie 10 punktów z kolokwiów (czyli 50%). Kolokwia będą wspólne dla wszystkich grup, będą się odbywać w poniedziałki przed południem.<br />
<br />
Wymagane będzie wykonanie dwóch projektów indywidualnych. Jeden z nich będzie wymagany w połowie semestru, a drugi do końca semestru. Nie będzie możliwości oddawania projektów po zakończeniu zajęć semestru zimowego. Z każdego projektu będzie można uzyskać 10 punktów, z czego do zaliczenia przedmiotu konieczne jest zdobycie łącznie 10 punktów z projektów (czyli 50%).<br />
<br />
Na większości (10) zajęć będą miały miejsce kartkówki na podstawie materiału z zajęć poprzednich. Pytania nie będą wymagające, będą raczej sprawdzać uwagę i systematyczną pracę. Łącznie za kartkówki można uzyskać 10 punktów (maksymalnie 1 punkt za każdą), natomiast przy zaliczeniu przedmiotu dozwolone jest niezaliczenie maksymalnie trzech kartkówek (nieusprawiedliwiona nieobecność również pociąga za sobą niezaliczenie kartkówki). Punktów z kartkówek nie można poprawiać.<br />
<br />
Na podstawie łącznej liczby punktów (maksymalnie do zdobycia jest 50) zostanie obliczona ocena z ćwiczeń:<br />
# [25–30) pkt: 3.0 (dst)<br />
# [30–35) pkt: 3.5 (dst+)<br />
# [35–40) pkt: 4.0 (db)<br />
# [40–45) pkt: 4.5 (db+)<br />
# [45–50] pkt: 5.0 (bdb)<br />
<br />
Dozwolone są maksymalnie dwie nieusprawiedliwione nieobecności.</div>Prozanskihttp://brain.fuw.edu.pl/edu/index.php?title=ZasadyZaliczenia&diff=4163ZasadyZaliczenia2015-10-01T17:05:02Z<p>Prozanski: pierwsza wersja zasad zaliczenia</p>
<hr />
<div>W trakcie semestru odbędą się dwa kolokwia. Z każdego z nich będzie można uzyskać 10 punktów, z czego do zaliczenia przedmiotu konieczne będzie zdobycie łącznie 10 punktów z kolokwiów (czyli 50%). Kolokwia będą wspólne dla wszystkich grup, będą się odbywać w poniedziałki przed południem.<br />
<br />
Wymagane będzie wykonanie dwóch projektów indywidualnych. Jeden z nich będzie wymagany w połowie semestru, a drugi do końca semestru. Nie będzie możliwości oddawania projektów po zakończeniu zajęć semestru zimowego. Z każdego projektu będzie można uzyskać 10 punktów, z czego do zaliczenia przedmiotu konieczne jest zdobycie łącznie 10 punktów z projektów (czyli 50%).<br />
<br />
Na większości (10) zajęć będą miały miejsce kartkówki na podstawie materiału z zajęć poprzednich. Pytania nie będą wymagające, będą raczej sprawdzać uwagę i systematyczną pracę. Łącznie za kartkówki można uzyskać 10 punktów (maksymalnie 1 punkt za każdą), natomiast przy zaliczeniu przedmiotu dozwolone jest niezaliczenie maksymalnie trzech kartkówek (nieusprawiedliwiona nieobecność również pociąga za sobą niezaliczenie kartkówki). Punktów z kartkówek nie można poprawiać.<br />
<br />
Na podstawie łącznej liczby punktów (maksymalnie do zdobycia jest 50) zostanie obliczona ocena z ćwiczeń:<br />
# [25–30) pkt: 3.0 (dst)<br />
# [30–35) pkt: 3.5 (dst+)<br />
# [35–40) pkt: 4.0 (db)<br />
# [40–45) pkt: 4.5 (db+)<br />
# [45–50] pkt: 5.0 (bdb)<br />
<br />
Dozwolone są maksymalnie dwie nieusprawiedliwione nieobecności.</div>Prozanskihttp://brain.fuw.edu.pl/edu/index.php?title=Analiza_sygna%C5%82%C3%B3w_-_%C4%87wiczenia&diff=4162Analiza sygnałów - ćwiczenia2015-10-01T16:46:18Z<p>Prozanski: link do zasad zaliczenia</p>
<hr />
<div>[[Category:Przedmioty specjalizacyjne]]<br />
<!--<br />
#[[Systemy liniowe niezmiennicze w czasie | Systemy liniowe niezmiennicze w czasie]]<br />
[[Ćwiczenia 8|Filtrowanie obrazów]]<br />
[[Ćwiczenia 9|Analiza czas-częstość]]<br />
[[Ćwiczenia 10|Analiza czas-częstość &mdash; STFT i transformata falkowa]]<br />
[[Ćwiczenia 11|Analiza czas-częstość &mdash; reprezentacje energetyczne]]<br />
--><br />
<br />
[[ZasadyZaliczenia|Zasady zaliczenia ćwiczeń]]<br />
<br />
#[[Ćwiczenia 1|Sygnały]] -- wstęp, generacja sygnaów testowych, aliasing, przypomnienie Pythona, wprowadzenie Svaroga; jak wczytać do Svaroga z Pythona i w drugą stronę<br />
#[[Ćwiczenia 2|Transformata Fouriera]] i [[Ćwiczenia 2_2|Transformata Fouriera cd]] -- piszemy w Pythonie, porównujemy ze Svarogiem<br />
#[[Ćwiczenia 3|Okienkowanie sygnału i transformata Fouriera]], [[Nieparametryczne widmo mocy |Estymacja widma mocy w oparciu o transformatę Fouriera]] -- widmo średniej vs. średnie widmo<br />
#[[Ćwiczenia 6|Filtry I]] -- w Pythonie i w Svarogu<br />
#[[Ćwiczenia 7|Filtry II]] -- oglądanie charakterystyk filtrów, ogładanie widma sygnałów przefiltrowanych, efektywność filtrów, przesuwanie fazy<br />
#[[Ćwiczenia 4|Funkcja autokorelacji i procesy AR]]<br />
#[[Ćwiczenia 5|Procesy AR]] -- AR w Pyhonie + DTF w Svarogu (MK)<br />
#[[AS cwiczenia ICA| ICA]] -- montaże, ICA<br />
#[[AS cwiczeniaTF|Czas-częstość]] STFT i falki<br />
#[[AS cwiczeniaMP|Matching pursuit]] MP w Svarogu, zabawa parametrami dekompozycji, zabawa filtrowaniem map w Svarogu, postprocessing w Pythonie<br />
# uśrednianie gęstości energii vs gęstość energii uśrednionego sygnału: symulacje, ERD/S<br />
#[[AS cwiczeniaUNIFIKACJA]] -- ostatnie ćwiczenia na ktrych porównujemy na tych samych sygnałach rzeczywistych i symulowanych działanie różnych metod -- MMP vs DTF vs ICA, STFT vs WT vs MP itp, potrzebne fajne przykłady<br />
<br />
autor dr Jarosław Żygierewicz<br />
<br />
do uzupełnienia:<br />
# jakiś eksport wyników prcessingu ze Svaroga do plików czytanych w Pythonie (poza MP)?<br />
# funkcjonalności pluginów - coś jeszcze potrzebne?<br />
<br />
<!--<br />
== INFORMACJE DODATKOWE ==<br />
Materiały 2014/2015<br />
# [[kolokwia2013_2014|Informacje o zaliczeniu ćwiczeń]] <br />
# [[kolokwia2014_2015_kol1|Zagadnienia przygotowawcze do 1 kolokwium]] <br />
# [[Projekt2014|Projekt zaliczeniowy]]<br />
<br />
Materiały 2013/2014<br />
# [[Plik:Zadania_powtorzeniowe.pdf|Zadania powtorzeniowe do 2 kolokwium]]<br />
--></div>Prozanski