# Advanced NI Data Analysis Term Projects
> 2. Preprocess the Go/NoGo EEG data and write a program using wavelet or filter-
Hilbert transform to find the instantaneous phase of different channels. Then, try to compute the difference between channels for finding the phase-based connectivity by computing phase vector differences and testing the non-uniformity of the differences so as to test (1) any pair(s) of channels manifesting significant connectivity differences across different proportions of Go:NoGo trials (8:2, 5:5, and 2:8) (120 points max.).

## Result

sig_pairs = {
    'O2', 'VEOG';
    'OZ', 'VEOG';
    'PZ', 'VEOG';
    'C4', 'CZ';
    'TP8', 'F3';
    'P7', 'FP2';
    'CPZ', 'T7';
    'FC4', 'F3';
    'F4', 'F7';
    'T7', 'HEOG';
    'FP2', 'VEOG'
};

These pairs show significant variation in PLV (phase-locking value) between Go and NoGo conditions across task probability contexts (2:8, 5:5, 8:2), with effect sizes (Kendall’s W) above 0.3, indicating at least moderate effect strength.

1. Fronto-Ocular Interactions (VEOG Involvement)
O2–VEOG, OZ–VEOG, PZ–VEOG, FP2–VEOG
These suggest tight coupling between posterior visual/parietal areas and ocular movements.
May reflect:
Enhanced top-down visual attention during inhibitory demands.
Eye-movement suppression associated with response inhibition.
2. Central-Motor Integration
C4–CZ
Both sites are over sensorimotor cortex.
Reflects coordinated motor preparation or inhibition.
3. Fronto-Temporal Coupling
TP8–F3, CPZ–T7, FC4–F3
Involves executive and auditory/attentional integration areas.
TP8 and T7 (temporal-parietal junction) often support attention reorienting; F3 and FC4 (left frontal) support inhibitory control.
4. Fronto-Parietal and Prefrontal Interactions
P7–FP2, F4–F7
Frontal and parietal regions are core to executive attention and response regulation.
F4–F7 may reflect lateralized frontal control dynamics.
FP2 (right prefrontal) is especially relevant to top-down control over attention and inhibition.
5. Lateral Fronto-Ocular Connection
T7–HEOG
T7 (left temporal) and HEOG (horizontal eye movement).
May reflect coordinated visual scanning or inhibition of gaze shifts under differing expectations.
Could also relate to speech-related motor planning under increased load.

## ANALYSIS PROCESS

- [ ] s1_28 55 82 = Go:Nogo = 2:8; 5:5; 8:2
- [ ] **Sample: SUBJ 1-1(s1_28)**

### EEG Loading
- s1_28 = Go:Nogo = 2:8


% Add EEGLAB to your MATLAB path. This line is fine.
addpath(genpath('/Users/nanxi/Documents/MATLAB/eeglab2025.0.0'));

% Launch EEGLAB GUI and load EEGLAB functions. This is also fine.
eeglab

% Load your EEG data. This line is fine.
EEG = pop_loadcnt('/Users/nanxi/Desktop/!!NYCU/。課（待更新讀與筆記）/。神經科學學程/神經資訊/advanced neuro-image /GNG_EEG/prob28/s1_28_mo.cnt', 'dataformat', 'auto');

% --- NEW CODE FOR CHANNEL LOCATIONS ---

% 1. Define the base path to your EEGLAB installation
%    Make sure this path is accurate for your setup.
eeglab_base_path = '/Users/nanxi/Documents/MATLAB/eeglab2025.0.0';

% 2. Construct the full path to the standard 10-20 system channel file
%    This file is typically found within the 'plugins/dipfit/standard_BEM/elec/' folder
chanlocs_file = fullfile(eeglab_base_path, 'plugins', 'dipfit', 'standard_BEM', 'elec', 'standard_1020.elc');

% 3. IMPORTANT: Verify the channel location file exists
if ~exist(chanlocs_file, 'file')
    error(['Channel location file not found: ', chanlocs_file, ...
           newline, 'Please check if EEGLAB is correctly installed and the path is accurate.']);
end

% 4. Load the channel locations into your EEG dataset
%    'loadfile', 'on' instructs pop_readlocs to actually load the specified file.
%    'filetype', 'autodetect' helps EEGLAB recognize the file format.
EEG = pop_chanedit(EEG, 'lookup', '/Users/nanxi/Documents/MATLAB/eeglab2025.0.0/sample_locs/Standard-10-20-Cap81.locs');


