Aktualna wersja na dzień 11:05, 14 kwi 2026

Spis treści

1 Plan zajęć
2 Sesja 1: Ślepa separacja źródeł i CSP — wprowadzenie
3 Sesja 2: CSP dla P300 — dane realne
4 Sesja 3: CSP dla P300 — interpretacja i separacja cech
5 Sesja 4: Filtry przestrzenne dla SSEP — cosSinCSP
- 5.1 Filtry przestrzenne dla SSEP
6 Sesja 5: ICA — teoria i wydobywanie interesującego komponentu
- 6.1 ICA jako filtr przestrzenny
7 Sesja 6: ICA — artefakty i synteza metod
- 7.1 ZADANIE: Identyfikacja artefaktów
  - 7.1.1 W raporcie:

Plan zajęć

Materiał realizowany jest w ciągu 3 tygodni (6 spotkań po ~2,5 h).

Sesja	Temat	Kluczowe sekcje na tej stronie
1	Wykład BSS + CSP; ćwiczenie symulacyjne	Ślepa separacja źródeł, Common Spatial Pattern, Ćwiczenie symulacyjne
2	Dane P300: wczytanie, cięcie epok, wizualizacja ERP, implementacja CSP	Analiza wstępna, Analiza CSP
3	Mapki topograficzne; separacja cech	Analiza CSP (cd.), Wybór i separacja cech
4	cosSinCSP jako uogólnienie CSP; dane SSVEP	Filtry przestrzenne dla SSVEP
5	Wykład ICA; komponenty alfa	ICA jako filtr przestrzenny, Komponenty alfa
6	Artefakty (ICLabel/MARA); synteza CSP/cosSinCSP/ICA	Identyfikacja artefaktów

Sesja 1: Ślepa separacja źródeł i CSP — wprowadzenie

Czas: ~60 min wykład + ~90 min ćwiczenie symulacyjne

Zajęcia wprowadzają do problemu filtracji przestrzennej sygnałów EEG. Punktem wyjścia jest ogólny problem ślepej separacji źródeł (BSS), z którego CSP wyłania się jako szczególny przypadek dla dwóch warunków eksperymentalnych. Filtry przestrzenne omawiane na tych zajęciach stanowią fundament klasycznych metod BCI (P300, SSVEP, motor imagery) oraz są bezpośrednim odpowiednikiem warstw przestrzennych w sieciach neuronowych takich jak ShallowConvNet i EEGNet — do tej analogii wrócimy w sesji 6.

Slajdy do wykładu (PDF)

Zakres wykładu

model generatywny EEG: [math]x(t) = As(t)[/math], macierz kowariancji, idea diagonalizacji
BSS jako ogólny problem separacji źródeł
CSP jako szczególny przypadek: dwa warunki eksperymentalne, iloraz Rayleigha, uogólnione zagadnienie własne
interpretacja filtrów i topografii źródeł

Ćwiczenie symulacyjne

Zob. Ćwiczenie symulacyjne w sekcji CSP poniżej.

Ćwiczenie symulacyjne

Zadania do samodzielnej realizacji

Poniższy kod dostarcza gotowej infrastruktury: generacji sygnałów, macierzy mieszającej i wizualizacji. Zadaniem jest samodzielna implementacja kluczowych kroków analizy CSP.

Zadanie 1 (obowiązkowe): Obliczenie macierzy kowariancji

Napisz funkcję srednia_kowariancja(X, ind), która dla danej tablicy sygnałów X (wymiary: powtórzenia × kanały × próbki) i wektora indeksów czasowych ind zwraca średnią macierz kowariancji uśrednioną po powtórzeniach. Dla każdego powtórzenia zastosuj funkcję cov i znormalizuj wynik przez jego ślad (trace). Sprawdź wynik porównując z macierzami R_B i R_E obliczonymi w dostarczonym kodzie.

Zadanie 2 (obowiązkowe): Implementacja CSP

Korzystając z funkcji eig oraz macierzy kowariancji obliczonych w zadaniu 1, wyznacz macierz filtrów przestrzennych W i odpowiadające im wartości własne. Odtwórz sygnały źródłowe dla wszystkich powtórzeń. Narysuj chmury punktów (amplituda źródła 1 vs amplituda źródła 2) osobno dla warunków baseline i ERP — przed transformacją CSP i po niej. Co zmieniło się w kształcie chmur?

Wykreśl także wizualizację przebiegów czasowych sygnałów X i sygnałów S uzyskanych po rozmieszaniu.

Zadanie 3 (obowiązkowe): Interpretacja wartości własnych

Wartości własne znajdują się na przekątnej macierzy Lambda. Który filtr najbardziej różnicuje warunki baseline i ERP? Jaki jest związek między wartością własną a separacją widoczną na wykresach punktowych?

Zadanie 4 (dla chętnych): Wpływ macierzy mieszającej

W dostarczonym kodzie macierz P (filtr przestrzenny) wynosi:

P = [1 2; 1.5 1.3]

Zmień ją kolejno na macierz bliską jednostkowej (np. P = [1 0.1; 0.1 1]) oraz na macierz z dużymi elementami pozadiagonalnymi (np. P = [1 5; 4 1]). Jak zmienia się kształt chmury punktów przed transformacją CSP? Czy CSP poprawnie oddziela źródła w obu przypadkach?

Zadanie 5 (dla chętnych): Wpływ szumu

W kodzie szum jest wyłączony: n = 0*randn(size(tmp)). Usuń mnożnik 0* i zbadaj jak poziom szumu wpływa na separację źródeł. Przy jakim stosunku sygnału do szumu CSP przestaje skutecznie rozdzielać warunki?

Pytanie do dyskusji

Dostarczony kod używa eig(R_E, R_B). Alternatywą jest eig(R_E, (R_E+R_B)/2). Czym różnią się te dwie wersje? W jakich sytuacjach eksperymentalnych druga wersja mogłaby być bardziej uzasadniona?

Szkielet kodu do uzupełnienia

Poniższy szkielet zawiera gotową infrastrukturę symulacji. Uzupełnij brakujące fragmenty oznaczone komentarzem % TUTAJ zgodnie z zadaniami 1–3. Pełny kod referencyjny dostępny jest poniżej — zajrzyj do niego dopiero po samodzielnej próbie.

%% Parametry symulacji — nie modyfikuj tej sekcji
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));

P = [1 2; 1.5 1.3];   % filtr przestrzenny (prawdziwy)
A = P^(-1);            % topografia źródeł

for r = 1:N_rep
    s1 = sin(2*pi*11*t + pi/2) + 0.02*randn(size(t));
    s2 = exp(-((t-0.8)/0.05).^2) + 0.01*randn(size(t));
    s(r,1,:) = s1;
    s(r,2,:) = s2;
    tmp = squeeze(s(r,:,:));
    n = 0*randn(size(tmp));   % Zadanie 5: zmień 0* na 1* aby włączyć szum
    X(r,:,:) = A*tmp + n;
end

baseline_ind = find(t < 0.5);
ERP_ind      = find(t >= 0.5);

R_B = zeros(N_chan, N_chan);
R_E = zeros(N_chan, N_chan);

%% Zadanie 1: Oblicz średnie macierze kowariancji
% korzystając z funkcji srednia_kowariancja(X, ind)

% Dla każdego powtórzenia r:
%   - wytnij fragment sygnału X(r,:,baseline_ind) i analogicznie ERP
%   - oblicz macierz kowariancji funkcją cov() — pamiętaj o transpozycji
%   - znormalizuj przez ślad (funkcja trace())
%   - dodaj do sumy
% Na końcu podziel przez N_rep







%% Zadanie 2: Rozwiąż uogólnione zagadnienie własne
% Użyj funkcji eig(A,B) gdzie A=R_E, B=R_B
% W wyniku otrzymasz macierz wektorów własnych W (w kolumnach)
% i macierz Lambda z wartościami własnymi na przekątnej

% TUTAJ: [W, Lambda] = ...

disp('Wartości własne (przekątna Lambda):')
disp(diag(Lambda))

%% Zadanie 2 cd.: Odtwórz sygnały źródłowe
% Dla każdego powtórzenia r zastosuj filtr: S(r,:,:) = W' * X(r,:,:)

S = zeros(N_rep, N_chan, length(t));

for r = 1:N_rep
    % TUTAJ: S(r,:,:) = ...
end

%% Zadanie 3: Wizualizacja — chmury punktów przed i po CSP
x_base_k1 = X(:,1,baseline_ind);
x_base_k2 = X(:,2,baseline_ind);
x_ERP_k1  = X(:,1,ERP_ind);
x_ERP_k2  = X(:,2,ERP_ind);

s_base_k1 = squeeze(S(:,1,baseline_ind));
s_base_k2 = squeeze(S(:,2,baseline_ind));
s_ERP_k1  = squeeze(S(:,1,ERP_ind));
s_ERP_k2  = squeeze(S(:,2,ERP_ind));

figure(1); clf
subplot(1,2,1)
    plot(x_base_k1(:), x_base_k2(:), 'b.'); hold on
    plot(x_ERP_k1(:),  x_ERP_k2(:),  'r.');
    axis equal; xlim([-2 2]); ylim([-2 2])
    xlabel('Kanał 1'); ylabel('Kanał 2')
    title('Przed CSP')
    legend('baseline','ERP')

subplot(1,2,2)
    plot(s_base_k1(:), s_base_k2(:), 'b.'); hold on
    plot(s_ERP_k1(:),  s_ERP_k2(:),  'r.');
    axis equal
    xlabel('Źródło 1'); ylabel('Źródło 2')
    title('Po CSP')
    legend('baseline','ERP')
%% Wizualizacja przebiegów czasowych X i S
% TUTAJ:

%% Zadanie 4 (dla chętnych): zmień macierz P powyżej i powtórz analizę
%% Zadanie 5 (dla chętnych): włącz szum (zmień 0* na 1* w linii z n=...)

Rozwiązanie — pokaż dopiero po samodzielnej próbie

% symulowany eksperyment składa się z sinusoidy udającej alfę spoczynkową i
% funkcji Gaussa udającego potencjał wywołany
% źródła te są symulowane niezależnie a potem mieszane przez macierz A
% 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)); 
    tmp =cov(B');
    R_B = R_B + tmp./trace(tmp);%B*B' ;
 
    E = squeeze(X(r,:,ERP_ind));
    tmp = cov(E');
    R_E = R_E + tmp./trace(tmp);%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ż optymalizacja 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));

