Tekst źródłowy strony Laboratorium EEG/CSP

[[Laboratorium_EEG]]/BSS
=Ślepa separacja źródeł=
Rozważmy N - kanałowy sygnał EEG.
Próbkę tego sygnału możemy przedstawić jako punkt w przestrzeni rozpiętej przez osie, z których każda reprezentuje wartość potencjału w jednym kanale. Cały sygnał tworzy w tej przestrzeni chmurę punktów. Rozciągłość tej chmury w danym kierunku mówi nam o wariancji (zmienności) sygnału w tym kierunku. 

Taki zbiór punktów wygodniej jest analizować w układzie współrzędnych zgodnym z osiami głównymi macierzy kowariancji.
W dalszej części rozważań założymy, że te przestrzenie, w których rozważamy sygnały są przestrzeniami wektorowymi, a pojedyncze próbki wielokanałowego sygnału są wektorami. 
[[Plik:Kowariancja.png|200px|center]]

==Filtry przestrzenne i ślepa separacja źródeł==
Sygnał EEG jest superpozycją aktywności elektrycznej wielu źródeł.
Jak można estymować aktywność samych źródeł?
[[Plik:Mieszanie.png|200px|center]]
Niech:
: <math>s(t)</math> - aktywność niezależnych źródeł,
: <math>x(t)</math> mierzony sygnał
: <math>A</math> macierz przejścia taka, że:
::<math>x(t) = A s(t)</math> (*)
:<math>s(t) = A^{-1}x(t) = P x(t)</math>
Macierz kowariancji dla sygnałów <math>x(t)</math> estymujemy tak:
:<math> C_x = E[x(t)x(t)^T]</math>
Podstawiając (*) mamy:
:<math> C_x = E[x x^T] = E[As(As)^T] = A E[s s^T] A^T = A C_s A^T</math>
Z założenia, że źródła są niezależne wynika, że macierz <math>C_s</math> jest diagonalna.
Przekształcając powyższe równanie możemy zapisać:
:<math>A^{-1} C_x (A^T)^{-1} = P C_x P^T = C_s</math>
Odwzorowanie <math>P = A^{-1}</math> diagonalizuje macierz <math>C_x</math>. 

Powyższe rozumowanie jest słuszne w przypadku gdy mamy do czynienia z sygnałem stacjonarnym, tzn. jego macierz kowariancji jest niezależna od czasu -> caly czas aktywna jest ta sama konfiguracja źródeł niezależnych.
W przypadku gdy tak nie jest to konstrukcję filtra przestrzennego można oprzeć o  jednoczesną diagonalizację macierzy kowariancji odpowiadających różnym stanom osoby badanej.

[[Plik:Diagonalizacja.png|200px|center]]



==Common Spatial Pattern ==
===Koncepcja===
Dla ustalenia uwagi możemy myśleć o eksperymencie wywołującym potencjał P300. Mamy w nim dwie sytuacje eksperymentalne. Oznaczmy (<math>T</math> - target) próby, w których pojawił się oczekiwany bodziec, zaś  (<math>NT</math> - non-target) gdy pojawił się bodziec standardowy.
Chcielibyśmy znaleźć taki montaż, czyli taką kombinację liniową kanałów, które maksymalizuje stosunek mocy (wariancji) sygnałów rejestrowanych w dwóch rożnych warunkach eksperymentalnych. 

===Formalizm===
Metoda ta polega na znalezieniu takiego kierunku <math>w</math> w przestrzeni sygnałów, że sygnał z warunku <math>T</math> rzutowany na ten kierunek ma dużą wariancje a sygnał z warunku <math>NT</math> ma wariancję małą. 