% Display structure field
whos EEG

% Key fields of interest:
EEG.srate         % Sampling rate
EEG.nbchan        % Number of channels
EEG.trials        % Number of trials (1 if continuous)
EEG.pnts          % Number of time points per trial
EEG.chanlocs      % Channel labels and coordinates
EEG.event         % Event markers

figure;
ch = 1;
plot(EEG.times, EEG.data(ch,:));  % channel 1
xlabel('Time (ms)');
ylabel('Amplitude (μV)');
title(['Channel: ' EEG.chanlocs(ch).labels]);

figure;
ch = 5;
plot(EEG.times, EEG.data(ch,:));  % channel 1
xlabel('Time (ms)');
ylabel('Amplitude (μV)');
title(['Channel: ' EEG.chanlocs(ch).labels]);

figure;
ch = 10;
plot(EEG.times, EEG.data(ch,:));  % channel 1
xlabel('Time (ms)');
ylabel('Amplitude (μV)');
title(['Channel: ' EEG.chanlocs(ch).labels]);

{EEG.chanlocs.labels}

### Preprocessing
- DOWN Sr to 250Hz
- Band Passing: pop_eegfilt (eeglab)
- ICA Label: erase the **red** ones; **red** >0.6


% downsampling
EEG = pop_resample(EEG, 250)

% bandpass
EEG = pop_eegfiltnew(EEG, 3, 20, [], 0, [], 0); % FIR bandpass, zero-phase

% ICA
EEG = pop_runica(EEG, 'extended',1,'interupt','on');
% --- NEW: Automatic Component Detection using ICLabel ---
% 1. Run ICLabel to classify components
EEG = iclabel(EEG); % This adds classification probabilities to EEG.etc.ic_classification.ICLabel
% 2. Flag components for rejection based on ICLabel probabilities
%    You set thresholds for what to reject. Common thresholds:
%    - Reject if 'Eye' probability > 0.6 (60%) AND 'Brain' probability < 0.4 (40%)
%    - Reject if 'Muscle' probability > 0.6
%    - Reject if 'Heart' probability > 0.6
%    - Reject if 'Line Noise' probability > 0.6
%
%    These thresholds are customizable. You might need to adjust them based on your data.
%    pop_icflag will flag components based on these criteria and store them in EEG.reject.gcompreject
%
%    The syntax for pop_icflag is:
%    pop_icflag(EEG, class_names, thresholds, flags_to_reject, 'defaultreject')
%    flags_to_reject: 'auto' means it will reject components flagged by ICLabel.
%    'defaultreject' means it will use ICLabel's recommended default thresholds.
%    Let's try with default settings first:
EEG = pop_icflag(EEG, [NaN NaN; % Brain: no threshold for flagging brain as artifact
                       0.6 1;  % Eye: if Eye probability > 80% (and Brain < 20% by default if not specified)
                       0.6 1;  % Muscle: if Muscle probability > 80%
                       0.6 1;  % Heart: if Heart probability > 80%
                       0.6 1;  % Line Noise: if Line Noise probability > 80%
                       0.6 1;  % Channel Noise: if Channel Noise probability > 80%
                       NaN NaN;% Other: no threshold for flagging 'Other' as artifact
                      ]);

% 3. Review flagged components (Highly Recommended)
%    Although automated, it's good practice to visually review the flagged components.
%    pop_selectcomps will now open with the automatically flagged components already highlighted (usually in red).
%    You can then manually unflag any components you believe are brain activity, or flag others that ICLabel missed.
EEG = pop_selectcomps(EEG, 1:EEG.nbchan); % Ensure to display all existing components (e.g., 1:32)
% 4. Remove the selected components
EEG = pop_subcomp(EEG, [], 0); % remove selected components (those flagged manually or by ICLabel)



figure;
ch = 10;
plot(EEG.times, EEG.data(ch,:));  % channel 1
xlabel('Time (ms)');
ylabel('Amplitude (μV)');
title(['Channel: ' EEG.chanlocs(ch).labels]);

### Event Label Checking