%%%%%%%%%%%%%% Ilustracje %%%%%%%%%%%%%%%%%%%%%%%
% 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 sytuacji pomiarowej -\newline znane są potencjały na elektrodach w dwóch warunkach eksperymentalnych')
subplot(2,2,3);
    plot(t(baseline_ind),(squeeze(X(:,2,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 odpowiadający największej wartości własnej jest
    % kierunkiem najbardziej różnicującym warunki eksperymentalne
    disp('wartości własne znajdują się na przekątnej 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łasny 1')
    line([0, w2(1) ],[0,w2(2)],'Color',[1,0,1])
    text(w2(1),w2(2),'wektor własny 2')
    
    
    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')

Sesja 2: CSP dla P300 — dane realne

Czas: ~150 min (jedno spotkanie)

Pracujemy z danymi kalibracyjnymi z eksperymentu P300. Celem sesji jest przejście pełnego pipeline'u: od surowych danych do odtworzonych źródeł CSP.

Zakres

wczytanie danych kalibracyjnych, identyfikacja znaczników T i NT
cięcie epok (−200 do +800 ms), usuwanie trendu liniowego
montaż względem połączonych uszu
wizualizacja ERP w układzie topograficznym
obliczenie macierzy kowariancji R_T i R_NT, normalizacja przez ślad
rozwiązanie uogólnionego zagadnienia własnego: [W, Lambda] = eig(R_T, R_NT)
odtworzenie sygnałów źródłowych: S = W'*EEG
przegląd estymowanych źródeł w dwóch warunkach

Materiały

Zob. Analiza wstępna i Analiza CSP poniżej.

Sesja 3: CSP dla P300 — interpretacja i separacja cech

Czas: ~150 min (jedno spotkanie)

Kontynuacja analizy CSP. Celem sesji jest interpretacja wyników w kategoriach fizjologicznych oraz ocena separowalności klas w przestrzeni cech opartej na mocy sygnału.

Zakres

mapki topograficzne filtra przestrzennego i topografii źródła (funkcja topoplot z pakietu EEGLAB)
związek wartości własnych z siłą różnicowania warunków T i NT
wykresy punktowe mocy: oś X — moc składowej z największą wartością własną, oś Y — moc kolejnej składowej; jeden punkt = jedno powtórzenie
uśrednianie po 2, 4, 6, 8 i 10 realizacjach — obserwacja poprawy separacji wraz z liczbą uśrednień
dyskusja: utożsamianie komponentów CSP z fizjologicznymi źródłami jest ryzykowne; komponenty maksymalizują kontrast między warunkami, niekoniecznie odpowiadają izolowanym generatorom

Materiały

Zob. Analiza CSP i Wybór i separacja cech poniżej.

Analiza wstępna

%% KROK 1: Wczytanie danych
% ReadManager to pomocnicza klasa do odczytu plików w formacie OpenBCI.
% Przyjmuje trzy pliki: XML z metadanymi, RAW z próbkami i TAG ze znacznikami zdarzeń.
nazwaPliku = 'p_6301423_calibration_p300.obci';
nameOfXMLFile  = strcat(nazwaPliku, '.xml');
nameOfTagFile  = strcat(nazwaPliku, '.tag');
namesOfDataFiles = strcat(nazwaPliku, '.raw');

rm = ReadManager(nameOfXMLFile, namesOfDataFiles, nameOfTagFile);

% Pobieramy podstawowe parametry nagrania
numberOfChannels  = rm.get_param('number_of_channels');
namesOfChannels   = rm.get_param('channels_names');
samplingFrequency = rm.get_param('sampling_frequency');

% Pobieramy listę wszystkich znaczników zdarzeń z nagrania
tagsStruct = rm.get_tags();

%% KROK 2: Identyfikacja znaczników Target i Non-Target
% Każdy znacznik typu 'blink' odpowiada jednemu mignięciu w eksperymencie P300.
% Pole 'index' mówi który symbol mignął, pole 'target' który symbol był oczekiwany.
% Jeśli index == target to było to mignięcie docelowe (Target), w przeciwnym razie
% Non-Target.
targetTimeStamps    = [];
NonTargetTimeStamps = [];

for structNumber = 1:length(tagsStruct)
    if strcmp(tagsStruct(structNumber).name, 'blink')
        index  = tagsStruct(structNumber).children.index;
        target = tagsStruct(structNumber).children.target;
        if index == target
            targetTimeStamps    = [targetTimeStamps    tagsStruct(structNumber).start_timestamp];
        else
            NonTargetTimeStamps = [NonTargetTimeStamps tagsStruct(structNumber).start_timestamp];
        end
    end
end

%% KROK 3: Wczytanie i filtrowanie sygnału
samples = double(rm.get_samples());
samples = samples(1:8,:);   % odrzucamy kanały bez EEG
numberOfChannels = 8;

% Filtr dolnoprzepustowy Czebyszewa 2. rodzaju — odcięcie 25 Hz
[b, a] = cheby2(6, 80, 25/(samplingFrequency/2), 'low');
for ch = 1:numberOfChannels
    samples(ch,:) = filtfilt(b, a, samples(ch,:));
end

%% KROK 4: Montaż względem średniej (common average reference)
% Odejmujemy od każdego kanału średnią ze wszystkich kanałów w danej chwili.
% Jest to prostsze i bardziej ogólne niż montaż względem połączonych uszu.
M = -ones(numberOfChannels, numberOfChannels) / numberOfChannels;
M = M + eye(numberOfChannels) * (numberOfChannels+1) / numberOfChannels;
samples = 0.0715 * M * samples;

%% KROK 5: Cięcie epok
PRE  = -0.2;
POST =  0.8;
wycinek = floor(PRE*samplingFrequency : POST*samplingFrequency);

TargetSignal    = zeros(length(targetTimeStamps),    numberOfChannels, length(wycinek));
NonTargetSignal = zeros(length(NonTargetTimeStamps), numberOfChannels, length(wycinek));

for trialNumber = 1:length(targetTimeStamps)
    trigerOnset = floor(targetTimeStamps(trialNumber) * samplingFrequency);
    tenWycinek  = wycinek + trigerOnset;
    if tenWycinek(1) > 0 && tenWycinek(end) <= size(samples,2)
        tmpSignal = samples(:, tenWycinek);
        tmpSignal = detrend(tmpSignal')';   % usuwanie trendu liniowego
        TargetSignal(trialNumber,:,:) = tmpSignal;
    end
end

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

% Oś czasu w milisekundach — przyda się do rysowania
times = (wycinek / samplingFrequency) * 1000;

%% KROK 6: Wizualizacja ERP w układzie topograficznym
% Oblicz potencjały wywołane (ERP) jako średnie po powtórzeniach
ERP_T  = squeeze(mean(TargetSignal,    1));   % wymiary: kanały × czas
ERP_NT = squeeze(mean(NonTargetSignal, 1));

% TUTAJ: narysuj ERP_T i ERP_NT dla każdego kanału w subplotach
% ułożonych w siatce odpowiadającej przybliżonym pozycjom elektrod.
% Na każdym subplocie nałóż oba warunki różnymi kolorami.
% Oś X: czas w ms (użyj wektora times), oś Y: amplituda w µV.
% Zadbaj o legendę i tytuły z nazwami kanałów (namesOfChannels).

%% KROK 7: Obliczenie średnich macierzy kowariancji
% Dla każdego powtórzenia oblicz macierz kowariancji funkcją cov(),
% znormalizuj przez jej ślad (trace()), a następnie uśrednij po powtórzeniach.
% Pamiętaj: cov() oczekuje macierzy próbki × kanały, czyli potrzebujesz transpozycji.

R_T  = zeros(numberOfChannels, numberOfChannels);
R_NT = zeros(numberOfChannels, numberOfChannels);

for r = 1:length(targetTimeStamps)
    % TUTAJ: wytnij r-te powtórzenie T, oblicz kowariancję, znormalizuj, dodaj do R_T
end

for r = 1:length(NonTargetTimeStamps)
    % TUTAJ: wytnij r-te powtórzenie NT, oblicz kowariancję, znormalizuj, dodaj do R_NT
end

% TUTAJ: podziel R_T i R_NT przez odpowiednie liczby powtórzeń

%% KROK 8: Rozwiązanie uogólnionego zagadnienia własnego
% Szukamy filtrów przestrzennych W maksymalizujących stosunek mocy
% w warunku T względem NT: eig(R_T, R_NT)
% W kolumnach W znajdują się filtry, na przekątnej Lambda — wartości własne.

% TUTAJ: [W, Lambda] = ...

disp('Wartości własne (przekątna Lambda):')
disp(diag(Lambda))

%% KROK 9: Odtworzenie sygnałów źródłowych
% Dla każdego powtórzenia zastosuj filtr przestrzenny:
% s(r,:,:) = W' * x(r,:,:)
% Pamiętaj o wymiarach — squeeze() może być pomocne.

S_T  = zeros(size(TargetSignal));
S_NT = zeros(size(NonTargetSignal));

for r = 1:size(TargetSignal, 1)
    % TUTAJ: S_T(r,:,:) = ...
end

for r = 1:size(NonTargetSignal, 1)
    % TUTAJ: S_NT(r,:,:) = ...
end

%% KROK 10: Wizualizacja estymowanych źródeł
% Narysuj średnie przebiegi źródeł (ERP w przestrzeni CSP) dla obu warunków.
% Ułóż subploty w prostokątnej siatce — jeden subplot na każde źródło.
% Dla którego źródła różnica między T i NT jest największa?
% Czy odpowiada to źródłu z największą wartością własną?

ERP_S_T  = squeeze(mean(S_T,  1));
ERP_S_NT = squeeze(mean(S_NT, 1));

% TUTAJ: narysuj ERP_S_T i ERP_S_NT w subplotach

ZADANIE: Analiza CSP

Przegląd badań o P300: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2715154/

Link do Read menager [1]

Wykonać 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.

Wybór i separacja cech

Przedstaw na rysunkach nałożone na siebie pojedyncze realizacje z warunków target i non-target po rzutowaniu na wektor [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 (suma kwadratów wartości próbek w wybranym zakresie czasu) odpowiadającego największej wartości własnej, a na drugiej osi kolejnej (mniejszej) wartości własnej; jeden punkt reprezentuje jedno powtórzenie.
Wykonaj serię wykresów jak w poprzednim punkcie dla uśrednień sygnałów kolejno po 2, 4, 6, 8 i 10 realizacjach:
- Liczymy potencjał wywołany dla danej liczby powtórzeń.
- Następnie podnosimy wartości próbek do kwadratu
- i sumujemy je w wybranym zakresie czasu.
Zaobserwuj jak zmienia się separacja w grupach target i non-target.

Sesja 4: Filtry przestrzenne dla SSEP — cosSinCSP

Czas: ~150 min (jedno spotkanie)

Sesja wprowadza metodę cosSinCSP jako naturalne uogólnienie CSP o wiedzę dziedzinową dotyczącą sygnału SSEP. W klasycznym CSP dwa warunki eksperymentalne definiują dwie macierze kowariancji. W cosSinCSP rolę jednej z macierzy przejmuje projekcja sygnału na przestrzeń sinusów i cosinusów częstości stymulacji i jej harmonicznych. Formalizm ilorazu Rayleigha i uogólnionego zagadnienia własnego pozostaje identyczny.

Dane używane na tych zajęciach pochodzą z eksperymentu SSEP z archiwum — bez sesji rejestracyjnej. Na wykresach widma sygnału przed transformacją CSP dominuje artefakt sieciowy 50 Hz, który praktycznie zasłania odpowiedź na stymulację. Zadaniem jest sprawdzenie czy cosSinCSP potrafi ten sygnał wydobyć.

Zakres

idea cosSinCSP: macierz S złożona z sinusów i cosinusów harmonicznych, projekcja sygnału na przestrzeń SSEP i resztę Y
implementacja funkcji cosSinCSP w MATLAB
analiza danych archiwalnych: widma Welcha dla kanałów oryginalnych i źródeł CSP — porównanie SNR przed i po transformacji
prezentacja filtrów przestrzennych i topografii źródeł

Materiały

Zob. Filtry przestrzenne dla SSEP poniżej.

Filtry przestrzenne dla SSEP

Teoria

Ciekawa koncepcja filtra przestrzennego dla SSEP zaprezentowana jest tu: http://www.eurasip.org/Proceedings/Eusipco/Eusipco2009/contents/papers/1569193209.pdf

Pokrótce można ją rozumieć podobnie do tego co robiliśmy rozważając filtry przestrzenne CSP z tym, że dla SSEP oraz innych potencjałów wywołanych stanu ustalonego możemy skorzystać z dodatkowych informacji dotyczących poszukiwanych źródeł. Wiemy mianowicie, że powinny one oscylować z częstością bodźca i być może jej harmonicznych.

Przyda nam się macierz [math]S[/math] zbudowana tak, że w kolejnych kolumnach znajdują się sinusy i cosinusy kolejnych częstości harmonicznych. Wektory te unormujemy, żeby miały energię równą 1. Innymi słowy macierz [math]S[/math] zbudowana jest z wersorów rozpinających przestrzeń, w której powinien znajdować się sygnał SSEP.

W Matlabie możemy taką macierz zbudować tak:

% Fs                - częstość próbkowania
% numberOfSamples   - długość sygnału w próbkach
% numberOfHarmonics - liczba harmonicznych, które chcemy włączyć do analizy
t = (0:1:numberOfSamples - 1) / Fs;
S = zeros(numberOfSamples, 2*numberOfHarmonics);

for harmonicNumber = 1:numberOfHarmonics
    c = cos(2*pi*stimulationFrequency*harmonicNumber*t);
    s = sin(2*pi*stimulationFrequency*harmonicNumber*t);
    S(:,(harmonicNumber - 1)*2 + 1) = c/norm(c);
    S(:,(harmonicNumber - 1)*2 + 2) = s/norm(s);
end

Aby w badanym sygnale znaleźć składowe odpowiadające SSEP musimy rzutować sygnał [math]X[/math] (macierz sygnałów kanały × próbki) na przestrzeń rozpiętą przez [math]S[/math]:

[math]A = X S[/math]

Macierz [math]A[/math] zawiera współczynniki będące iloczynami skalarnymi sygnałów i wersorów — mówią one o tym jak dużo sinusa bądź cosinusa o danej częstości jest w pierwotnym sygnale. Komponenty SSEP zawarte w sygnale [math]X[/math] odzyskujemy tak:

[math]\mathrm{SSEP} = A S^T[/math]

Modelujemy rejestrowany sygnał jako:

[math]X = \mathrm{SSEP} + Y[/math]

gdzie:

[math]Y = X - \mathrm{SSEP}[/math]

to wszystkie komponenty sygnału, które nas nie interesują.

Filtr przestrzenny, który chcemy zbudować powinien maksymalizować stosunek wariancji [math]\mathrm{SSEP}[/math] do wariancji [math]Y[/math]. Macierz kowariancji powinna być uśredniona po powtórzeniach, a kowariancja sygnału w każdym powtórzeniu powinna być znormalizowana przez jej ślad. Dalej stosujemy technikę znaną z CSP: maksymalizację ilorazu Rayleigha przez rozwiązanie uogólnionego zagadnienia własnego dla macierzy kowariancji [math]\mathrm{SSEP}[/math] i [math]Y[/math].

Zadanie: implementacja funkcji cosSinCSP

Dane do zadania: Plik:PrzykladoweDaneSSVEP.mat.gz. Pełny przykład z danymi z eksperymentu SSEP: Plik:SSVEP demo csp.tar.gz

Zaimplementuj funkcję cosSinCSP zgodnie z poniższym szkieletem. Prawidłowo zaimplementowana funkcja wraz ze skryptem demonstracyjnym poniżej powinna generować rysunek:

function [W, Lambda] = cosSinCSP(signal, stimulationFrequency, numberOfHarmonics, Fs)
% cosSinCSP - filtr przestrzenny dla sygnałów SSEP
%
% Wejście:
%   signal                - dane EEG, wymiary: powtórzenia × kanały × próbki
%   stimulationFrequency  - częstość stymulacji w Hz
%   numberOfHarmonics     - liczba harmonicznych do uwzględnienia
%   Fs                    - częstość próbkowania w Hz
%
% Wyjście:
%   W      - macierz filtrów przestrzennych (kanały × kanały), filtry w kolumnach
%   Lambda - macierz diagonalna wartości własnych

[numberOfTrials, numberOfChannels, numberOfSamples] = size(signal);

%% KROK 1: Zbuduj macierz referencyjną S (sinusy i cosinusy harmonicznych)
% S ma wymiary: próbki × 2*numberOfHarmonics
% Każda kolumna to jeden wersor (unormowany sinus lub cosinus).
% Wróć do sekcji Teoria powyżej i przełóż tamten kod na tę funkcję.

% TUTAJ

%% KROK 2: Oblicz średnie macierze kowariancji C_SSEP i C_Y
% Dla każdego powtórzenia wytnij sygnał x (kanały × próbki),
% rozłóż go na składową SSEP i resztę Y zgodnie z równaniami z sekcji Teoria,
% oblicz macierze kowariancji obu składowych, znormalizuj przez ślad
% i uśrednij po powtórzeniach.

C_SSEP = zeros(numberOfChannels, numberOfChannels);
C_Y    = zeros(numberOfChannels, numberOfChannels);

for trial = 1:numberOfTrials
    x = squeeze(signal(trial, :, :));   % kanały × próbki

    % TUTAJ: oblicz SSEP i Y

    % TUTAJ: oblicz kowariancje, znormalizuj, dodaj do C_SSEP i C_Y
end

% TUTAJ: podziel C_SSEP i C_Y przez numberOfTrials

%% KROK 3: Rozwiąż uogólnione zagadnienie własne
% Znajdź filtry przestrzenne maksymalizujące stosunek mocy SSEP do mocy Y.
% Analogia z CSP: jakie macierze wstawiłeś do eig() tam?

% TUTAJ: [W, Lambda] = ...

end

Wskazówka do interpretacji wyników

Filtr związany z największą wartością własną daje komponent o największym stosunku mocy SSEP do mocy tła. Sprawdź na widmach Welcha czy odpowiada on składowej z wyraźnym pikiem przy częstości stymulacji.

Rozwiązanie — pokaż dopiero po samodzielnej próbie

function [W, Lambda] = cosSinCSP(signal, stimulationFrequency, numberOfHarmonics, Fs)

[numberOfTrials, numberOfChannels, numberOfSamples] = size(signal);

%% KROK 1: Macierz referencyjna S
t = (0:1:numberOfSamples-1) / Fs;
S_ref = zeros(numberOfSamples, 2*numberOfHarmonics);

for harmonicNumber = 1:numberOfHarmonics
    c = cos(2*pi*stimulationFrequency*harmonicNumber*t);
    s = sin(2*pi*stimulationFrequency*harmonicNumber*t);
    S_ref(:,(harmonicNumber-1)*2 + 1) = c/norm(c);
    S_ref(:,(harmonicNumber-1)*2 + 2) = s/norm(s);
end

%% KROK 2: Średnie macierze kowariancji
C_SSEP = zeros(numberOfChannels, numberOfChannels);
C_Y    = zeros(numberOfChannels, numberOfChannels);

for trial = 1:numberOfTrials
    x    = squeeze(signal(trial, :, :));   % kanały × próbki
    A    = x * S_ref;                      % projekcja na przestrzeń SSEP
    SSEP = A * S_ref';                     % rekonstrukcja składowej SSEP
    Y    = x - SSEP;                       % reszta

    tmp = cov(SSEP');
    C_SSEP = C_SSEP + tmp/trace(tmp);

    tmp = cov(Y');
    C_Y = C_Y + tmp/trace(tmp);
end

C_SSEP = C_SSEP / numberOfTrials;
C_Y    = C_Y    / numberOfTrials;

%% KROK 3: Uogólnione zagadnienie własne
[W, Lambda] = eig(C_SSEP, C_Y);

end

Skrypt demonstracyjny i analiza wyników

% wczytujemy dane
load('PrzykladoweDaneSSVEP.mat');

[numberOfTrials numberOfChannels numberOfSamples] = size(X.data);
namesOfChannels = X.channels;

numberOfHarmonics = 3;
signal = X.data;   % powtórzenie × kanał × próbki

[W, Lambda] = cosSinCSP(signal, X.stimulation, numberOfHarmonics, X.sampling);

%% Odtworzenie sygnałów źródłowych
S = zeros(size(signal));
for powt = 1:size(signal,1)
    S(powt,:,:) = W' * squeeze(signal(powt,:,:));
end

%% Widma Welcha: kanały oryginalne vs źródła CSP
figure('Name', ['Stymulacja: ', num2str(X.stimulation), ' Hz'])
for i = 1:numberOfChannels
    subplot(2, numberOfChannels, i)
    PP = 0;
    for rep = 1:numberOfTrials
        x = squeeze(signal(rep,i,:));
        [Pxx, ff] = pwelch(x, X.sampling, 1, X.sampling, X.sampling);
        PP = PP + Pxx;
    end
    plot(ff(ff<60), PP(ff<60))
    title(namesOfChannels{i})

    subplot(2, numberOfChannels, numberOfChannels+i)
    PP = 0;
    for rep = 1:numberOfTrials
        s = squeeze(S(rep,i,:));
        [Pss, ff] = pwelch(s, X.sampling, 1, X.sampling, X.sampling);
        PP = PP + Pss;
    end
    plot(ff(ff<60), PP(ff<60))
    title(['źródło CSP: ', num2str(i)])
    xlabel(['\lambda = ', num2str(Lambda(i,i), '%.2f')])
end

%% KROK 4: Mapki topograficzne
% Zidentyfikuj komponent z największą wartością własną.
% Narysuj dwie mapki przy użyciu funkcji topoplot() z pakietu EEGLAB:
%
% (a) Filtr przestrzenny: wagi z jakimi kanały EEG składają się na ten komponent.
%     Wskazówka: filtr to kolumna macierzy W.
%
% (b) Topografia źródła: z jakimi wagami komponent dociera do poszczególnych elektrod.
%     Wskazówka: topografia zawarta jest w wierszach macierzy odwrotnej do W —
%     przypomnij sobie analogiczny krok z analizy CSP dla P300.
%
% Czym różnią się te dwie mapki? Którą łatwiej interpretować fizjologicznie?

% TUTAJ: [~, idx] = max(diag(Lambda));
% TUTAJ: topoplot(...)  % filtr
% TUTAJ: topoplot(...)  % topografia

Sesja 5: ICA — teoria i wydobywanie interesującego komponentu

Czas: ~150 min (jedno spotkanie)

Sesja wprowadza analizę składowych niezależnych (ICA) jako trzecie podejście do problemu ślepej separacji źródeł — obok BSS/CSP i cosSinCSP. Kluczowa różnica: CSP szuka kierunków maksymalizujących kontrast wariancji między warunkami (kryterium drugiego rzędu), ICA szuka składowych niezależnych statystycznie, maksymalizując niegaussowość (kryterium wyższego rzędu). ICA nie wymaga dwóch warunków eksperymentalnych — działa na jednym stacjonarnym sygnale.

Zakres wykładu (~40 min)

model generatywny ICA: [math]\mathbf{x} = \mathbf{D}\mathbf{s}[/math], założenie niegaussowości składowych
negoentropia i algorytm FastICA
porównanie z CSP: kiedy ICA, kiedy CSP
ograniczenia: liczba próbek a liczba kanałów ([math]kN^2[/math]), niejednoznaczność kolejności i skali komponentów

Zakres ćwiczenia (~80 min)

wczytanie danych do EEGLAB, edycja lokalizacji elektrod
filtrowanie górnoprzepustowe (FIR, 0,5 Hz), zmiana referencji
obliczenie ICA na całym sygnale
identyfikacja komponentów zawierających rytm alfa: topografia, widmo, przebieg czasowy
usunięcie komponentów niezawierających alfy, odtworzenie sygnału na elektrodach
porównanie wyników między co najmniej trzema uruchomieniami ICA: czy komponenty są powtarzalne?

Materiały

Zob. ICA jako filtr przestrzenny i Wydobywanie interesujących komponentów poniżej.

ICA jako filtr przestrzenny

Wstęp teoretyczny do ICA

Definicja

Independent Component Analysis (ICA) jest metodą statystycznej analizy sygnałów, która dokonuje dekompozycji wielokanałowych zapisów na składowe niezależne w sensie statystycznym. Dwie składowe s₁ i s₂ są niezależne jeżeli wiedza o wartości s₁ nie daje żadnych informacji o możliwych wartościach s₂. ICA może być wyrażona przez prosty model generatywny:

x = Ds

gdzie x = {x¹, x², ..., xⁿ} jest zmierzonym n kanałowym sygnałem, D jest macierzą mieszającą zaś s = {s¹, s², ..., sⁿ} jest aktywnością n źródeł. Podstawowym założeniem dotyczącym s jest to, że sⁱ są statystycznie niezależne. Aby wyestymować model musimy też założyć, że składowe mają niegaussowskie rozkłady wartości (Hyvärinen, 2000).

Dodatkowo model ten zakłada następujące fakty:

Sygnał jest liniową mieszaniną aktywności źródeł
Sygnały pochodzące z każdego ze źródeł są niezależne od pozostałych
Źródła oraz proces ich mieszania są stacjonarne, tzn, ich momenty statystyczne nie zależą od czasu
Energie (wariancje) źródeł nie mogą być wyznaczone jednoznacznie. Dzieje się tak ponieważ pomnożenie amplitudy i-tego źródła może być uzyskane poprzez przemnożenie albo sⁱ albo przez przemnożenie i-tej kolumny macierzy D. Naturalnym rozwiązaniem tej niejednoznaczności jest wprowadzenie konwencji, że komponenty są normowane tak aby miały wariancję 1: E[(sⁱ)²] = 1.
Kolejność komponentów jest dowolna. Bo jeśli w ten sam sposób zmienimy kolejność komponentów w s i kolumn w D to dostaniemy dokładnie ten sam sygnał x.

Głównym wyzwaniem w analizie ICA jest estymacja macierzy mieszającej D. Gdy jest ona znana to komponenty mogą być wyliczone w następujący sposób:

s = D⁻¹x

Estymacja

Znalezienie niezależnych komponentów może być rozważane w świetle Centralnego Twierdzenia Granicznego jako poszukiwanie komponentów o możliwie nie gaussowskim rozkładzie. Aby zrozumieć to podejście prześledźmy heurystykę zaproponowaną przez (Hyvärinen, 2000). Dla prostoty załóżmy, że poszukiwane źródła niezależne mają identyczne rozkłady. Zdefiniujmy

y = w^Tx.

Zauważmy, że jeśli

w

jest jedną z kolumn macierzy

D⁻¹,

to y jest jednym z poszukiwanych komponentów. Zamieniając zmienne

z = D^Tw

możemy napisać

y = w^Tx = w^TDs = z^Ts.

Uwidacznia to fakt, że y jest liniową kombinacją składowych sⁱ z wagami danymi przez z_i.

Z centralnego twierdzenia granicznego wynika, że suma niezależnych zmiennych losowych ma bardziej gaussowski charakter niż każda z tych zmiennych osobno. Liniowa kombinacja staje się najmniej gaussowska gdy z ma tylko jeden element niezerowy. W tym przypadku y jest proporcjonalny do sⁱ. Zatem problem estymacji modelu ICA może być sformułowany jako problem znalezienia takiego wektora w, który maksymalizuje niegaussowskość

y = w^Tx.

Maksymalizacja niegaussowskości y daje jeden niezależny komponent odpowiadający jednemu z 2n maksimów (bo mamy sⁱ i −sⁱ) w krajobrazie optymalizacyjnym. Aby znaleźć wszystkie niezależne komponenty musimy znaleźć wszystkie maksima. Ponieważ komponenty są nieskorelowane, to poszukiwania kolejnych komponentów można kontynuować w podprzestrzeni ortogonalnej do już znalezionych komponentów.

Obliczenia

Intuicyjna heurystyka poszukiwania najbardziej niegaussowskich składowych może być użyta do wyprowadzenia różnych funkcji kosztu, których optymalizacja daje model ICA, np. kurtoza.

[math]kurt(y) = E\{y^4\} - 3(E{y^2})^2[/math]

Inną miarą gassowskości jest neg-entropia, którą można wyprowadzić z entropii: Entropia jest miarą średniego zdziwienia wynikiem obserwacji zmiennej losowej: [math]H(Y) = - \sum_i P(Y= a_i) \log(P(Y=a_i)) [/math]

Negentropia jest zdefiniowana:

[math] J(y) = H(y_{gauss}) -H(y)[/math] gdzie [math] y_{gauss} [/math] jest gassuwską zmienną losową o takiej samej kowaiancji jak [math] y [/math].

Negentropia jest skomplikowana obliczeniow, więc w praktyce używana jest formuła przybliżona:

[math] J(y) \varpropto [E\{G (y)\} - E\{G(\nu)\}][/math]

[math]\nu [/math] jest zmienną losową ze standardowego rozkładu normalnego , a G są pewnymi niekwadratowymi funkcjami.

W algorytmie FastICA extremum negentropii jest znajdowane w procedurze bazującej na optymalizacji Netwona. (szczegóły np.: sekcja 6 w https://www.cs.helsinki.fi/u/ahyvarin/papers/NN00new.pdf)

Procedura wykorzystywana w eeglabie („runica”, Makeig 1996) dąży do minimalizacji informacji wzajemnej. Oba podejścia są w przybliżeniu równoważne (Hyvärinen, 2000), chociaż owo przybliżenie dla sygnałów elektrofizjologicznych nie zostało to jeszcze w pełni wyeksplorowane. Dla sygnałów o niskiej wymiarowości i spełniających dokładnie założenia ICA wszystkie powszechnie wykorzystywane algorytmy dają niemal identyczne wyniki.

Bardzo ważna uwaga: ogólną zasadą jest, że jeśli estymujemy N stabilnych komponentów (z N-kanałowych danych) to musimy dysponować kN² punktami danych w każdym kanale, gdzie N² jest liczbą elementów w macierzy D, którą ICA próbuje wyestymować, k jest liczbą całkowitą. Nie ma dobrych oszacowań teoretycznych na wielkość k, z praktycznych obserwacji wynika, że rośnie ona z liczbą kanałów.

Możliwe zastosowania

Najczęściej ICA jest stosowana jako narzędzie do:

usuwania artefaktów z sygnałów EEG (ruchy oczu i mięśnie)
wydobywania składowych do dalszej analizy (Onton, 2006)
jako analiza wstępna do lokalizacji źródeł (Grau, 2007).
ICA jest także stosowana w analize sygnałów EKG i EMG.

Bibliografia

Bazowa praca:

A. Hyvärinen. Fast and Robust Fixed-Point Algorithms for Independent Component Analysis. IEEE Transactions on Neural Networks 10(3):626-634, 1999 http://www.cs.helsinki.fi/u/ahyvarin/papers/TNN99_reprint.pdf

Nieco prościej opisana wersja z przykładami:

Hyvärinen, A. and Oja, E. (2000). Independent component analysis: Algorithms and applications. Neural Networks, 13(4-5):411–430.

https://www.cs.helsinki.fi/u/ahyvarin/papers/NN00new.pdf

Grau, C., Fuentemilla, L., Marco-Pallars, J. (2007). Functional neural dynamics underlying auditory event-related n1 and n1 suppression response. Neuroimage, 36(6):522–31.
Makeig, S., Bell, A., Jung, T.-P., Sejnowski,T. (1996). Independent component analysis of electroencephalographic data. W: Touretzky, D., Mozer, M., and Hasselmo, M., editors, Advances in Neural Information Processing Systems, volume 8, pages 145–151. MIT Press, Cambridge, MA.
Onton, J., Makeig, S. (2006). Information-based modeling of event-related brain dynamics. Prog Brain Res., 159:99–120.
Tutorial: http://sccn.ucsd.edu/wiki/Chapter_09:_Decomposing_Data_Using_ICA
http://sccn.ucsd.edu/~arno/indexica.html
http://cis.legacy.ics.tkk.fi/aapo/papers/IJCNN99_tutorialweb/

Wydobywanie interesujących komponentów

ZADANIE: Wydobywanie interesujących komponentów

Dane do tej części ćwiczeń proszę pobrać i rozpakować w swoim katalogu: http://www.fuw.edu.pl/~jarekz/LabEEG/Dane_do_ICA_alfa.tar.gz

Pochodzą one z eksperymentu, w którym osoba badana siedziała z zamkniętymi oczami słuchając nagrania czytanego spokojnym głosem. Metadane opisujące sygnał znajdują się w pliku Miro.xml, zaś lokalizacje elektrod w pliku Miro-10-20-Cap.locs.

Proszę:

wczytać dane do eeglaba
wyedytować lokalizację elektrod
usunąć kanały nie zawierające EEG
zmienić referencje na średnią z kanałów A1 i A2
przefiltrować filtrem FIR górnoprzepustowym z częstością odcięcia 0,5 Hz
obejrzeć wstępnie przygotowane dane
policzyć ICA na całym sygnale
obejrzeć właściwości otrzymanych komponentów
- Czy są wśród nich takie, które zawierają znaczny udział rytmu alfa?
- Jaka jest ich topografia?
usunąć wszystkie komponenty nie zawierające alfy
odtworzyć z tych komponentów sygnał na elektrodach
wykonać dekompozycję ICA kilkukrotnie (co najmniej 3) i porównać wyniki
- Czy uzyskiwane komponenty są powtarzalne?
- Swoje wyniki porównać też z sąsiednimi grupami.

Sesja 6: ICA — artefakty i synteza metod

Czas: ~150 min (jedno spotkanie)

Pierwsza część sesji poświęcona jest praktycznemu zastosowaniu ICA do usuwania artefaktów z danych EEG przy użyciu automatycznych klasyfikatorów komponentów (ICLabel, MARA). Druga część syntetyzuje całość materiału z bloku filtrów przestrzennych: BSS/CSP, cosSinCSP i ICA są trzema odpowiedziami na ten sam ogólny problem separacji źródeł, różniącymi się założeniami i kryterium optymalizacji.

Zakres ćwiczenia — artefakty (~75 min)

wczytanie danych Arousal do EEGLAB, usunięcie kanałów nieEEG
obliczenie ICA, przegląd topografii komponentów
identyfikacja artefaktów ocznych i mięśniowych

 przy użyciu wtyczek ICLabel i MARA

porównanie sygnału przed i po usunięciu komponentów artefaktowych

Zakres syntezy (~45 min)

tabela porównawcza CSP / cosSinCSP / ICA:

 założenia, dane wejściowe, kryterium optymalizacji, kiedy stosować

związek filtrów przestrzennych z architekturami sieci neuronowych:

 warstwa depthwise conv w ShallowConvNet i EEGNet jako
 wyuczony odpowiednik macierzy W z CSP —
 szczegóły w CSP a sieci neuronowe poniżej

omówienie raportów, pytania

Materiały

Zob. Identyfikacja artefaktów poniżej.

ZADANIE: Identyfikacja artefaktów

Proszę pobrać dane:

Pochodzą one z eksperymentu w którym osoba badana czytała słowa o różnych właściwościach wzbudzania emocji.

wczytaj je do eeglaba
wczytaj lokalizację kanałów z pliku Arousal-10-20-Cap.locs
obejrzyj przebiegi czasowe
odrzuć kanał z diodą (21) i z GSR (20)
zrób dekompozycję ICA
obejrzyj topografię komponentów
zidentyfikuj komponenty odpowiadające mruganiu i aktywności mięśniowej.

UWAGA: Aktualnie do wykrywania komponentów artefaktowych warto posłużyć się wtyczkami do eeglaba dostępnymi przez stronę:

https://sccn.ucsd.edu/eeglab/plugin_uploader/plugin_list_all.php

ICLabel
MARA

W raporcie:

zaprezentuj fragmenty sygnału zawierającego artefakty oczne i mięśniowe przed i po zastosowaniu czyszczenia poprzez usuwanie komponentów zdominowanych przez artefakty.
zaprezentuj topografię i przebiegi czasowe komponentów zidentyfikowanych jako artefakty oczne i mięśniowe.

-->

Anonimowy

Szukaj

Laboratorium EEG/CSP: Różnice pomiędzy wersjami

Przestrzenie nazw

Więcej

Działania na stronie

Aktualna wersja na dzień 11:05, 14 kwi 2026

Spis treści

Plan zajęć

Sesja 1: Ślepa separacja źródeł i CSP — wprowadzenie

Ćwiczenie symulacyjne

Zadania do samodzielnej realizacji

Szkielet kodu do uzupełnienia

Sesja 2: CSP dla P300 — dane realne

Sesja 3: CSP dla P300 — interpretacja i separacja cech

Analiza wstępna

ZADANIE: Analiza CSP

Wybór i separacja cech

Sesja 4: Filtry przestrzenne dla SSEP — cosSinCSP

Filtry przestrzenne dla SSEP

Teoria

Zadanie: implementacja funkcji cosSinCSP

Skrypt demonstracyjny i analiza wyników

Sesja 5: ICA — teoria i wydobywanie interesującego komponentu

ICA jako filtr przestrzenny

Definicja

Estymacja

Obliczenia

Możliwe zastosowania

Bibliografia

ZADANIE: Wydobywanie interesujących komponentów

Sesja 6: ICA — artefakty i synteza metod

ZADANIE: Identyfikacja artefaktów

W raporcie:

Nawigacja

Nawigacja

Narzędzia Wiki

Narzędzia Wiki

@@ Linia 1: / Linia 1: @@
 [[Laboratorium_EEG]]/BSS
-=Prezentacja=
+=Plan zajęć=
-[https://brain.fuw.edu.pl/edu/images/2/2f/BSS.pdf slajdy z prezentacji]
+Materiał realizowany jest w ciągu 3 tygodni (6 spotkań po ~2,5 h).
+{| class="wikitable"
+|-
+! Sesja !! Temat !! Kluczowe sekcje na tej stronie
+|-
+| '''1''' || Wykład BSS + CSP; ćwiczenie symulacyjne
+| [[#Ślepa separacja źródeł (BSS)|Ślepa separacja źródeł]], [[#Common Spatial Pattern|Common Spatial Pattern]], [[#Ćwiczenie symulacyjne|Ćwiczenie symulacyjne]]
+|-
+| '''2'''|| Dane P300: wczytanie, cięcie epok, wizualizacja ERP, implementacja CSP
+| [[#Analiza wstępna|Analiza wstępna]], [[#ZADANIE: Analiza CSP|Analiza CSP]]
+|-
+| '''3'''|| Mapki topograficzne; separacja cech
+| [[#ZADANIE: Analiza CSP|Analiza CSP]] (cd.), [[#Wybór i separacja cech|Wybór i separacja cech]]
+|-
+| '''4'''|| cosSinCSP jako uogólnienie CSP; dane SSVEP
+| [[#Filtry przestrzenne dla SSEP|Filtry przestrzenne dla SSVEP]]
+|-
+| '''5'''|| Wykład ICA; komponenty alfa
+| [[#ICA jako filtr przestrzenny|ICA jako filtr przestrzenny]], [[#ZADANIE: Wydobywanie interesujących komponentów|Komponenty alfa]]
+|-
+| '''6'''|| Artefakty (ICLabel/MARA); synteza CSP/cosSinCSP/ICA
+| [[#ZADANIE: Identyfikacja artefaktów|Identyfikacja artefaktów]]
+|}
+=Sesja 1: Ślepa separacja źródeł i CSP — wprowadzenie=
+;Czas: ~60 min wykład + ~90 min ćwiczenie symulacyjne
+Zajęcia wprowadzają do problemu filtracji przestrzennej sygnałów EEG.
+Punktem wyjścia jest ogólny problem ślepej separacji źródeł (BSS),
+z którego CSP wyłania się jako szczególny przypadek dla dwóch warunków
+eksperymentalnych. Filtry przestrzenne omawiane na tych zajęciach stanowią
+fundament klasycznych metod BCI (P300, SSVEP, motor imagery) oraz są
+bezpośrednim odpowiednikiem warstw przestrzennych w sieciach neuronowych
+takich jak ShallowConvNet i EEGNet — do tej analogii wrócimy w sesji 6.
+[https://drive.google.com/file/d/1WrXgWWDJ7g5ZdYltm_8s2NXbGamw4R8N/view?usp=sharing Slajdy do wykładu (PDF)]
+;Zakres wykładu:
+* model generatywny EEG: <math>x(t) = As(t)</math>, macierz kowariancji, idea diagonalizacji
+* BSS jako ogólny problem separacji źródeł
+* CSP jako szczególny przypadek: dwa warunki eksperymentalne, iloraz Rayleigha, uogólnione zagadnienie własne
+* interpretacja filtrów i topografii źródeł
+;Ćwiczenie symulacyjne:
+Zob. [[#Ćwiczenie symulacyjne|Ćwiczenie symulacyjne]] w sekcji CSP poniżej.
+====Ćwiczenie symulacyjne ====
+====Zadania do samodzielnej realizacji====
+Poniższy kod dostarcza gotowej infrastruktury: generacji sygnałów, macierzy mieszającej i wizualizacji.
+Zadaniem jest samodzielna implementacja kluczowych kroków analizy CSP.
+;Zadanie 1 (obowiązkowe): Obliczenie macierzy kowariancji
+Napisz funkcję <code>srednia_kowariancja(X, ind)</code>, która dla danej tablicy sygnałów <code>X</code>
+(wymiary: powtórzenia × kanały × próbki) i wektora indeksów czasowych <code>ind</code>
+zwraca średnią macierz kowariancji uśrednioną po powtórzeniach.
+Dla każdego powtórzenia zastosuj funkcję <code>cov</code> i znormalizuj wynik przez jego ślad (<code>trace</code>).
+Sprawdź wynik porównując z macierzami <code>R_B</code> i <code>R_E</code> obliczonymi w dostarczonym kodzie.
+;Zadanie 2 (obowiązkowe): Implementacja CSP
+Korzystając z funkcji <code>eig</code> oraz macierzy kowariancji obliczonych w zadaniu 1,
+wyznacz macierz filtrów przestrzennych <code>W</code> i odpowiadające im wartości własne.
+Odtwórz sygnały źródłowe dla wszystkich powtórzeń.
+Narysuj chmury punktów (amplituda źródła 1 vs amplituda źródła 2) osobno dla warunków
+baseline i ERP — przed transformacją CSP i po niej. Co zmieniło się w kształcie chmur?
+Wykreśl także wizualizację przebiegów czasowych sygnałów X i sygnałów S uzyskanych po rozmieszaniu.
+;Zadanie 3 (obowiązkowe): Interpretacja wartości własnych
+Wartości własne znajdują się na przekątnej macierzy <code>Lambda</code>.
+Który filtr najbardziej różnicuje warunki baseline i ERP?
+Jaki jest związek między wartością własną a separacją widoczną na wykresach punktowych?
+;Zadanie 4 (dla chętnych): Wpływ macierzy mieszającej
+W dostarczonym kodzie macierz <code>P</code> (filtr przestrzenny) wynosi:
+<pre>
+P = [1 2; 1.5 1.3]
+</pre>
+Zmień ją kolejno na macierz bliską jednostkowej (np. <code>P = [1 0.1; 0.1 1]</code>)
+oraz na macierz z dużymi elementami pozadiagonalnymi (np. <code>P = [1 5; 4 1]</code>).
+Jak zmienia się kształt chmury punktów przed transformacją CSP?
+Czy CSP poprawnie oddziela źródła w obu przypadkach?
+;Zadanie 5 (dla chętnych): Wpływ szumu
+W kodzie szum jest wyłączony: <code>n = 0*randn(size(tmp))</code>.
+Usuń mnożnik <code>0*</code> i zbadaj jak poziom szumu wpływa na separację źródeł.
+Przy jakim stosunku sygnału do szumu CSP przestaje skutecznie rozdzielać warunki?
+;Pytanie do dyskusji:
+Dostarczony kod używa <code>eig(R_E, R_B)</code>.
+Alternatywą jest <code>eig(R_E, (R_E+R_B)/2)</code>.
+Czym różnią się te dwie wersje? W jakich sytuacjach eksperymentalnych
+druga wersja mogłaby być bardziej uzasadniona?
+====Szkielet kodu do uzupełnienia====
+Poniższy szkielet zawiera gotową infrastrukturę symulacji.
+Uzupełnij brakujące fragmenty oznaczone komentarzem <code>% TUTAJ</code>
+zgodnie z zadaniami 1–3. Pełny kod referencyjny dostępny jest poniżej
+— zajrzyj do niego dopiero po samodzielnej próbie.
+<source lang="matlab">
+%% Parametry symulacji — nie modyfikuj tej sekcji
+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));
+P = [1 2; 1.5 1.3];   % filtr przestrzenny (prawdziwy)
+A = P^(-1);            % topografia źródeł
+for r = 1:N_rep
+    s1 = sin(2*pi*11*t + pi/2) + 0.02*randn(size(t));
+    s2 = exp(-((t-0.8)/0.05).^2) + 0.01*randn(size(t));
+    s(r,1,:) = s1;
+    s(r,2,:) = s2;
+    tmp = squeeze(s(r,:,:));
+    n = 0*randn(size(tmp));   % Zadanie 5: zmień 0* na 1* aby włączyć szum
+    X(r,:,:) = A*tmp + n;
+end
+baseline_ind = find(t < 0.5);
+ERP_ind      = find(t >= 0.5);
+R_B = zeros(N_chan, N_chan);
+R_E = zeros(N_chan, N_chan);
+%% Zadanie 1: Oblicz średnie macierze kowariancji
+% korzystając z funkcji srednia_kowariancja(X, ind)
+% Dla każdego powtórzenia r:
+%   - wytnij fragment sygnału X(r,:,baseline_ind) i analogicznie ERP
+%   - oblicz macierz kowariancji funkcją cov() — pamiętaj o transpozycji
+%   - znormalizuj przez ślad (funkcja trace())
+%   - dodaj do sumy
+% Na końcu podziel przez N_rep
-=Ślepa separacja źródeł (BSS)=
-{{hidden begin|title=Wstęp teoretyczny do BSS}}
-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.
+%% Zadanie 2: Rozwiąż uogólnione zagadnienie własne
-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.
+% Użyj funkcji eig(A,B) gdzie A=R_E, B=R_B
-[[Plik:Kowariancja.png|200px|center]]
+% W wyniku otrzymasz macierz wektorów własnych W (w kolumnach)
+% i macierz Lambda z wartościami własnymi na przekątnej
-==Filtry przestrzenne i ślepa separacja źródeł==
+% TUTAJ: [W, Lambda] = ...
-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, czyli przez cały czas aktywna jest ta sama konfiguracja źródeł niezależnych.
+disp('Wartości własne (przekątna Lambda):')
-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.
+disp(diag(Lambda))
-[[Plik:Diagonalizacja.png|200px|center]]
+%% Zadanie 2 cd.: Odtwórz sygnały źródłowe
+% Dla każdego powtórzenia r zastosuj filtr: S(r,:,:) = W' * X(r,:,:)
-==Common Spatial Pattern ==
+S = zeros(N_rep, N_chan, length(t));
-===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óry maksymalizuje stosunek mocy (wariancji) sygnałów rejestrowanych w dwóch rożnych warunkach eksperymentalnych.
-===Formalizm===
+for r = 1:N_rep
-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łą.
+    % TUTAJ: S(r,:,:) = ...
+end
-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>:
+%% Zadanie 3: Wizualizacja — chmury punktów przed i po CSP
-:<math> s_w(t) = w^T x(t)</math>
+x_base_k1 = X(:,1,baseline_ind);
-Wariancja tego rzutowanego sygnału to:
+x_base_k2 = X(:,2,baseline_ind);
-:<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>
+x_ERP_k1  = X(:,1,ERP_ind);
-Zatem znalezienie właściwego kierunku rzutowania można wyrazić jako szukanie maksimum wyrażenia <math> J(w) </math> (jest to tzw. iloraz Rayleigha):
+x_ERP_k2  = X(:,2,ERP_ind);
-: <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{ 2 C_{T} w \left(w^T  C_{NT} w\right)  -2C_{NT} w \left(w^T  C_{T} w \right)}{\left(w^T  C_{NT} w\right)^2} = \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 do rozwiązania tzw. uogólnione zagadnienie własne:
-:<math>       C_{T}w  =\frac{w^T  C_{T} w}{w^T  C_{NT} w}   C_{NT} w   </math>
-We wzorze tym liczba <math>\lambda=\frac{w^T  C_{T} w}{w^T  C_{NT} w}</math> spełniająca to równanie jest uogólnioną wartością własną, wtedy <math>w</math> jest uogólnionym wektorem własnym odpowiadającym tej wartości.
-Aby znaleźć <math> \lambda</math> i <math>w</math> możemy wykorzystać w Matlabie funkcję <tt>eig</tt>. Funkcja ta rozwiązuje (również) uogólnione zagadnienia własne postaci ''Aw''=&lambda;''Bw'' dostarczając w wyniku macierz wektorów własnych (w kolumnach) oraz macierz zawierającą na przekątnej odpowiadające im wartości własne.
+s_base_k1 = squeeze(S(:,1,baseline_ind));
-{{hidden end}}
+s_base_k2 = squeeze(S(:,2,baseline_ind));
+s_ERP_k1  = squeeze(S(:,1,ERP_ind));
+s_ERP_k2  = squeeze(S(:,2,ERP_ind));
-<!--
+figure(1); clf
-Odwzorowanie to można przedstawić w postaci macierzy <math>P</math>, której każdy wiersz zawiera wagi dla odpowiednich kanałów.
+subplot(1,2,1)
-Macierz zawierająca sygnał <math>X^{\pm}(t)</math>
+    plot(x_base_k1(:), x_base_k2(:), 'b.'); hold on
-ma wymiary <math>C \times  N</math>, gdzie
+    plot(x_ERP_k1(:),  x_ERP_k2(:),  'r.');
-<math>C</math> to liczba kanałów EEG, natomiast
+    axis equal; xlim([-2 2]); ylim([-2 2])
-<math>N</math> to liczba próbek dla każdego z kanałów.
+    xlabel('Kanał 1'); ylabel('Kanał 2')
-Macierz <math>P</math> przekształca sygnał <math>X^{\pm}(t)</math>  zgodnie ze wzorem:
+    title('Przed CSP')
-:<math>X^{\pm}_{CSP}(t)=P^T  X^{\pm}(t) </math>
+    legend('baseline','ERP')
-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 przekątnej można również zapisać za pomocą wzoru:
-:<math>k^{+}_{i} = \vec{p}^{T}_{i} R^{+} \vec{p}_{i}    </math> (2.6)
+subplot(1,2,2)
-:<math>k^{-}_{i} = \vec{p}^T_{i} R^{-}\vec{p}_{i} </math> (2.7)
+    plot(s_base_k1(:), s_base_k2(:), 'b.'); hold on
-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
+    plot(s_ERP_k1(:),  s_ERP_k2(:),  'r.');
-otrzymać równanie:
+    axis equal
+    xlabel('Źródło 1'); ylabel('Źródło 2')
+    title('Po CSP')
+    legend('baseline','ERP')
+%% Wizualizacja przebiegów czasowych X i S
+% TUTAJ:
-:<math>R^{+} \vec{p}_{i} = \frac{k^{+}_{i}}{k^{-}_{i}} R^{-} \vec{p}_{i} </math>  (2.8)
+%% Zadanie 4 (dla chętnych): zmień macierz P powyżej i powtórz analizę
+%% Zadanie 5 (dla chętnych): włącz szum (zmień 0* na 1* w linii z n=...)
-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.
+</source>
-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 ====
+{{hidden begin|title=Rozwiązanie — pokaż dopiero po samodzielnej próbie}}
-{{hidden begin|title=kod przykładowy}}
 <source lang  = matlab>
 % symulowany eksperyment składa się z sinusoidy udającej alfę spoczynkową i
@@ Linia 239: / Linia 354: @@
 {{hidden end}}
-==Zastosowanie filtra CSP do detekcji potencjału P300==
+=Sesja 2: CSP dla P300 — dane realne=
-{{hidden begin|title=Eksperyment}}
-===Eksperyment===
+;Czas: ~150 min (jedno spotkanie)
-* Proszę zapoznać się z instrukcją: http://laboratorium-eeg.braintech.pl/rozdz10.html
-* Proszę wczytać i uruchomić (na sucho) demo Demos->EEG_P300wz
+Pracujemy z danymi kalibracyjnymi z eksperymentu P300.
-* Wspólne omówienie konstrukcji i potencjalnych modyfikacji tego scenariusza
+Celem sesji jest przejście pełnego pipeline'u: od surowych danych
+do odtworzonych źródeł CSP.
+;Zakres:
+* wczytanie danych kalibracyjnych, identyfikacja znaczników T i NT
+* cięcie epok (−200 do +800 ms), usuwanie trendu liniowego
+* montaż względem połączonych uszu
+* wizualizacja ERP w układzie topograficznym
+* obliczenie macierzy kowariancji R<sub>T</sub> i R<sub>NT</sub>, normalizacja przez ślad
+* rozwiązanie uogólnionego zagadnienia własnego: <code>[W, Lambda] = eig(R_T, R_NT)</code>
+* odtworzenie sygnałów źródłowych: <code>S = W'*EEG</code>
+* przegląd estymowanych źródeł w dwóch warunkach
+;Materiały:
+Zob. [[#Analiza wstępna|Analiza wstępna]] i [[#ZADANIE: Analiza CSP|Analiza CSP]] poniżej.
+=Sesja 3: CSP dla P300 — interpretacja i separacja cech=
+;Czas: ~150 min (jedno spotkanie)
+Kontynuacja analizy CSP. Celem sesji jest interpretacja wyników
+w kategoriach fizjologicznych oraz ocena separowalności klas
+w przestrzeni cech opartej na mocy sygnału.
+;Zakres:
+* mapki topograficzne filtra przestrzennego i topografii źródła (funkcja <code>topoplot</code> z pakietu EEGLAB)
+* związek wartości własnych z siłą różnicowania warunków T i NT
+* wykresy punktowe mocy: oś X — moc składowej z największą wartością własną, oś Y — moc kolejnej składowej; jeden punkt = jedno powtórzenie
+* uśrednianie po 2, 4, 6, 8 i 10 realizacjach — obserwacja poprawy separacji wraz z liczbą uśrednień
+* dyskusja: utożsamianie komponentów CSP z fizjologicznymi źródłami jest ryzykowne; komponenty maksymalizują kontrast między warunkami, niekoniecznie odpowiadają izolowanym generatorom
+;Materiały:
+Zob. [[#ZADANIE: Analiza CSP|Analiza CSP]] i [[#Wybór i separacja cech|Wybór i separacja cech]] poniżej.
-====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 &bdquo;P300&rdquo;
-** &bdquo;P-Brain2013 Signal (with ID)&rdquo; jako np. &bdquo;Sygnal&rdquo;
-** &bdquo;P-Brain2013 Calibration p300&rdquo; jako &bdquo;kalibracjaP300&rdquo;
-** &bdquo;P-Brain 2013 p300&rdquo; jako &bdquo;Labirynt&rdquo;
-* 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:
-*** <tt>active_channels</tt>
-*** <tt>channel_names</tt>
-** dodajemy też linijkę: <tt>sampling_rate = 256</tt>
-* wchodzimy po kolei do katalogów: Sygnal_configs, kalibracjaP300_configs, Labirynt_configs i odnajdujemy plik switch_backup.ini. 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ę
-::: <tt>ignore_channels = DriverSaw;AmpSaw;PO7;PO8</tt>
-::: na
-::: <tt>ignore_channels = DriverSaw;AmpSaw;A1;A2</tt>
-::* oraz linię:
-::: <tt>montage_channels = PO7;PO8</tt>
-::: na
-::: <tt>montage_channels = A1;A2</tt>
-====Przeprowadzenie badania:====
-# Uruchom scenariusz &bdquo;Sygnał&rdquo;.
-# 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 &bdquo;Sygnal&rdquo;. 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 &bdquo;Sygnal&rdquo; 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: <tt>test1</tt>.
-# Uruchamiamy scenariusz &bdquo;kalibracjaP300&rdquo;. 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 &bdquo;Labirynt&rdquo;.
-# Danych z kalibracji potrzebować będziemy kilka zestawów.  Proszę powtórzyć kilkukrotnie scenariusz &bdquo;kalibracjaP300&rdquo;. Przed każdym uruchomieniem trzeba zmienić string w pliku <tt>file_id_name</tt> np. na <tt>test???</tt> gdzie <tt>???</tt> oznacza kolejne numery.
--->
-{{hidden end}}
 ===Analiza wstępna===
-Poszczególne etapy analizy proszę kodować w osobnych funkcjach. Funkcje te powinny być wywoływane z nadrzędnego skryptu, który powinien umożliwic wykoanie całości analiz.
-* Wczytać dane kalibracyjne do Matlaba i pociąć je na realizacje typu T &mdash; &bdquo;target&rdquo; (związane z wystąpieniami litery &bdquo;B&rdquo;) i NT &mdash; &bdquo;non-target&rdquo; (pozostałe litery) o długości &minus;200 do +800 ms wokół triggerów. Dla każdej realizacji odjąć trend liniowy.
-* Sygnał zmontować wzgl. &bdquo;połączonych uszu&rdquo; i wyświetlić średnie przebiegi dla warunku T i NT w układzie topograficznym &mdash; można wykorzystać w tym celu poniższy fragment kodu.
-{{hidden begin|title= fragment kodu do rysowania dwóch ERPów -- zamiast funkcji plottopo}}
-<source lang = matlab>
-%% stworzenie osi
-      % zakładam, że dane do rysowania są w dwóch strukturach EEG1 i EEG2 (struktury eeglab),
-      % te dane są podzielone na realizacje i mają kształt (kanały , czas, epoki)
-      % i że mamy już załadowane położenia elektrod we współrzędnych
-      % biegunowych
-      w = 0.125;
-      h = 0.125;
-      sc = 0.8;
-      figure()
-     for chanNum = 1:length(EEG1.chanlocs)
-         r =EEG1.chanlocs(chanNum).radius;
-         theta = EEG1.chanlocs(chanNum).theta;
-         x = r* sin(theta/180*pi)*sc+0.46;
-         y = r* cos(theta/180*pi)*sc+0.46;
-         ax(chanNum) = axes('Position', [x, y, w, h]);
-     end
-%%uśredniam EEG po powtórzeniach aby otrzymać ERP
-  ERP1 =squeeze( mean( EEG1.data,3)) ;
-  ERP2 =squeeze( mean( EEG2.data,3)) ;
-%% rysowanie uśrednionych potencjałów z dwóch struktur EEG
- for chanNum = 1:length(EEG1.chanlocs)
-     hold(ax(chanNum),'on');
-     plot(ax(chanNum), EEG1.times, ERP1(chanNum ,:));
-     plot(ax(chanNum), EEG2.times, ERP2(chanNum ,:));
-     title(ax(chanNum),EEG1.chanlocs(chanNum).labels );
-     hold(ax(chanNum),'off');
- end
-</source>
-{{hidden end}}
-<small>Poniżej zaprezentowany jest przykładowy skrypt do cięcia danych wokół znaczników. Działa on z plikami zawartymi w archiwum:
+<source lang="matlab">
-: [[Plik:KalibracjaP300.tar.gz]]
+%% KROK 1: Wczytanie danych
-Korzysta z funkcji pomocniczych dostępnych w dystrybucji obci w katalogu
+% ReadManager to pomocnicza klasa do odczytu plików w formacie OpenBCI.
-: /usr/share/openbci/analysis/matlab_obci_signal_processing
+% Przyjmuje trzy pliki: XML z metadanymi, RAW z próbkami i TAG ze znacznikami zdarzeń.
-Openbci można pobrać z https://github.com/BrainTech/openbci</small>
-{{hidden begin|title=przykładowy skrypt}}
-<source lang = matlab>
-% ustalamy nzawy plików z danymi
 nazwaPliku = 'p_6301423_calibration_p300.obci';
-nameOfXMLFile = strcat(nazwaPliku,'.xml');
+nameOfXMLFile  = strcat(nazwaPliku, '.xml');
-nameOfTagFile = strcat(nazwaPliku,'.tag'); %tagi = znaczniki zdarzeń
+nameOfTagFile  = strcat(nazwaPliku, '.tag');
-namesOfDataFiles = strcat(nazwaPliku,'.raw');
+namesOfDataFiles = strcat(nazwaPliku, '.raw');
-% inicjujemy obiekt rm
+rm = ReadManager(nameOfXMLFile, namesOfDataFiles, nameOfTagFile);
-rm = ReadManager(nameOfXMLFile,namesOfDataFiles,nameOfTagFile);
-% obieramy przydatne parametry i znaczniki
+% Pobieramy podstawowe parametry nagrania
 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
+% Pobieramy listę wszystkich znaczników zdarzeń z nagrania
-numberOfStruct = length(tagsStruct);
+tagsStruct = rm.get_tags();
-targetTimeStamps = [];
+%% KROK 2: Identyfikacja znaczników Target i Non-Target
+% Każdy znacznik typu 'blink' odpowiada jednemu mignięciu w eksperymencie P300.
+% Pole 'index' mówi który symbol mignął, pole 'target' który symbol był oczekiwany.
+% Jeśli index == target to było to mignięcie docelowe (Target), w przeciwnym razie
+% Non-Target.
+targetTimeStamps    = [];
 NonTargetTimeStamps = [];
-for structNumber = 1:numberOfStruct % iterujemy się przez tagi
-     if(strcmp(tagsStruct(structNumber).name,'blink')) % szukamy tagów o nazwie 'blink'
+for structNumber = 1:length(tagsStruct)
-         index = tagsStruct(structNumber).children.index; % tu jest numer pola stymulacji
+     if strcmp(tagsStruct(structNumber).name, 'blink')
-         target= tagsStruct(structNumber).children.target;% tu jest numer pola na którym wyświetlany jest target
+         index  = tagsStruct(structNumber).children.index;
-         if index == target % warunek na to, że mamy do czynienia z tagiem target
+         target = tagsStruct(structNumber).children.target;
-             targetTimeStamps = [targetTimeStamps tagsStruct(structNumber).start_timestamp]; %dodajemy timeStamp do listy targetów
+         if index == target
+             targetTimeStamps    = [targetTimeStamps    tagsStruct(structNumber).start_timestamp];
          else
-             NonTargetTimeStamps = [NonTargetTimeStamps tagsStruct(structNumber).start_timestamp];%dodajemy timeStamp do listy non-targetów
+             NonTargetTimeStamps = [NonTargetTimeStamps tagsStruct(structNumber).start_timestamp];
          end
      end
 end
+%% KROK 3: Wczytanie i filtrowanie sygnału
+samples = double(rm.get_samples());
+samples = samples(1:8,:);   % odrzucamy kanały bez EEG
+numberOfChannels = 8;
+% Filtr dolnoprzepustowy Czebyszewa 2. rodzaju — odcięcie 25 Hz
-% pobieramy próbki
+[b, a] = cheby2(6, 80, 25/(samplingFrequency/2), 'low');
-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,:));
+     samples(ch,:) = filtfilt(b, a, samples(ch,:));
 end
-% montujemy dane do wspólnej średniej (common average)
+%% KROK 4: Montaż względem średniej (common average reference)
-M = -ones(8,8)/8;
+% Odejmujemy od każdego kanału średnią ze wszystkich kanałów w danej chwili.
-M=M+eye(8,8)*9/8;
+% Jest to prostsze i bardziej ogólne niż montaż względem połączonych uszu.
-samples = 0.0715*M*samples;
+M = -ones(numberOfChannels, numberOfChannels) / numberOfChannels;
+M = M + eye(numberOfChannels) * (numberOfChannels+1) / numberOfChannels;
+samples = 0.0715 * M * samples;
+%% KROK 5: Cięcie epok
+PRE  = -0.2;
+POST =  0.8;
+wycinek = floor(PRE*samplingFrequency : POST*samplingFrequency);
-% wycinamy dane wokół znaczników
+TargetSignal    = zeros(length(targetTimeStamps),    numberOfChannels, length(wycinek));
-PRE = -0.2; % czas przed tagiem w sek.
+NonTargetSignal = zeros(length(NonTargetTimeStamps), numberOfChannels, length(wycinek));
-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);
+     trigerOnset = floor(targetTimeStamps(trialNumber) * samplingFrequency);
-     tenWycinek = wycinek + trigerOnset;
+     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ł
+     if tenWycinek(1) > 0 && tenWycinek(end) <= size(samples,2)
-         tmpSignal = samples(:,tenWycinek);
+         tmpSignal = samples(:, tenWycinek);
-         tmpSignal = detrend(tmpSignal')'; % usuwanie liniowego trendu - przy krótkich wycinkach działa lepiej niż filtrowanie górnoprzepustowe
+         tmpSignal = detrend(tmpSignal')';   % usuwanie trendu liniowego
-         TargetSignal(trialNumber, :,:) = tmpSignal;
+         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);
+     trigerOnset = floor(NonTargetTimeStamps(trialNumber) * samplingFrequency);
-     tenWycinek = wycinek + trigerOnset;
+     tenWycinek  = wycinek + trigerOnset;
-     if tenWycinek(1)>0 && tenWycinek(end)<=size(samples,2)
+     if tenWycinek(1) > 0 && tenWycinek(end) <= size(samples,2)
-         tmpSignal = samples(:,tenWycinek);
+         tmpSignal = samples(:, tenWycinek);
          tmpSignal = detrend(tmpSignal')';
-         NonTargetSignal(trialNumber, :,:) = tmpSignal;
+         NonTargetSignal(trialNumber,:,:) = tmpSignal;
      end
 end
-%
-% dla ilustracji podglądamy średnie po powtórzeniach ze wszystkich targetów
+% Oś czasu w milisekundach — przyda się do rysowania
-% i non-targetów
+times = (wycinek / samplingFrequency) * 1000;
-plot(squeeze(mean(TargetSignal,1))','r');
-hold on
+%% KROK 6: Wizualizacja ERP w układzie topograficznym
-plot(squeeze(mean(NonTargetSignal,1))','b')
+% Oblicz potencjały wywołane (ERP) jako średnie po powtórzeniach
+ERP_T  = squeeze(mean(TargetSignal,    1));   % wymiary: kanały × czas
+ERP_NT = squeeze(mean(NonTargetSignal, 1));
+% TUTAJ: narysuj ERP_T i ERP_NT dla każdego kanału w subplotach
+% ułożonych w siatce odpowiadającej przybliżonym pozycjom elektrod.
+% Na każdym subplocie nałóż oba warunki różnymi kolorami.
+% Oś X: czas w ms (użyj wektora times), oś Y: amplituda w µV.
+% Zadbaj o legendę i tytuły z nazwami kanałów (namesOfChannels).
+%% KROK 7: Obliczenie średnich macierzy kowariancji
+% Dla każdego powtórzenia oblicz macierz kowariancji funkcją cov(),
+% znormalizuj przez jej ślad (trace()), a następnie uśrednij po powtórzeniach.
+% Pamiętaj: cov() oczekuje macierzy próbki × kanały, czyli potrzebujesz transpozycji.
+R_T  = zeros(numberOfChannels, numberOfChannels);
+R_NT = zeros(numberOfChannels, numberOfChannels);
+for r = 1:length(targetTimeStamps)
+    % TUTAJ: wytnij r-te powtórzenie T, oblicz kowariancję, znormalizuj, dodaj do R_T
+end
+for r = 1:length(NonTargetTimeStamps)
+    % TUTAJ: wytnij r-te powtórzenie NT, oblicz kowariancję, znormalizuj, dodaj do R_NT
+end
+% TUTAJ: podziel R_T i R_NT przez odpowiednie liczby powtórzeń
+%% KROK 8: Rozwiązanie uogólnionego zagadnienia własnego
+% Szukamy filtrów przestrzennych W maksymalizujących stosunek mocy
+% w warunku T względem NT: eig(R_T, R_NT)
+% W kolumnach W znajdują się filtry, na przekątnej Lambda — wartości własne.
+% TUTAJ: [W, Lambda] = ...
+disp('Wartości własne (przekątna Lambda):')
+disp(diag(Lambda))
+%% KROK 9: Odtworzenie sygnałów źródłowych
+% Dla każdego powtórzenia zastosuj filtr przestrzenny:
+% s(r,:,:) = W' * x(r,:,:)
+% Pamiętaj o wymiarach — squeeze() może być pomocne.
+S_T  = zeros(size(TargetSignal));
+S_NT = zeros(size(NonTargetSignal));
+for r = 1:size(TargetSignal, 1)
+    % TUTAJ: S_T(r,:,:) = ...
+end
+for r = 1:size(NonTargetSignal, 1)
+    % TUTAJ: S_NT(r,:,:) = ...
+end
+%% KROK 10: Wizualizacja estymowanych źródeł
+% Narysuj średnie przebiegi źródeł (ERP w przestrzeni CSP) dla obu warunków.
+% Ułóż subploty w prostokątnej siatce — jeden subplot na każde źródło.
+% Dla którego źródła różnica między T i NT jest największa?
+% Czy odpowiada to źródłu z największą wartością własną?
+ERP_S_T  = squeeze(mean(S_T,  1));
+ERP_S_NT = squeeze(mean(S_NT, 1));
+% TUTAJ: narysuj ERP_S_T i ERP_S_NT w subplotach
 </source>
-{{hidden end}}
 ===ZADANIE: Analiza CSP===
@@ Linia 430: / Linia 560: @@
 --->
 * Wykonać 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?
+* 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 (do wykonania tych mapek wykorzystać funkcję <tt>topoplot</tt> z pakietu <tt>eeglab</tt>) wektorów odpowiadających:
+* Dla kanału najbardziej różnicującego wykonać mapki topograficzne wektorów odpowiadających:
 ** filtrowi przestrzennemu
 ** rzutu topograficznego źródła na elektrody.
-<!---
-* Zbadać powtarzalność topografii pomiędzy plikami konfiguracyjnymi.
---->
-{{hidden begin|title=Wybór i separacja cech}}
 ===Wybór i separacja cech===
@@ Linia 451: / Linia 576: @@
 {{hidden end}}
-==Filtry przestrzenne dla SSEP ==
+=Sesja 4: Filtry przestrzenne dla SSEP — cosSinCSP=
-{{hidden begin|title=Filtry przestrzenne dla SSEP}}
-=== Teoria===
-Ciekawa koncepcja filtra przestrzennego dla SSVEP zaprezentowana jest  tu: http://www.eurasip.org/Proceedings/Eusipco/Eusipco2009/contents/papers/1569193209.pdf
-Pokrótce można ją rozumieć podobnie do tego co robiliśmy rozważając filtry przestrzenne CSP z tym, że dla SSVEP oraz innych potencjałów wywołanych stanu ustalonego możemy skorzystać z dodatkowych informacji dotyczących poszukiwanych źródeł. Wiemy mianowicie, że powinny one oscylować z częstością bodźca, i być może jej harmonicznych.
+;Czas: ~150 min (jedno spotkanie)
-Przyda nam się macierz <math>S</math> zbudowana tak, że w kolejnych kolumnach znajdują się sinusy i cosinusy kolejnych częstości harmonicznych. Wektory te unormujemy, żeby miały energię równą 1. Innymi słowy macierz <math>S</math> zbudowana jest z wersorów rozpinających przestrzeń, w której powinien znajdować się sygnał SSVEP.
+Sesja wprowadza metodę cosSinCSP jako naturalne uogólnienie CSP
+o wiedzę dziedzinową dotyczącą sygnału SSEP. W klasycznym CSP
+dwa warunki eksperymentalne definiują dwie macierze kowariancji.
+W cosSinCSP rolę jednej z macierzy przejmuje projekcja sygnału
+na przestrzeń sinusów i cosinusów częstości stymulacji i jej harmonicznych.
+Formalizm ilorazu Rayleigha i uogólnionego zagadnienia własnego
+pozostaje identyczny.
-W matlabie możemy taką macierz zbudować tak:
+Dane używane na tych zajęciach pochodzą z eksperymentu SSEP
-<source lang = matlab>
+z archiwum — bez sesji rejestracyjnej. Na wykresach widma
-% Fs - częstość próbkowania
+sygnału przed transformacją CSP dominuje artefakt sieciowy 50 Hz,
-% numberOfSamples - długość sygnału w próbkach
+który praktycznie zasłania odpowiedź na stymulację. Zadaniem
+jest sprawdzenie czy cosSinCSP potrafi ten sygnał wydobyć.
+;Zakres:
+* idea cosSinCSP: macierz S złożona z sinusów i cosinusów harmonicznych, projekcja sygnału na przestrzeń SSEP i resztę Y
+* implementacja funkcji <code>cosSinCSP</code> w MATLAB
+* analiza danych archiwalnych: widma Welcha dla kanałów oryginalnych i źródeł CSP — porównanie SNR przed i po transformacji
+* prezentacja filtrów przestrzennych i topografii źródeł
+;Materiały:
+Zob. [[#Filtry przestrzenne dla SSEP|Filtry przestrzenne dla SSEP]] poniżej.
+==Filtry przestrzenne dla SSEP==
+===Teoria===
+Ciekawa koncepcja filtra przestrzennego dla SSEP zaprezentowana jest tu: http://www.eurasip.org/Proceedings/Eusipco/Eusipco2009/contents/papers/1569193209.pdf
+Pokrótce można ją rozumieć podobnie do tego co robiliśmy rozważając filtry przestrzenne CSP z tym, że dla SSEP oraz innych potencjałów wywołanych stanu ustalonego możemy skorzystać z dodatkowych informacji dotyczących poszukiwanych źródeł. Wiemy mianowicie, że powinny one oscylować z częstością bodźca i być może jej harmonicznych.
+Przyda nam się macierz <math>S</math> zbudowana tak, że w kolejnych kolumnach znajdują się sinusy i cosinusy kolejnych częstości harmonicznych. Wektory te unormujemy, żeby miały energię równą 1. Innymi słowy macierz <math>S</math> zbudowana jest z wersorów rozpinających przestrzeń, w której powinien znajdować się sygnał SSEP.
+W Matlabie możemy taką macierz zbudować tak:
+<source lang="matlab">
+% Fs                - częstość próbkowania
+% numberOfSamples   - długość sygnału w próbkach
 % numberOfHarmonics - liczba harmonicznych, które chcemy włączyć do analizy
-t = (0:1:numberOfSamples - 1)/Fs;
+t = (0:1:numberOfSamples - 1) / Fs;
 S = zeros(numberOfSamples, 2*numberOfHarmonics);
@@ Linia 476: / Linia 630: @@
 </source>
+Aby w badanym sygnale znaleźć składowe odpowiadające SSEP musimy rzutować sygnał <math>X</math> (macierz sygnałów ''kanały &times; próbki'') na przestrzeń rozpiętą przez <math>S</math>:
+:<math>A = X S</math>
+Macierz <math>A</math> zawiera współczynniki będące iloczynami skalarnymi sygnałów i wersorów — mówią one o tym jak dużo sinusa bądź cosinusa o danej częstości jest w pierwotnym sygnale. Komponenty SSEP zawarte w sygnale <math>X</math> odzyskujemy tak:
+:<math>\mathrm{SSEP} = A S^T</math>
+Modelujemy rejestrowany sygnał jako:
+:<math>X = \mathrm{SSEP} + Y</math>
+gdzie:
+:<math>Y = X - \mathrm{SSEP}</math>
+to wszystkie komponenty sygnału, które nas nie interesują.
+Filtr przestrzenny, który chcemy zbudować powinien maksymalizować stosunek wariancji <math>\mathrm{SSEP}</math> do wariancji <math>Y</math>. Macierz kowariancji powinna być uśredniona po powtórzeniach, a kowariancja sygnału w każdym powtórzeniu powinna być znormalizowana przez jej ślad. Dalej stosujemy technikę znaną z CSP: maksymalizację ilorazu Rayleigha przez rozwiązanie uogólnionego zagadnienia własnego dla macierzy kowariancji <math>\mathrm{SSEP}</math> i <math>Y</math>.
+===Zadanie: implementacja funkcji cosSinCSP===
+Dane do zadania: [[Plik:PrzykladoweDaneSSVEP.mat.gz]].
+Pełny przykład z danymi z eksperymentu SSEP: [[Plik:SSVEP_demo_csp.tar.gz]]
+Zaimplementuj funkcję <code>cosSinCSP</code> zgodnie z poniższym szkieletem.
+Prawidłowo zaimplementowana funkcja wraz ze skryptem demonstracyjnym poniżej
+powinna generować rysunek: [[Plik:Rys_SSVEP_demo.png|400px|podpis grafiki]]
+<source lang="matlab">
+function [W, Lambda] = cosSinCSP(signal, stimulationFrequency, numberOfHarmonics, Fs)
+% cosSinCSP - filtr przestrzenny dla sygnałów SSEP
+%
+% Wejście:
+%   signal                - dane EEG, wymiary: powtórzenia × kanały × próbki
+%   stimulationFrequency  - częstość stymulacji w Hz
+%   numberOfHarmonics     - liczba harmonicznych do uwzględnienia
+%   Fs                    - częstość próbkowania w Hz
+%
+% Wyjście:
+%   W      - macierz filtrów przestrzennych (kanały × kanały), filtry w kolumnach
+%   Lambda - macierz diagonalna wartości własnych
+[numberOfTrials, numberOfChannels, numberOfSamples] = size(signal);
+%% KROK 1: Zbuduj macierz referencyjną S (sinusy i cosinusy harmonicznych)
+% S ma wymiary: próbki × 2*numberOfHarmonics
+% Każda kolumna to jeden wersor (unormowany sinus lub cosinus).
+% Wróć do sekcji Teoria powyżej i przełóż tamten kod na tę funkcję.
+% TUTAJ
+%% KROK 2: Oblicz średnie macierze kowariancji C_SSEP i C_Y
+% Dla każdego powtórzenia wytnij sygnał x (kanały × próbki),
+% rozłóż go na składową SSEP i resztę Y zgodnie z równaniami z sekcji Teoria,
+% oblicz macierze kowariancji obu składowych, znormalizuj przez ślad
+% i uśrednij po powtórzeniach.
+C_SSEP = zeros(numberOfChannels, numberOfChannels);
+C_Y    = zeros(numberOfChannels, numberOfChannels);
+for trial = 1:numberOfTrials
+    x = squeeze(signal(trial, :, :));   % kanały × próbki
+    % TUTAJ: oblicz SSEP i Y
+    % TUTAJ: oblicz kowariancje, znormalizuj, dodaj do C_SSEP i C_Y
+end
+% TUTAJ: podziel C_SSEP i C_Y przez numberOfTrials
+%% KROK 3: Rozwiąż uogólnione zagadnienie własne
+% Znajdź filtry przestrzenne maksymalizujące stosunek mocy SSEP do mocy Y.
+% Analogia z CSP: jakie macierze wstawiłeś do eig() tam?
+% TUTAJ: [W, Lambda] = ...
+end
+</source>
+;Wskazówka do interpretacji wyników:
+Filtr związany z największą wartością własną daje komponent
+o największym stosunku mocy SSEP do mocy tła.
+Sprawdź na widmach Welcha czy odpowiada on składowej
+z wyraźnym pikiem przy częstości stymulacji.
+{{hidden begin|title=Rozwiązanie — pokaż dopiero po samodzielnej próbie}}
+<source lang="matlab">
+function [W, Lambda] = cosSinCSP(signal, stimulationFrequency, numberOfHarmonics, Fs)
+[numberOfTrials, numberOfChannels, numberOfSamples] = size(signal);
+%% KROK 1: Macierz referencyjna S
+t = (0:1:numberOfSamples-1) / Fs;
+S_ref = zeros(numberOfSamples, 2*numberOfHarmonics);
+for harmonicNumber = 1:numberOfHarmonics
+    c = cos(2*pi*stimulationFrequency*harmonicNumber*t);
+    s = sin(2*pi*stimulationFrequency*harmonicNumber*t);
+    S_ref(:,(harmonicNumber-1)*2 + 1) = c/norm(c);
+    S_ref(:,(harmonicNumber-1)*2 + 2) = s/norm(s);
+end
+%% KROK 2: Średnie macierze kowariancji
+C_SSEP = zeros(numberOfChannels, numberOfChannels);
+C_Y    = zeros(numberOfChannels, numberOfChannels);
+for trial = 1:numberOfTrials
+    x    = squeeze(signal(trial, :, :));   % kanały × próbki
+    A    = x * S_ref;                      % projekcja na przestrzeń SSEP
+    SSEP = A * S_ref';                     % rekonstrukcja składowej SSEP
+    Y    = x - SSEP;                       % reszta
+    tmp = cov(SSEP');
+    C_SSEP = C_SSEP + tmp/trace(tmp);
+    tmp = cov(Y');
+    C_Y = C_Y + tmp/trace(tmp);
+end
-Aby w badanym sygnale znaleźć składowe odpowiadające SSVEP musimy rzutować sygnał <math>X</math> (macierz sygnałów ''kanały &times; próbki'')  na przestrzeń rozpiętą przez <math>S</math>:
+C_SSEP = C_SSEP / numberOfTrials;
-:<math>A = X*S</math>
+C_Y    = C_Y    / numberOfTrials;
-Macierz <math>A</math> zawiera współczynniki będące iloczynami skalarnymi sygnałów i wersorów. Mówią one o tym &bdquo;jak dużo&rdquo; jest sinusa bądź cosinusa o danej częstości w pierwotnym sygnale. Komponenty SSVEP zawarte w sygnale <math>X</math> odzyskujemy tak:
-:<math>\mathrm{SSVEP} = A S^T</math>
-Modelujemy rejestrowany sygnał jako:
+%% KROK 3: Uogólnione zagadnienie własne
-:<math>X = \mathrm{SSVEP} + Y </math>
+[W, Lambda] = eig(C_SSEP, C_Y);
-gdzie:
-:<math>Y = X-\mathrm{SSVEP}</math>
-: to wszystkie komponenty sygnału, które nas nie interesują.
-Filtr przestrzenny, który chcemy zbudować powinien maksymalizować stosunek wariancji <math>\mathrm{SSVEP} = A S^T</math> do wariancji <math>Y = X-\mathrm{SSVEP}</math>. Macierz kowariancji powinna być uśredniona po powtórzeniach a kowariancja sygnału w każdym powtórzeniu powinna być znormalizowana poprzez podzielenie przez jej ślad (Matlabowa funkcja <tt>cov</tt> już wykonuje tę operację).
+end
-Dalej możemy zastosować technikę znaną z konstrukcji filtrów CSP, tzn. maksymalizacji ilorazu Rayleigha za pomocą rozwiązania uogólnionego zagadnienia własnego dla macierzy kowariancji <math>\mathrm{SSVEP} </math> i <math>Y </math>.
+</source>
+{{hidden end}}
-===Poniżej prosta demonstracja dla danych zebranych EEG przy stymulacji SSVEP z częstotliwością 38 Hz.===
+===Skrypt demonstracyjny i analiza wyników===
-Spakowane dane: [[Plik:PrzykladoweDaneSSVEP.mat.gz]].
-W oparciu o powyższy opis proszę zaimplementować  funkcję <tt>cosSinCSP</tt>. Prawidłowo zaimplementowana funkcja wraz z poniższym kodem powinna generować rysunek:
-[[Plik:Rys_SSVEP_demo.png|400px|podpis grafiki]]
-Przykładowy skrypt i dane prezentujący konstrukcję i działanie tego typu filtrów przestrzennych dla pełnych danych z eksperymentu SSVEP: [[Plik:SSVEP_demo_csp.tar.gz]]
-{{hidden begin|title=przykładowy skrypt}}
+<source lang="matlab">
-<source lang = matlab>
 % wczytujemy dane
 load('PrzykladoweDaneSSVEP.mat');
@@ Linia 504: / Linia 761: @@
 [numberOfTrials numberOfChannels numberOfSamples] = size(X.data);
 namesOfChannels = X.channels;
-% numberOfChannels numberOfSamples numberOfTrials
-W = zeros(numberOfChannels,numberOfChannels);
 numberOfHarmonics = 3;
-signal = X.data; % (  powtórzenie,  kanał, próbki)
+signal = X.data;   % powtórzenie × kanał × próbki
+[W, Lambda] = cosSinCSP(signal, X.stimulation, numberOfHarmonics, X.sampling);
+%% Odtworzenie sygnałów źródłowych
 S = zeros(size(signal));
-W = cosSinCSP(signal,X.stimulation,numberOfHarmonics,X.sampling);
 for powt = 1:size(signal,1)
-     S(powt,:,:) = W'*squeeze(signal(powt,:,:));
+     S(powt,:,:) = W' * squeeze(signal(powt,:,:));
 end
-figure('Name',['Stymulacja: ',num2str(X.stimulation),' Hz'])
+%% Widma Welcha: kanały oryginalne vs źródła CSP
-for i =1:numberOfChannels
+figure('Name', ['Stymulacja: ', num2str(X.stimulation), ' Hz'])
-    % rysujemy widma uśrednione po realizacjach dla danych
+for i = 1:numberOfChannels
-    % z oryginalnych kanałów EEG
+     subplot(2, numberOfChannels, i)
-     subplot(2,8,i)
+     PP = 0;
-     PP=0;
      for rep = 1:numberOfTrials
          x = squeeze(signal(rep,i,:));
-         [Pxx,ff] = pwelch(x, X.sampling, 1, X.sampling, X.sampling);
+         [Pxx, ff] = pwelch(x, X.sampling, 1, X.sampling, X.sampling);
-         PP =PP + Pxx;
+         PP = PP + Pxx;
      end
-     plot(ff(ff<60),PP(ff<60))
+     plot(ff(ff<60), PP(ff<60))
      title(namesOfChannels{i})
-    % rysujemy widma uśrednione po realizacjach dla danych
+     subplot(2, numberOfChannels, numberOfChannels+i)
-    % z estymowanych źródeł CSP
+     PP = 0;
-     subplot(2,8,8+i)
-     PP=0;
      for rep = 1:numberOfTrials
          s = squeeze(S(rep,i,:));
-         [Pss,ff]=pwelch(s, X.sampling, 1, X.sampling, X.sampling);
+         [Pss, ff] = pwelch(s, X.sampling, 1, X.sampling, X.sampling);
-         PP =PP + Pss;
+         PP = PP + Pss;
      end
-     plot(ff(ff<60),PP(ff<60))
+     plot(ff(ff<60), PP(ff<60))
      title(['źródło CSP: ', num2str(i)])
+    xlabel(['\lambda = ', num2str(Lambda(i,i), '%.2f')])
 end
-</source>
-{{hidden end}}
+%% KROK 4: Mapki topograficzne
+% Zidentyfikuj komponent z największą wartością własną.
+% Narysuj dwie mapki przy użyciu funkcji topoplot() z pakietu EEGLAB:
+%
+% (a) Filtr przestrzenny: wagi z jakimi kanały EEG składają się na ten komponent.
+%     Wskazówka: filtr to kolumna macierzy W.
+%
+% (b) Topografia źródła: z jakimi wagami komponent dociera do poszczególnych elektrod.
+%     Wskazówka: topografia zawarta jest w wierszach macierzy odwrotnej do W —
+%     przypomnij sobie analogiczny krok z analizy CSP dla P300.
+%
+% Czym różnią się te dwie mapki? Którą łatwiej interpretować fizjologicznie?
-====ZADANIE: Analiza danych z eksperymentu własnego ====
+% TUTAJ: [~, idx] = max(diag(Lambda));
-# Przefiltruj sygnały EEG w paśmie 1-45 Hz za pomocą procedury <tt>filtfilt</tt>.
+% TUTAJ: topoplot(...)  % filtr
-# Na podstawie sygnału trigger oraz danych zapisanych w pliku wyodrębnij sygnały EEG zarejestrowane w trakcie stymulacji z odpowiednimi częstościami.
+% TUTAJ: topoplot(...)  % topografia
-# Uśrednij sygnały odpowiadające stymulacji tą samą częstością.
+</source>
-# Obejrzyj uśrednione sygnały. (Zaprezentuj je).
-#* Sposób I:
+=Sesja 5: ICA — teoria i wydobywanie interesującego komponentu=
-#**Dla każdej realizacji wyestymuj przy pomocy metody Welcha widmo mocy sygnału EEG.
-#**Dla każdej częstości stymulacji wyznacz poziom tła na podstawie widm pochodzących ze stymulacji innymi częstościami. Np. dla stymulacji częstością 10 Hz poziom tła można wyznaczyć jako 95 centyl ze zbioru wartości widma w częstości 10 Hz dla stymulacji częstościami {4, 7, 13, 16, 20, 25, 30, 35, 40} Hz.
-#**Dla każdej częstości stymulacji wyznacz miarę odpowiedzi SSVEP (amplitudę widma odpowiadającą częstości stymulacji powyżej poziomy tła dla pojedynczej próby).
-#**Zaprezentuj widma otrzymane przy stymulacjach różnymi częstościami wraz z korytarzem odpowwiadajacym 95% przedziałowi ufności wyznaczonemu dla widm z pojedynczych realizacji.
-#**Sporządź wykres odpowiedzi SSVEP od częstości z zaznaczeniem przedziałów ufności i poziomu tła.
-#* Sposób II:
-#**Wyestumuj filtr CSP-SSVEP wspólny dla wszystkich częstości stymulacji
-#**Powtórz kroki ze sposobu I dla uzyskanych komponentów.
-#**Zaprezentuj filtry przestrzenne i topografie dla poszczególnych komponentów
-#**Dla przypomnienia:
-#***Filtr to zestaw współczynników z jakimi należy zsumować sygnały z poszczególnych kanałów EEG aby dostać komponenty odpowiadające hipotetycznym źródłom nieskorelowanym. Filtr można zilustrować na głowie, przypisując poszczególnym pozycjom elektrod wagi równe współrzędnym wektora w (kolumna macierzy W).
-#***Topografia źródła to zestaw współczynników z jakimi docierają one do poszczególnych kanałów EEG. Topografia zawarta jest w wierszach macierzy odwrotnej do W.
-#**Można rysunki wykonać jako macierze 5x5 i w pozycji elektrody kolorujemy proporcjonalnie do współczynnika.
-===SSVEP-BCI===
+;Czas: ~150 min (jedno spotkanie)
-W zajęciach tych przydadzą nam się informacje z:
-https://brain.fuw.edu.pl/edu/index.php/Laboratorium_EEG/Wprowadzenie_do_syg_online
-====Wstęp====
-Naszym celem będzie stworzenie prostego interfejsu wykorzystującego zjawisko SSVEP.
-Po dwóch stronach monitora zamocujemy diody. Każda będzie migać ze swoją ustaloną częstością (np. 16 i 22 Hz) - warto zadbać aby nie były to częstości powiązane ze sobą harmonicznie, a z drugiej strony aby, biorąc pod uwagę wyniki poprzedniego zadania, dawały dobra odpowiedź SSVEP.
-Eksperyment będzie miał dwie cześci: sesję kalibracje i sesję on-line. Pomiędzy tymi sesjami będziemy uczyć kalsyfikator. Na podstawie części kalibracyjnej ustalimy jakie parametry przetwarznego on-line sygnału świadczą o patrzeniu się na diodę z lewej a jakie na tą z prawej strony ekranu.
+Sesja wprowadza analizę składowych niezależnych (ICA) jako trzecie
+podejście do problemu ślepej separacji źródeł — obok BSS/CSP
+i cosSinCSP. Kluczowa różnica: CSP szuka kierunków maksymalizujących
+kontrast wariancji między warunkami (kryterium drugiego rzędu),
+ICA szuka składowych niezależnych statystycznie, maksymalizując
+niegaussowość (kryterium wyższego rzędu). ICA nie wymaga
+dwóch warunków eksperymentalnych — działa na jednym stacjonarnym sygnale.
-Potem w sesji online będziemy porównywać (za pomocą predykcji klasyfikatora) rejestrowany sygnał z tymi warościami kalibracyjnymi i na tej podstawie system będzie zwracał informację o wyborze lewej lub prawej diody. To powinno umożliwić prostą komunikację na zasadzie pytanie i odpowiedź TAK/NIE.
+;Zakres wykładu (~40 min):
+* model generatywny ICA: <math>\mathbf{x} = \mathbf{D}\mathbf{s}</math>, założenie niegaussowości składowych
+* negoentropia i algorytm FastICA
+* porównanie z CSP: kiedy ICA, kiedy CSP
+* ograniczenia: liczba próbek a liczba kanałów (<math>kN^2</math>),  niejednoznaczność kolejności i skali komponentów
-==== Sesja kalibracycjna====
+;Zakres ćwiczenia (~80 min):
-* Zakładamy czepek
+* wczytanie danych do EEGLAB, edycja lokalizacji elektrod
-* Częstość próbkowania ustawiamy na 256Hz
+* filtrowanie górnoprzepustowe (FIR, 0,5 Hz), zmiana referencji
-* Na podstawie kodu: https://brain.fuw.edu.pl/edu/index.php/Laboratorium_EEG/Wprowadzenie_do_syg_online#.C4.86wiczenie:_Wykorzystanie_pomiaru_EMG_do_sterowania_on-line
+* obliczenie ICA na całym sygnale
-tworzymy programik, który zamiast pętli while tworzy odpowiednie pętle for aby:
+* identyfikacja komponentów zawierających rytm alfa:  topografia, widmo, przebieg czasowy
-* pięciokrotnie zarejestrować sekwencję trzech warunków eksperymentalnych:
+* usunięcie komponentów niezawierających alfy,  odtworzenie sygnału na elektrodach
-** patrz 5s na diodę z lewej (warunek l)
+* porównanie wyników między co najmniej trzema uruchomieniami ICA:  czy komponenty są powtarzalne?
-** patrz 5s na środek ekranu (warunek s)
-** patrz 5s na diodę z prawej (warunek p)
-* w czasie tego patrzenia:
-** odbieramy próbki w pakietach o długości 0.5s
-** każdy pakiet filtrujem (technika filtrowania on-line) w pasmach dookoła wybranych dwóch częstości (PASMO_LEWE / PASMO_PRAWE)
-** przefiltrowany pakiet przeliczamy na RMS
-** zbieramy sześć zbiorów wyników -  zestawy RMSów związane z każdym z waunków kalibracyjnych (LEWY/SPOCZYNEK/PRAWY) dla PASMO_LEWE i PASMO_PRAWE.
-* Normalizujemy RMSy. Dla pasma lewego obliczamy (RMS_(l/p/s) - np.mean(RMS_s)) / np.sdt(RMS_s) i analogicznie dla pasma prawego. Małe indeksy oznaczaają tu warunki. Zapamiętujemy współczynniki normalizacyjne.
-** Ogladamy rozkłady uzyskanych wielkości (znormalizowanych RMSów).
-** Uczymy klasyfikator np. regresję logistyczną (https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html) do rozpoznawania, z którym warunkiem patrzenia mamy do czynienia.
-==== Sesja online====
-* wczytujemy wyuczony model i kalibracyjne wartości np.mean(RMS_s)) i np.sdt(RMS_s) dla pasm lewego i prawego
-* Odbieramy próbki w pakietach o długości 0.5s
-** każdy pakiet filtrujem (technika filtrowania on-line) w pasmach dookoła wybranych dwóch częstości (PASMO_LEWE / PASMO_PRAWE)
-** przefiltrowany pakiet przeliczamy na RMS
-** Normalizujemy RMSy. Dla pasma lewego obliczamy (RMS_(l/p/s) - np.mean(RMS_s)) / np.sdt(RMS_s) i analogicznie dla pasma prawego. Małe indeksy oznaczaają tu warunki.
-** robimy predykcję klasyfikatora, z którym warunkiem patrzenia mamy do czynienia. Wyświetlamy na ekranie komunikat.
-* Testujemy czy powyższy schemat analizy pozwala na komunikację.
-{{hidden end}}
-<!--
-===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&Omega; i różnice pomiędzy oporami różnych elektrod nie przekraczały 20%.
-# Oklejamy kwadrat 3&times;3 elektrod na korze słuchowej z lewej strony (elektrody FT7, FC5, FC3, T7, C5, T3, TP7, CP5, CP3), 3&times;3 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 być 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 (Czebyszewa 2 rodzaju) zgodne z wyznaczonym pasmem. Zbadaj funkcje przenoszenia i odpowiedzi impulsowej.
-# Powycinaj sygnały od &minus;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 &minus;4 do &minus;2 s. 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.
--->
+;Materiały:
+Zob. [[#ICA jako filtr przestrzenny|ICA jako filtr przestrzenny]]
+i [[#ZADANIE: Wydobywanie interesujących komponentów|Wydobywanie interesujących komponentów]] poniżej.
 ==ICA jako filtr przestrzenny==
 {{hidden begin|title=Wstęp teoretyczny do ICA}}
@@ Linia 772: / Linia 976: @@
 {{hidden end}}
+=Sesja 6: ICA — artefakty i synteza metod=
+;Czas: ~150 min (jedno spotkanie)
+Pierwsza część sesji poświęcona jest praktycznemu zastosowaniu ICA
+do usuwania artefaktów z danych EEG przy użyciu automatycznych
+klasyfikatorów komponentów (ICLabel, MARA). Druga część syntetyzuje
+całość materiału z bloku filtrów przestrzennych: BSS/CSP, cosSinCSP
+i ICA są trzema odpowiedziami na ten sam ogólny problem separacji
+źródeł, różniącymi się założeniami i kryterium optymalizacji.
+;Zakres ćwiczenia — artefakty (~75 min):
+* wczytanie danych Arousal do EEGLAB, usunięcie kanałów nieEEG
+* obliczenie ICA, przegląd topografii komponentów
+* identyfikacja artefaktów ocznych i mięśniowych
+  przy użyciu wtyczek ICLabel i MARA
+* porównanie sygnału przed i po usunięciu komponentów artefaktowych
+;Zakres syntezy (~45 min):
+* tabela porównawcza CSP / cosSinCSP / ICA:
+  założenia, dane wejściowe, kryterium optymalizacji, kiedy stosować
+* związek filtrów przestrzennych z architekturami sieci neuronowych:
+  warstwa depthwise conv w ShallowConvNet i EEGNet jako
+  wyuczony odpowiednik macierzy W z CSP —
+  szczegóły w [[#CSP a sieci neuronowe|CSP a sieci neuronowe]] poniżej
+* omówienie raportów, pytania
+;Materiały:
+Zob. [[#ZADANIE: Identyfikacja artefaktów|Identyfikacja artefaktów]] poniżej.
 ===ZADANIE: Identyfikacja artefaktów ===
 Proszę pobrać dane:
@@ Linia 840: / Linia 1073: @@
 Przelicz potencjały z elektrod, w których występuję odpowiedź ASSR na montaż Hjortha i powtórz analizę opisaną powyżej.
 --->
+<!--
+STARA WERSJA:
+<!--
+=Prezentacja=
+[https://brain.fuw.edu.pl/edu/images/2/2f/BSS.pdf slajdy z prezentacji]
+=Ślepa separacja źródeł (BSS)=
+{{hidden begin|title=Wstęp teoretyczny do BSS}}
+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, czyli przez cały 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óry 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 Rayleigha):
+: <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{ 2 C_{T} w \left(w^T  C_{NT} w\right)  -2C_{NT} w \left(w^T  C_{T} w \right)}{\left(w^T  C_{NT} w\right)^2} = \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 do rozwiązania tzw. uogólnione zagadnienie własne:
+:<math>       C_{T}w  =\frac{w^T  C_{T} w}{w^T  C_{NT} w}   C_{NT} w   </math>
+We wzorze tym liczba <math>\lambda=\frac{w^T  C_{T} w}{w^T  C_{NT} w}</math> spełniająca to równanie jest uogólnioną wartością własną, wtedy <math>w</math> jest uogólnionym wektorem własnym odpowiadającym tej wartości.
+Aby znaleźć <math> \lambda</math> i <math>w</math> możemy wykorzystać w Matlabie funkcję <tt>eig</tt>. Funkcja ta rozwiązuje (również) uogólnione zagadnienia własne postaci ''Aw''=&lambda;''Bw'' dostarczając w wyniku macierz wektorów własnych (w kolumnach) oraz macierz zawierającą na przekątnej odpowiadające im wartości własne.
+{{hidden end}}
+-->
+<!--
+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 przekątnej 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>.
+-->
+<!--
+==Zastosowanie filtra CSP do detekcji potencjału P300==
+{{hidden begin|title=Eksperyment}}
+===Eksperyment===
+* Proszę zapoznać się z instrukcją: http://laboratorium-eeg.braintech.pl/rozdz10.html
+* Proszę wczytać i uruchomić (na sucho) demo Demos->EEG_P300wz
+* Wspólne omówienie konstrukcji i potencjalnych modyfikacji tego scenariusza
+====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 &bdquo;P300&rdquo;
+** &bdquo;P-Brain2013 Signal (with ID)&rdquo; jako np. &bdquo;Sygnal&rdquo;
+** &bdquo;P-Brain2013 Calibration p300&rdquo; jako &bdquo;kalibracjaP300&rdquo;
+** &bdquo;P-Brain 2013 p300&rdquo; jako &bdquo;Labirynt&rdquo;
+* 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:
+*** <tt>active_channels</tt>
+*** <tt>channel_names</tt>
+** dodajemy też linijkę: <tt>sampling_rate = 256</tt>
+* wchodzimy po kolei do katalogów: Sygnal_configs, kalibracjaP300_configs, Labirynt_configs i odnajdujemy plik switch_backup.ini. 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ę
+::: <tt>ignore_channels = DriverSaw;AmpSaw;PO7;PO8</tt>
+::: na
+::: <tt>ignore_channels = DriverSaw;AmpSaw;A1;A2</tt>
+::* oraz linię:
+::: <tt>montage_channels = PO7;PO8</tt>
+::: na
+::: <tt>montage_channels = A1;A2</tt>
+====Przeprowadzenie badania:====
+# Uruchom scenariusz &bdquo;Sygnał&rdquo;.
+# 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 &bdquo;Sygnal&rdquo;. 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 &bdquo;Sygnal&rdquo; 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: <tt>test1</tt>.
+# Uruchamiamy scenariusz &bdquo;kalibracjaP300&rdquo;. 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 &bdquo;Labirynt&rdquo;.
+# Danych z kalibracji potrzebować będziemy kilka zestawów.  Proszę powtórzyć kilkukrotnie scenariusz &bdquo;kalibracjaP300&rdquo;. Przed każdym uruchomieniem trzeba zmienić string w pliku <tt>file_id_name</tt> np. na <tt>test???</tt> gdzie <tt>???</tt> oznacza kolejne numery.
+-->
+<!--
+===SSVEP-BCI===
+W zajęciach tych przydadzą nam się informacje z:
+https://brain.fuw.edu.pl/edu/index.php/Laboratorium_EEG/Wprowadzenie_do_syg_online
+====Wstęp====
+Naszym celem będzie stworzenie prostego interfejsu wykorzystującego zjawisko SSVEP.
+Po dwóch stronach monitora zamocujemy diody. Każda będzie migać ze swoją ustaloną częstością (np. 16 i 22 Hz) - warto zadbać aby nie były to częstości powiązane ze sobą harmonicznie, a z drugiej strony aby, biorąc pod uwagę wyniki poprzedniego zadania, dawały dobra odpowiedź SSVEP.
+Eksperyment będzie miał dwie cześci: sesję kalibracje i sesję on-line. Pomiędzy tymi sesjami będziemy uczyć kalsyfikator. Na podstawie części kalibracyjnej ustalimy jakie parametry przetwarznego on-line sygnału świadczą o patrzeniu się na diodę z lewej a jakie na tą z prawej strony ekranu.
+Potem w sesji online będziemy porównywać (za pomocą predykcji klasyfikatora) rejestrowany sygnał z tymi warościami kalibracyjnymi i na tej podstawie system będzie zwracał informację o wyborze lewej lub prawej diody. To powinno umożliwić prostą komunikację na zasadzie pytanie i odpowiedź TAK/NIE.
+==== Sesja kalibracycjna====
+* Zakładamy czepek
+* Częstość próbkowania ustawiamy na 256Hz
+* Na podstawie kodu: https://brain.fuw.edu.pl/edu/index.php/Laboratorium_EEG/Wprowadzenie_do_syg_online#.C4.86wiczenie:_Wykorzystanie_pomiaru_EMG_do_sterowania_on-line
+tworzymy programik, który zamiast pętli while tworzy odpowiednie pętle for aby:
+* pięciokrotnie zarejestrować sekwencję trzech warunków eksperymentalnych:
+** patrz 5s na diodę z lewej (warunek l)
+** patrz 5s na środek ekranu (warunek s)
+** patrz 5s na diodę z prawej (warunek p)
+* w czasie tego patrzenia:
+** odbieramy próbki w pakietach o długości 0.5s
+** każdy pakiet filtrujem (technika filtrowania on-line) w pasmach dookoła wybranych dwóch częstości (PASMO_LEWE / PASMO_PRAWE)
+** przefiltrowany pakiet przeliczamy na RMS
+** zbieramy sześć zbiorów wyników -  zestawy RMSów związane z każdym z waunków kalibracyjnych (LEWY/SPOCZYNEK/PRAWY) dla PASMO_LEWE i PASMO_PRAWE.
+* Normalizujemy RMSy. Dla pasma lewego obliczamy (RMS_(l/p/s) - np.mean(RMS_s)) / np.sdt(RMS_s) i analogicznie dla pasma prawego. Małe indeksy oznaczaają tu warunki. Zapamiętujemy współczynniki normalizacyjne.
+** Ogladamy rozkłady uzyskanych wielkości (znormalizowanych RMSów).
+** Uczymy klasyfikator np. regresję logistyczną (https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html) do rozpoznawania, z którym warunkiem patrzenia mamy do czynienia.
+==== Sesja online====
+* wczytujemy wyuczony model i kalibracyjne wartości np.mean(RMS_s)) i np.sdt(RMS_s) dla pasm lewego i prawego
+* Odbieramy próbki w pakietach o długości 0.5s
+** każdy pakiet filtrujem (technika filtrowania on-line) w pasmach dookoła wybranych dwóch częstości (PASMO_LEWE / PASMO_PRAWE)
+** przefiltrowany pakiet przeliczamy na RMS
+** Normalizujemy RMSy. Dla pasma lewego obliczamy (RMS_(l/p/s) - np.mean(RMS_s)) / np.sdt(RMS_s) i analogicznie dla pasma prawego. Małe indeksy oznaczaają tu warunki.
+** robimy predykcję klasyfikatora, z którym warunkiem patrzenia mamy do czynienia. Wyświetlamy na ekranie komunikat.
+* Testujemy czy powyższy schemat analizy pozwala na komunikację.
+{{hidden end}}
+<!--
+===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&Omega; i różnice pomiędzy oporami różnych elektrod nie przekraczały 20%.
+# Oklejamy kwadrat 3&times;3 elektrod na korze słuchowej z lewej strony (elektrody FT7, FC5, FC3, T7, C5, T3, TP7, CP5, CP3), 3&times;3 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 być 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 (Czebyszewa 2 rodzaju) zgodne z wyznaczonym pasmem. Zbadaj funkcje przenoszenia i odpowiedzi impulsowej.
+# Powycinaj sygnały od &minus;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 &minus;4 do &minus;2 s. 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.
+-->
+-->

Anonimowy

Szukaj

Laboratorium EEG/CSP: Różnice pomiędzy wersjami

Aktualna wersja na dzień 11:05, 14 kwi 2026

Spis treści

Plan zajęć

Sesja 1: Ślepa separacja źródeł i CSP — wprowadzenie

Ćwiczenie symulacyjne

Zadania do samodzielnej realizacji

Szkielet kodu do uzupełnienia

Sesja 2: CSP dla P300 — dane realne

Sesja 3: CSP dla P300 — interpretacja i separacja cech

Analiza wstępna

ZADANIE: Analiza CSP

Wybór i separacja cech

Sesja 4: Filtry przestrzenne dla SSEP — cosSinCSP

Filtry przestrzenne dla SSEP

Teoria

Zadanie: implementacja funkcji cosSinCSP

Skrypt demonstracyjny i analiza wyników

Sesja 5: ICA — teoria i wydobywanie interesującego komponentu

ICA jako filtr przestrzenny

Definicja

Estymacja

Obliczenia

Możliwe zastosowania

Bibliografia

ZADANIE: Wydobywanie interesujących komponentów

Sesja 6: ICA — artefakty i synteza metod

ZADANIE: Identyfikacja artefaktów

W raporcie:

Nawigacja

Narzędzia Wiki

Narzędzia dla stron