Rzutowanie sygnału <math>x(t)</math> na kierunek <math>w</math> odbywa się przez policzenie iloczynu skalarnego dla każdej chwili czasu <math>t</math>:
:<math> s_w(t) = w^T x(t)</math>
Wariancja tego rzutowanego sygnału to:
:<math> \mathrm{var}(s_w) = E[s_w s_w^T] = E[ w^T x (w^T x)^T] = w^T E[x x^T] w = w^T C_x w </math>
Zatem znalezienie właściwego kierunku rzutowania można wyrazić jako szukanie maksimum wyrażenia <math> J(w) </math>(jest to tzw. iloraz Rayleigh'a):
: <math>J(w) = \frac{w^T C_T w}{w^T C_{NT} w}  </math>
Ekstremum tego ilorazu można znaleźć poprzez policzenie gradientu <math>J(w)</math> i przyrównanie go do zera:
:<math> \nabla J(w) =  \frac{ 1  C_{T} w+w^T C_{T} 1}{w^T C_{NT} w}-\frac{w^T  C_{T} w\left( 1  C_{NT} w+w^T C_{NT} 1\right)}{\left(w^T  C_{NT} w\right)^2}</math>
ponieważ macierze kowariancji są symetryczne 
::<math>\nabla J(w) =   \frac{   1}{w^T  C_{NT} w}\left[    C_{T} w+ C_{T}w  -\frac{w^T  C_{T} w}{w^T  C_{NT} w} \left(   C_{NT} w+ C_{NT}w \right) \right]</math>
::<math>= \frac{   2}{w^T  C_{NT} w}\left[     C_{T}w  -\frac{w^T  C_{T} w}{w^T  C_{NT} w}   C_{NT} w \right]</math>
Przyrównując to wyrażenie do zera dostajemy:
:<math>       C_{T}w  =\frac{w^T  C_{T} w}{w^T  C_{NT} w}   C_{NT} w   </math>
Liczba <math> \lambda =    \frac{w^T  C_{T} w}{w^T  C_{NT} w}   C_{NT} w </math> jest uogólnioną wartością własną, zaś <math>w</math> jest uogólnionym wektorem własnym odpowiadającym tej wartości. 

Aby znaleźć <math> \lambda</math> i <math>w</math> wystarczy rozwiązać zagadnienie własne. W matlabie możemy w tym celu wykorzystać funkcję <tt>eig</tt>

<!--
Odwzorowanie to można przedstawić w postaci macierzy <math>P</math>, której każdy wiersz zawiera wagi dla odpowiednich kanałów. 
Macierz zawierająca sygnał <math>X^{\pm}(t)</math> 
ma wymiary <math>C \times  N</math>, gdzie 
<math>C</math> to liczba kanałów EEG, natomiast 
<math>N</math> to liczba próbek dla każdego z kanałów. 
Macierz <math>P</math> przekształca sygnał <math>X^{\pm}(t)</math>  zgodnie ze wzorem:
:<math>X^{\pm}_{CSP}(t)=P^T  X^{\pm}(t) </math>
Załóżmy dalej, że sygnały <math>X^{+}_{CSP} (t)</math> i <math>X^{-}_{CSP} (t)</math> są generowane przez niezależne procesy stochastyczne, tzn. spełnione są następujące warunki. 
# Syganły <math>X^{+}_{CSP} (t)</math> i <math>X^{-}_{CSP} (t)</math> są niezależne.
# Brak korelacji pomiędzy kanałami w sygnałach<math>X^{+}_{CSP} (t)</math> i <math>X^{-}_{CSP} (t)</math>.
# Przynajmniej dla jednego z kanałów wariancja przetransformowanego sygnału jest maksymalna przy wystąpieniu bodźca i minimalna przy jego braku.
Po przemnożeniu równania 2.3 przez <math>(X^{\pm}_{CSP} (t))^T</math> otrzymamy macierz kowariancji przetransformowanych sygnałów uśrednioną po realizacjach:
:<math>R^{\pm}_{CSP} (t) = X^{\pm}_{CSP} (t)(X^{\pm}_{CSP} (t))^T = P^T X^{\pm} (t) (X^{\pm}(t))^T P = P^T R^{\pm}P</math>  (2.4)
Gdzie <math>R^{\pm}</math> to macierz kowariancji sygnału uśredniona po realizacjach. Z warunków 1 i 2 wynika, że macierze <math>R^{+}_{CSP}</math> i <math>R^{-}_{CSP}</math> muszą być diagonalne, natomiast z warunku 3, że ich suma daje macierz jednostkową:
:<math>R^{+}_{CSP} + R^{-}_{CSP} = 1 </math>  (2.5)
Tak więc suma par diagonalnych wartości (<math>k^{+}_{i}</math> i <math>k^{-}_{i}</math>) w macierzach <math>R^{+}_{CSP}</math> i <math>R^{-}_{CSP}</math> musi być równa 1.
Korzystając z równania 2.5 wartości na diagonali można również zapisać za pomocą wzoru:

:<math>k^{+}_{i} = \vec{p}^{T}_{i} R^{+} \vec{p}_{i}    </math> (2.6)
:<math>k^{-}_{i} = \vec{p}^T_{i} R^{-}\vec{p}_{i} </math> (2.7)
gdzie <math>\vec{p}_{i}</math> to kolumnowy wektor macierzy <math>P</math>. Po przekształceniach ilorazu równań 2.6 i 2.7 można
otrzymać równanie:

:<math>R^{+} \vec{p}_{i} = \frac{k^{+}_{i}}{k^{-}_{i}} R^{-} \vec{p}_{i} </math>  (2.8)

Równanie to przedstawia ogólną formę zagadnienia wartości własnych. Takie przedstawienie problemu umożliwia zastosowanie do jego rozwiązania wydajnych metod algebraicznych.


Wektor własny <math>\vec{p}_i</math> jest interpretowany jako filtr przestrzenny. Dzięki transformacie <math>P</math> sygnał zostaje przeniesiony do przestrzeni, dla której różnica wariancji dla poszczególnych klas jest największa. Zgodnie z równaniem 2.5 najbardziej różniące się od siebie kanały są skorelowane z największą wartością własną <math>k^{+}</math>.
-->

===Ćwiczenie symulacyjne ===
<source lang  = matlab>
% symulowany eksperyment składa się z sinusoidy udającej alfę spoczynkową i
% gausa udającego potencjał wywołany
% źródła te są symulowane niezależnie a potem mieszane przez macierz L
% symulujemy źródła
% s1 - symuluje alfę
% s2 - symuluje "potencjał wywołany" (ERP)

%ustawiamy parametry do symulacji sygnałów
Fs = 100;
T = 1;
t = 0:1/Fs:T-1/Fs;
N_rep = 100;
N_chan = 2;
s = zeros(N_rep,N_chan, length(t));
X = zeros(N_rep,N_chan, length(t));

% filtr przestrzenny - z takimi wagami trzeba wziąść kanały EEG aby odzyskać sygnały źródłowe
P = [1 2
    1.5 1.3]; 


% topografie - z takimi wagami źródła dokładają się do poszczególnych elektrod
A = P^(-1); 


for r =1:N_rep % tworzymy kolejne realizacje "eksperymentu"
    s1 = sin(2*pi*11*t +pi/2+ 0*2*pi*rand())+ 0.02*randn(size(t));  % źródło alfa
    s2 = exp(-((t-0.8)/0.05).^2)+ 0.01*randn(size(t));              % źródło ERP
       
    s(r,1,:) = s1;
    s(r,2,:) = s2;
    tmp = squeeze(s(r,:,:));
    n = 0*randn(size(tmp));
    X(r,:,:) = A*tmp +n; % rzutujemy sygnały źródłowe na elektrody s -> x
    
end

% wycinamy warunki 
% baseline_ind   to indeksy pierwszej połowy każdego powtórzenia "baseline"
% ERP_ind        to indeksy drugiej połowy każdego powtórzenia zawierająca "ERP"
baseline_ind = find(t<0.5);
ERP_ind = find(t>=0.5);

x_baseline_kan_1 = X(:,1,baseline_ind);
x_baseline_kan_2 = X(:,2,baseline_ind);

x_ERP_kan_1 = X(:,1,ERP_ind);
x_ERP_kan_2 = X(:,2,ERP_ind);



% liczymy średnie macierze kowariancji:
R_E = zeros(N_chan,N_chan);
R_B = zeros(N_chan,N_chan);
for r =1:N_rep
    B = squeeze(X(r,:,baseline_ind)); 
    R_B = R_B + B*B' ;
    
    E = squeeze(X(r,:,ERP_ind));
    R_E = R_E + E*E' ;
end

R_B = R_B/N_rep;
R_E = R_E/N_rep;

% rozwiązujemy uogólnione zagadnienie własne
[W,Lambda]=eig(R_E,R_B); % możliwa jest też otymalizacja wzg. średniej macierzy kowariancji (R_B+R_A)/2);