% Safely loop through event types and display them
fprintf('Event types:\n');
for i = 1:length(EEG.event)
    t = EEG.event(i).type;
    if isnumeric(t)
        fprintf('%d\n', t);
    elseif ischar(t)
        fprintf('%s\n', t);
    elseif isstring(t)
        fprintf('%s\n', t);
    else
        disp(t);  % fallback
    end
end

%Epoching
for i = 1:length(EEG.event)
    EEG.event(i).type = num2str(EEG.event(i).type);
end

event_types = {EEG.event.type};
unique_types = unique(event_types);
disp(unique_types);

> {'1'} {'11'} {'12'} {'13'} {'21'} {'22'} {'23'} : 1=GO; 2=NOGO

### Epoching
- Windows = 2000ms(2s)
- Base line = (-500ms,0)
- Event = 0ms


% --- NEW CODE FOR GO/NOGO PLV MATRICES ---

% 1. Identify Go and NoGo trial indices
% You need to know your event codes for Go and NoGo trials.
% Let's assume Go trials are code 11, 12, 13 and NoGo are 21, 22, 23 (based on your earlier code)

go_event_codes = ['11' '12' '13'];
nogo_event_codes = ['21' '22' '23'];

go_trial_indices = [];
nogo_trial_indices = [];

% Loop through EEG.event to find the trial type for each epoch
% Assuming your 'EEG' variable here is the one after epoching for ALL events
% and 'EEG.event' contains the event information for each *epoch*.
% If you epoched separately earlier, you might need to combine them first or iterate.
% For simplicity, we assume EEG.event here corresponds to the current epoched EEG.
for i = 1:EEG.trials
    current_epoch_event_type = EEG.event(i).type; % Get the event type for the current epoch
    
    if ismember(current_epoch_event_type, go_event_codes)
        go_trial_indices = [go_trial_indices, i];
    elseif ismember(current_epoch_event_type, nogo_event_codes)
        nogo_trial_indices = [nogo_trial_indices, i];
    end
end




% --- Epoching and Plotting ---

% IMPORTANT: Convert event types to numeric BEFORE the loop
% This should be done once on the continuous data
for i = 1:length(EEG.event)
    if ischar(EEG.event(i).type) % Check if it's a string before converting
        EEG.event(i).type = str2num(EEG.event(i).type);
    end
end

% Store the original continuous EEG dataset BEFORE epoching
% We will make a copy of this for each iteration
ORIG_EEG = EEG; 

% Create a figure per event type
for code = [11 12 13 21 22 23]
    % Create a temporary EEG dataset for each epoching operation
    % This ensures pop_epoch always starts from the continuous data
    CURRENT_EEG = EEG; 
    
    % Epoch the data for the current event code
    CURRENT_EEG = pop_epoch(CURRENT_EEG, {num2str(code)}, [-0.5 1.5]);
    
    % Remove baseline from the epoched data
    CURRENT_EEG = pop_rmbase(CURRENT_EEG, [-500 0]);
    
    % Only plot if epochs were actually created for this code
    if ~isempty(CURRENT_EEG.epoch)
        figure;
        % Ensure pop_timtopo uses the correct time range for the epoched data
        pop_timtopo(CURRENT_EEG, [CURRENT_EEG.xmin*1000 CURRENT_EEG.xmax*1000], [NaN], sprintf('Event code %d', code));
        sgtitle(sprintf('Event code %d - Topography over time', code)); % Add a main title to the figure
    else
        fprintf('No epochs found for event code %d. Skipping plot.\n', code);
    end
end

### Analysis
#### wavelet

- spectrum : a bulge is manifested on 6-8Hz -> could be theta band

pop_spectopo(EEG, 1, [EEG.xmin*1000 EEG.xmax*1000], 'EEG','freqrange', [3 20],'electrodes', 'off');

- -> theta band (3-8Hz) -> 6Hz for central freq

current_EEG = EEG_epoched;

frequency = 6; % Hz, for theta band
srate = current_EEG.srate;
time = -0.5:1/srate:1.5; 
if length(time) ~= current_EEG.pnts
    warning('Time vector length does not match EEG.pnts. Adjusting time vector.');
    time = linspace(current_EEG.xmin, current_EEG.xmax, current_EEG.pnts);
end

#### instantaneous phase
- event 1 = 32 channels / once

s = (6/(2*pi*frequency))^2;
wavelet = exp(2*1i*pi*frequency.*time) .* exp(-time.^2./(2*s)/frequency);
n_wavelet = length(wavelet)+1;
n_data = current_EEG.pnts; % Number of time points in each epoch
n_convolution = n_wavelet + n_data - 1;

% Initialize phase_data for ALL trials in EEG_epoched
phase_data = zeros(current_EEG.nbchan, n_data, current_EEG.trials);

% Loop over ALL trials (Go and NoGo) and channels to get phase data
for trial = 1:current_EEG.trials
    for chan = 1:current_EEG.nbchan
        eeg_data = squeeze(current_EEG.data(chan,:,trial));
        fft_data = fft(eeg_data, n_convolution);
        fft_wavelet = fft(wavelet, n_convolution);
        conv_result = ifft(fft_wavelet .* fft_data, n_convolution); % Removed *sqrt(s) here, it's often not in definition
        half_wave = floor(n_wavelet/2);
        conv_result = conv_result(half_wave + 1 : end - half_wave);
        phase_data(chan,:,trial) = angle(conv_result); % Instantaneous phase
    end
end

#### PLV(phase locking value)
- evaluate connectivity : 32 x 32 Matrics
- Calculate PLV for Go trials and NoGo trials

% --- Identify Go and NoGo trial indices from the COMBINED epoched dataset ---

go_event_codes = [11 12 13];
nogo_event_codes = [21 22 23];

go_trial_indices = [];
nogo_trial_indices = [];

for i = 1:current_EEG.trials
    current_epoch_event_type = current_EEG.event(i).type; 
    if ismember(current_epoch_event_type, go_event_codes)
        go_trial_indices = [go_trial_indices, i];
    elseif ismember(current_epoch_event_type, nogo_event_codes)
        nogo_trial_indices = [nogo_trial_indices, i];
    end
end

% Extract phase data for Go and NoGo trials
phase_data_go = phase_data(:,:,go_trial_indices);
phase_data_nogo = phase_data(:,:,nogo_trial_indices);

num_go_trials = length(go_trial_indices);
num_nogo_trials = length(nogo_trial_indices);

disp(['Number of Go trials: ', num2str(num_go_trials)]);
disp(['Number of NoGo trials: ', num2str(num_nogo_trials)]);

% 2. Calculate PLV for Go trials
PLV_Go = zeros(current_EEG.nbchan, current_EEG.nbchan);
for ch1 = 1:current_EEG.nbchan
    for ch2 = 1:current_EEG.nbchan
        if ch1 == ch2
            PLV_Go(ch1, ch2) = 1; % Self-coherence is always 1
            continue;
        end
        instant_phase_diff_rad_go = squeeze(phase_data_go(ch1,:,:) - phase_data_go(ch2,:,:));
        complex_phase_vectors_go = exp(1i * instant_phase_diff_rad_go);
        all_complex_phase_vectors_go = complex_phase_vectors_go(:);
        PLV_Go(ch1,ch2) = abs(mean(all_complex_phase_vectors_go));
    end
end

figure;
imagesc(PLV_Go); colorbar; title('Go PLV'); xlabel('Channels'); ylabel('Channels');


% 3. Calculate PLV for NoGo trials
PLV_NoGo = zeros(current_EEG.nbchan, current_EEG.nbchan);
for ch1 = 1:current_EEG.nbchan
    for ch2 = 1:current_EEG.nbchan
        if ch1 == ch2
            PLV_NoGo(ch1, ch2) = 1; % Self-coherence is always 1
            continue;
        end
        instant_phase_diff_rad_nogo = squeeze(phase_data_nogo(ch1,:,:) - phase_data_nogo(ch2,:,:));
        complex_phase_vectors_nogo = exp(1i * instant_phase_diff_rad_nogo);
        all_complex_phase_vectors_nogo = complex_phase_vectors_nogo(:);
        PLV_NoGo(ch1,ch2) = abs(mean(all_complex_phase_vectors_nogo));
    end
end

figure;
imagesc(PLV_NoGo); colorbar; title('NoGo PLV'); xlabel('Channels'); ylabel('Channels');

### Statistic
- no significant under fdr
- without fdr, pairs of channels with p<0.05 and W>0.3(effect size)

% Statistical test
% Parameters
nSubjects = 5;
nChannels = size(plv_diff_2_8, 1);
alpha = 0.05;