% odzyskujemy sygnały źródeł
for r =1:N_rep
    S(r,:,:) = W'*squeeze(X(r,:,:));
end

% pobieramy wycinki odpowiadające obu częściom eksperymentu z estymowanych
% źródeł
s_baseline_estymowany_kan1 = squeeze(  S(:,1,baseline_ind));
s_baseline_estymowany_kan2 = squeeze(  S(:,2,baseline_ind));

s_ERP_estymowany_kan1 = squeeze(S(:,1,ERP_ind));
s_ERP_estymowany_kan2 = squeeze(S(:,2,ERP_ind));

%%%%%%%%%%%%%% Iloustracje %%%%%%%%%%%%%%%%%%%%%%%
% ilustracja sygnałów mierzonych
figure(1);clf
subplot(2,2,1);
    plot(t(baseline_ind),(squeeze(X(:,1,baseline_ind)))','b'); hold on
    plot(t(ERP_ind),(squeeze(  X(:,1,ERP_ind)))','r'); hold off
    xlabel('elektroda 1')
    title('ilustracja sutuacji pomiarowej -\newline znane są potencjały na elektrodach w dwóch warunkach eksperymentalnych')
subplot(2,2,3);
    plot(t(baseline_ind),(squeeze(X(:,1,baseline_ind)))','b'); hold on
    plot(t(ERP_ind),(squeeze(  X(:,2,ERP_ind)))','r'); hold off
    xlabel('elektroda 2')
subplot(1,2,2)
    plot(x_baseline_kan_1(:),x_baseline_kan_2(:),'b.');
    hold on
    plot(x_ERP_kan_1(:),x_ERP_kan_2(:),'r.');
    xlim([-2,2])
    ylim([-2,2])
    axis equal
    
    % wektor własny odpowiadajacy największej wartości własnej jest
    % kierunkiem najbardziej różnicującym warunki eksperymentalne
    disp('wartości własne znajdują się na diagonali macierzy Lambda')
    disp(Lambda)
    % rysujemy wersory jednostkowe w kierunkach wektorów własnych
    w1 = W(:,1); 
    w1 = w1/norm(w1);
    w2 = W(:,2); 
    w2 = w2/norm(w2);
    line([0, w1(1) ],[0,w1(2)],'Color',[0,0.3,0])
    text(w1(1),w1(2),'wektor własny1')
    line([0, w2(1) ],[0,w2(2)],'Color',[1,0,1])
    text(w2(1),w2(2),'wektor własny2')
    
    
    xlabel('Amplituda na elektrodzie 1')
    ylabel('Amplituda na elektrodzie 2')   
    legend('baseline','ERP')

% Ilustracja estymowanych źródeł
figure(2);clf
subplot(2,2,1);
    plot(t(baseline_ind),(squeeze(S(:,1,baseline_ind)))','b');hold on
    plot(t(ERP_ind),(squeeze(S(:,1,ERP_ind)))','r');hold off
    xlabel('estymowane zrodlo  1')
    title('ilustracja estymacji -\newline estymowane są potencjały żródeł w dwóch warunkach eksperymentalnych')    
subplot(2,2,3);
    plot(t(baseline_ind),(squeeze(S(:,2,baseline_ind)))','b');hold on
    plot(t(ERP_ind),(squeeze(S(:,2,ERP_ind)))','r');hold off
    xlabel('estymowane zrodlo  2');
subplot(1,2,2)
    plot(s_baseline_estymowany_kan1(:),s_baseline_estymowany_kan2(:),'b.');
    hold on
    plot(s_ERP_estymowany_kan1(:),s_ERP_estymowany_kan2(:),'r.');   
    
    xlabel('Amplituda estym. źródła 1')
    ylabel('Amplituda estym. źródła 2')   
    legend('baseline','ERP')
</source>

==Zastosowanie filtra CSP do detekcji potencjału P300==
===Eksperyment===
====Przygotowanie do badania:====
* założyć czepek z elektrodami w systemie 10-20
* elektrody referencyjne: M1 i M2
* elektroda GND w pozycji AFz

====Przygotowanie scenariuszy obci ====
* w terminalu uruchomić <tt>obci srv</tt>
* w terminalu uruchomić <tt>obci_gui --preset brain2013</tt>
* w interfejsie GUI zapisujemy scenariusze do własnego katalogu np "P300"
**"P-Brain2013 Signal (with ID)" jako np.  "Sygnal"  
** "P-Brain2013 Calibration p300 "  jako "kalibracjaP300" 
** "P-Brain 2013 p300" jako Labirynt
* w przeglądarce plików otwórz katalog ~/.obci/scenarios/P300. Powinien on zawierać pliki:  Sygnal.ini, kalibracjaP300.ini, Labirynt.ini oraz katalogi  Sygnal_configs, kalibracjaP300_configs, Labirynt_configs.
*  edytujemy parametry peerów.
** z katalogu "~/.obci/scenarios/P300/Sygnal_configs" kopiujemy plik "amplifier.ini" do katalogu "~/.obci/scenarios/P300" jako "global_amplifier.ini" . To pozwoli nam zmieniać ustawienia wzmacniacza dla wszystkich scenariuszy jednocześnie.
** w naszych scenariuszach zapisanych w plikach Sygnal.ini, kalibracjaP300.ini, Labirynt.ini podmieniamy ścieżkę <tt>config =</tt> w sekcji <tt>[peers.amplifier]</tt> tak aby 
pokazywała na ten skopiowany plik "global_amplifier.ini"
** w tym pliku "global_amplifier.ini" podmieniamy linijki (tak aby były to listy faktycznie wykorzystywanych kanałów) z:
*** active_channels
*** channel_names
** dodajemy też linijkę: sampling_rate = 256
* wchodzimy po kolei do katalogów: Sygnal_configs, kalibracjaP300_configs, Labirynt_configs i odnajdujemy plik <tt>switch_backup.ini</tt>. W tym pliku ustawiamy parametr <tt>new_scenario</tt> na pusty. To spowoduje, że scenariusze te nie będą uruchamiać kolejnych scenariuszy po zakończeniu działania.
** edytujemy plik "~/.obci/scenarios/P300/kalibracjaP300_configs/clasifier.ini" 
*** zmieniamy linię
:: ignore_channels = DriverSaw;AmpSaw;PO7;PO8
:: na 
:: ignore_channels = DriverSaw;AmpSaw;A1;A2
::* oraz linię:
:: montage_channels = PO7;PO8
:: na 
:: montage_channels = A1;A2

====Przeprowadzenie badania:====
# Uruchom scenariusz "Sygnał" 
# Tworzy on w katalogu domowym plik o nazwie <tt>file_id_name</tt>.
# Uruchamiamy svaroga z terminala poleceniem <tt>svarog</tt>. W zakładce sygnały on-line odnajdujemy nazwę naszego scenariusza "Sygnal". Podłączamy się do niego i  poprawiamy ewentualnie źle kontaktujące elektrody.
# Jak już jesteśmy zadowoleni z jakości sygnału to zatrzymujemy scenariusz "Sygnal" w obci.
# W pliku <tt>file_id_name</tt>  znajduje się string, który stanowi rdzeń do tworzenia nazw plików, z których korzystają nasze scenariusze. Proszę zmienić ten string np. na: "test1"
# Uruchamiamy scenariusz "kalibracjaP300". Badany będzie oglądał interfejs z trzema literami A B C migającymi w losowej kolejności. Zadaniem jest zliczanie mignięć litery B.
# Po zakończeniu kalibracji uruchamiamy scenariusz Labirynt.
# Danych z kalibracji potrzebować będziemy kilka zestawów.  Proszę powtórzyć kilku krotnie scenariusz kalibracjaP300. Przed każdym uruchomieniem trzeba zmienić string w pliku <tt>file_id_name</tt> np na "test???" gdzie ??? oznacza kolejne numery.

===Analiza wstępna===
* Wczytać dane kalibracyjne do Matlaba i pociąć je na realizacje typu T- target (związane z wystąpieniami litery "B") i NT- non-target (pozostałe litery) o długości -200 do +800 ms wokół trigerów. Dla każdej realizacji odjąć trend liniowy.
* Sygnał zmontować wzg. "średnich uszu" i wyświetlić średnie przebiegi dla warunku T i NT w układzie topograficznym - wykorzystać w tym celu funkcję <tt>topoplot</tt> z pakietu <tt>eeglab</tt>

Poniżej zaprezentowany jest przykładowy skrypt do cięcia danych wokół znaczników. Działa on z plikami zawartymi w archiwum:
: [[Plik:KalibracjaP300.tar.gz]] 
Korzysta z funkcji pomocniczych dostępnych w dystrybucji obci w katalogu 
: /usr/share/openbci/analysis/matlab_obci_signal_processing
Openbci można pobrać z https://github.com/BrainTech/openbci

<source lang = matlab>
% ustalamy nzawy plików z danymi
nazwaPliku = 'p_6301423_calibration_p300.obci';
nameOfXMLFile = strcat(nazwaPliku,'.xml');
nameOfTagFile = strcat(nazwaPliku,'.tag'); %tagi = znaczniki zdarzeń
namesOfDataFiles = strcat(nazwaPliku,'.raw');

% inicjujemy obiekt rm
rm = ReadManager(nameOfXMLFile,namesOfDataFiles,nameOfTagFile);

% obieramy przydatne parametry i znaczniki
numberOfChannels  = rm.get_param('number_of_channels');
namesOfChannels   = rm.get_param('channels_names');
samplingFrequency = rm.get_param('sampling_frequency');
tagsStruct        = rm.get_tags();

% tworzenie list znaczników Target i NonTarget
numberOfStruct = length(tagsStruct);
targetTimeStamps = [];
NonTargetTimeStamps = [];
for structNumber = 1:numberOfStruct % iterujemy się przez tagi
    if(strcmp(tagsStruct(structNumber).name,'blink')) % szukamy tagów o nazwie 'blink'
        index = tagsStruct(structNumber).children.index; % tu jest numer pola stymulacji
        target= tagsStruct(structNumber).children.target;% tu jest numer pola na którym wyświetlany jest target
        if index == target % warunek na to, że mamy do czynienia z tagiem target
            targetTimeStamps = [targetTimeStamps tagsStruct(structNumber).start_timestamp]; %dodajemy timeStamp do listy targetów
        else
            NonTargetTimeStamps = [NonTargetTimeStamps tagsStruct(structNumber).start_timestamp];%dodajemy timeStamp do listy non-targetów
        end
        
    end
end



% pobieramy próbki
samples = double(rm.get_samples()); % konwersja na double jest potrzebna żeby dobrze funkcjonowało filtrowanie
samples=samples(1:8,:); % odrzucamy kanały, które nie mają EEG
numberOfChannels =8;

% filtrujemy dolnoprzepustowo aby odrzucić artefakty sieci i część
% artefaktów mięśniowych
[b,a] = cheby2(6,80,25 /(samplingFrequency/2),'low');
for ch = 1:numberOfChannels
    samples(ch,:)=filtfilt(b,a,samples(ch,:));
end

% montujemy dane do wspólnej średniej (common average)
M = -ones(8,8)/8;
M=M+eye(8,8)*9/8;
samples = 0.0715*M*samples;

% wycinamy dane wokół znaczników
PRE = -0.2; % czas przed tagiem w sek.
POST = 0.8; % czas po tagu w sek.
wycinek = floor(PRE*samplingFrequency:POST*samplingFrequency); % tablica ze "standardowymi' indeksami do cięcia

% pobieramy targety
TargetSignal = zeros(length(targetTimeStamps),numberOfChannels, length(wycinek)); % tablica na sygnały target
for trialNumber = 1:length(targetTimeStamps)
    trigerOnset = floor(targetTimeStamps(trialNumber)*samplingFrequency);
    tenWycinek = wycinek + trigerOnset;
    if tenWycinek(1)>0 && tenWycinek(end)<=size(samples,2) % test czy wycinek który chcemy pobrać nie wystaje poza dostępny sygnał
        tmpSignal = samples(:,tenWycinek);
        tmpSignal = detrend(tmpSignal')'; % usuwanie liniowego trendu - przy krótkich wycinkach działa lepiej niż filtrowanie górnoprzepustowe
        TargetSignal(trialNumber, :,:) = tmpSignal;
    end
end

% pobieramy non-targety
NonTargetSignal = zeros(length(NonTargetTimeStamps),numberOfChannels, length(wycinek));% tablica na sygnały non-target
for trialNumber = 1:length(NonTargetTimeStamps)
    trigerOnset = floor(NonTargetTimeStamps(trialNumber)*samplingFrequency);
    tenWycinek = wycinek + trigerOnset;
    if tenWycinek(1)>0 && tenWycinek(end)<=size(samples,2)
        tmpSignal = samples(:,tenWycinek);
        tmpSignal = detrend(tmpSignal')';
        NonTargetSignal(trialNumber, :,:) = tmpSignal;
    end
end
%
% dla ilustracji podglądamy średnie po powtórzeniach ze wszystkich targetó∑
% i non-targetów
plot(squeeze(mean(TargetSignal,1))','r');
hold on
plot(squeeze(mean(NonTargetSignal,1))','b')
</source>

===Analiza CSP===
* Korzystając z danych kalibracyjnych wykonac analizę CSP wzmacniającą potencjał P300.
* Zaprezentować średnią ze wszystkich kanałów źródłowych z warunku target (jeden kolor) i non-target (inny kolor) w subplotach ułożonych w prostokątnej siatce. Zaobserwować dla którego kanału średnie różnią się najbardziej. Czy jest związek tego kanału z wartościami własnymi?

* Dla kanału najbardziej różnicującego wykonać mapki topograficzne wektorów odpowiadających:
** filtrowi przestrzennemu 
** rzutu topograficznego źródła na elektrody.
* Do wykonania tych mapek wykorzystać funkcję <tt>topoplot</tt> z pakietu <tt>eeglab</tt>
* Zbadać powtarzalność topografii pomiędzy plkiami konfiguracyjnymi.

===Wybór i separacja cech===
* przedstaw na rysunkach nałożone na siebie pojedyncze realizacje target i non-target po rzutowaniu na wetor <math>w</math> odpowiadający największej i kolejnej wartości  własnej. 
* przedstaw wykresy punktowe, takie że na jednej osi jest moc sygnału odpowiadającego największej wartości włanej a na drugiej osi kolejnej wartości własnej; jeden punkt reprezentuje jedno powtórzenie.
* wykonaj serię wykresów jak w poprzednim punkcie dla uśrednieni kolejno po 2,4,6,8,10 realizacjach. Zaobserwuj jak zmienia się separacja targetów i non-targetów.

==Filtry przestrzenne dla większej ilości warunków==
===FFDIAG===
===Analiza ERD/S z użyciem FFDIAG===

==Filtry przestrzenne dla SSEP ==

=== Teoria===
Proszę zapoznać się z koncepcją filtra przestrzennego dla SSVEP zaprezentowaną tu: http://www.eurasip.org/Proceedings/Eusipco/Eusipco2009/contents/papers/1569193209.pdf

===Eksperyment ASSR===
W eksprymencie wykorzystujemy układ do generacji potencjałów słuchowych stanu ustalonego (ASSR). Wejście układu ASSR typu mini-jack wkładamy w wyjście słuchawkowe w laptopie. Drugie wejście układu ASSR wkładamy do wyjścia triggera we wzmacniaczu. Uruchamiamy plik dźwiękowy MM40tr.wav. Można go znalezc w: http://www.fuw.edu.pl/~suffa/LabEEG/MM40tr.wav

Stymulacja dźwiękowa składa sie z fali nośnej o częstości 400 Hz modulowanej z częstością 40 Hz. Plik dźwiękowy zawiera 5 sekund ciszy i 5 sekund stymulacji, powtórzone 40 razy.

====Rejestracja sygnału====
# Zakładamy czepek i elektrody w systemie 10-10, dbamy o to by opory pomiędzy elektrodami były poniżej 5 k G i różnice pomiędzy oporami różnych elektrod nie przekraczały 20%.
# Oklejamy kwadrat 3x3 elektrod na korze słuchowej z lewej strony (elektrody FT7, FC5, FC3, T7, C5, T3, TP7, CP5, CP3), 3x3 elektrod na korze słuchowej z prawej strony (elektrody FT8, FC6, FC4, T8, C6, T4, TP8, CP6, CP4), elektrody Fz, Cz, Pz i Oz, elektrody referencyjne A1 i A2. W sumie powinno byc 24 elektrody.
# Elektrodę GND mocujemy na pozycji AFz.
# Sygnał rejestrujemy z częstością 2048 Hz.
# Do rejestracji stosujemy scenariusz 'ASSR' w interfejsie obci_gui

====Analiza====
 JZ: zmieniłbym analizę na czas-częstość i zrobił porównanie montażu usznego do filtra G.G. Moliny

Początek stymulacji dźwiękowej oznaczymy jako 0. Poniższą analizę zastosuj dla sygnałów w referencji do uśrednionych odprowadzeń usznych A1 i A2.
Wyznaczenie pasma częstości odpowiedzi ASSR
# Z sygnału wycinamy fragmenty od 0 do 5 sek. dla wszystkich elektrod położone nad korą słuchową.
# Dla każdej realizacji obliczamy widma metodą Welcha.
# Otrzymane zespolone widma uśredniamy po realizacjach.
# Sprawdzamy czy w uśrednionym widmie występuję maksimum w częstości modulacji tj. 40 Hz.

====Wyznaczenie przebiegu czasowego ERD i ERS====
# Zaprojektuj filtry pasmowo przepustowe (Chebyszewa 2 rodzaju) zgodne z wyznaczonym pasmem. Zbadaj funkcje przenoszenia i odpowiedzi impulsowej.
# Powycinaj sygnały od -5 do +10 sekund (wszystkie kanały). Przefiltruj każdą realizację.
# Oblicz moc chwilową za pomocą transformaty Hilberta (kwadrat modułu transformaty Hilberta).
# Uśrednij moc chwilową po realizacjach.
# Oblicz względną zmianę mocy chwilowej względem czasu -4 do -2. W ten sposób otrzymasz przebieg ERD i ERS w czasie.
# Wykreśl ERD i ERS w układzie topograficznym. (Rozmieść subploty tak, aby z w przybliżeniu odpowiadały pozycjom elektrod).
====Transformacja Hjortha====
Transformacja Hjortha jest przybliżeniem numerycznym transformacji Laplace'a, czyli drugiej pochodnej przestrzennej. Obliczamy ją jako różnicę potencjału pomiędzy daną elektrodą i średnią z czterech sąsiednich elektrod.
Przelicz potencjały z elektrod, w których występuję odpowiedź ASSR na montaż Hjortha i powtórz analizę ERD/ERS opisaną powyżej.

==ICA jako filtr przestrzenny==