% Store results
p_friedman = NaN(nChannels, nChannels);
effect_size_W = NaN(nChannels, nChannels);
posthoc_p = NaN(nChannels, nChannels, 3);  % [5v5, 2v8, 5v8]

% Loop through upper triangle (channel pairs)
for ch1 = 1:nChannels
    for ch2 = ch1+1:nChannels
        
        % Extract PLV diffs for each condition
        vals_2_8 = squeeze(plv_diff_2_8(ch1, ch2, :));
        vals_5_5 = squeeze(plv_diff_5_5(ch1, ch2, :));
        vals_8_2 = squeeze(plv_diff_8_2(ch1, ch2, :));
        
        % Combine into subjects × conditions matrix
        data = [vals_2_8, vals_5_5, vals_8_2];
        
        % Friedman test
        p = friedman(data, 1, 'off');
        p_friedman(ch1, ch2) = p;
        
        % Effect size: Kendall's W
        ranked_data = tiedrank(data);  % rank each row (subject) across conditions
        SS_total = sum((mean(ranked_data, 2) - mean(ranked_data(:))).^2);
        W = 12 * SS_total / (nSubjects^2 * (3^2 - 1));
        effect_size_W(ch1, ch2) = W;

        % Post-hoc Wilcoxon signed-rank tests
        [~, p1] = signrank(vals_2_8, vals_5_5);  % 2:8 vs 5:5
        [~, p2] = signrank(vals_2_8, vals_8_2);  % 2:8 vs 8:2
        [~, p3] = signrank(vals_5_5, vals_8_2);  % 5:5 vs 8:2
        posthoc_p(ch1, ch2, :) = [p1, p2, p3];
    end
end

% Flatten p-values for correction
all_p = p_friedman(~isnan(p_friedman));
[~, ~, adj_p_friedman] = fdr_bh(all_p, alpha, 'pdep', 'yes');

% Reinsert into matrix form
p_friedman_fdr = NaN(nChannels, nChannels);
mask = ~isnan(p_friedman);
p_friedman_fdr(mask) = adj_p_friedman;


% Optional: also FDR-correct post-hoc
posthoc_p_flat = posthoc_p(~isnan(posthoc_p));
[~, ~, adj_posthoc_p] = fdr_bh(posthoc_p_flat, alpha, 'pdep', 'yes');
posthoc_p_adj = NaN(size(posthoc_p));
posthoc_p_adj(~isnan(posthoc_p)) = adj_posthoc_p;

% Save results
save('PLV_stats_results.mat', 'p_friedman_fdr', 'effect_size_W', 'posthoc_p_adj');


% Thresholds
alpha = 0.05;
effect_size_thresh = 0.3;

% Initialize results array
results = [];

% Loop through upper triangle of matrices
for ch1 = 1:nChannels
    for ch2 = ch1+1:nChannels
        if p_friedman(ch1, ch2) < alpha && effect_size_W(ch1, ch2) > effect_size_thresh %it should be p_friedman_fdr, however, there will be no  significant results. 
            % Extract post-hoc p-values
            p1 = posthoc_p_adj(ch1, ch2, 1);
            p2 = posthoc_p_adj(ch1, ch2, 2);
            p3 = posthoc_p_adj(ch1, ch2, 3);

            % Append to results
            results = [results; ch1, ch2, ...
                                 p_friedman(ch1, ch2), ...
                                 effect_size_W(ch1, ch2), ...
                                 p1, p2, p3];
        end
    end
end

% Check if results are empty
if ~isempty(results)
    % Create readable table
    results_table = array2table(results, ...
        'VariableNames', {'Chan1', 'Chan2', ...
                          'Friedman_p', 'EffectSize_W', ...
                          'Posthoc_p_28_vs_55', ...
                          'Posthoc_p_28_vs_82', ...
                          'Posthoc_p_55_vs_82'});

    disp(results_table);
    writetable(results_table, 'Significant_Channel_Pairs_no_fdr.csv');
else
    disp('⚠️ No significant channel pairs found with current thresholds.');
end

chanlocs = EEG.chanlocs;

figure;
topoplot([], chanlocs, 'style', 'blank', 'electrodes', 'labelpoint');

for i = 1:size(results, 1)
    ch1 = results(i,1);
    ch2 = results(i,2);
    disp(chanlocs(ch1).labels)
    disp(chanlocs(ch2).labels)
end

## Result
see top