%% Fig 2A
load('s4.mat')
% Meta information
n_days = 1;
motor_electrode = 69;
aud_electrode = 13;
fs = 1200;

%% Motor
cdata_all = zeros;
total_trials = 0;
% Loop through all days
figure;
for a = 1:n_days
    
    day_struct = c71.day(a);
    n_run = size(day_struct.data_out, 2);    
    
    for b = 1:n_run
        
        temp_data = double(day_struct.data_out{1, b});
        
        % Load indices from specific run
        temp_stim_idxs = day_struct.ECoG.trigger_index{1, b}.stim_idxs;
        temp_cue_idxs = day_struct.ECoG.trigger_index{1, b}.cue_idxs;   
        temp_button_idxs = day_struct.ECoG.trigger_index{1, b}.button_idxs;
        temp_stim_schedule = day_struct.ECoG.stimulus_schedule{1, b};
        
        % Get total number of trials
        temp_n_trials = size(temp_stim_idxs, 2);
        
        % Get trials w/ specific cue/stimulus modalities, 
        aud_stim_trials = find(temp_stim_schedule(:, 5) == 1); 
        
        % Loop across all trials and get time segments for PSD      
        cue_segment = zeros(temp_n_trials, 2);
        bp_segment = zeros(temp_n_trials, 2);
        post_bp_segment = zeros(temp_n_trials, 2);
   
        for c = 1:temp_n_trials
            
            cue_segment(c, 1) = temp_cue_idxs(c);
            cue_segment(c, 2) = temp_cue_idxs(c) + 1200;

            bp_segment(c, 1) = temp_button_idxs(c) - 600;
            bp_segment(c, 2) = temp_button_idxs(c) + 600;
            
            post_bp_segment(c, 1) = temp_button_idxs(c) + 1200;
            post_bp_segment(c, 2) = temp_button_idxs(c) + 2400;  
            
        end
        
        % Loop across all trials to calculate spectrogram

        for f = 1:temp_n_trials
            
            cue_start_t = cue_segment(f, 1);
            cue_end_t = cue_segment(f, 2);
            bp_start_t = bp_segment(f, 1);
            bp_end_t = bp_segment(f, 2);
            post_bp_start_t = post_bp_segment(f, 1);
            post_bp_end_t = post_bp_segment(f, 2); 
            
            % Get Spectrogram
            h1 = PlotSpectrogram(temp_data(bp_start_t - 2400:bp_end_t + 2399, motor_electrode), fs);
            h = findobj(gca, 'Type', 'Image');
            cdata = h.CData;
            frq = h.YData;
            cdata_all = cdata_all + cdata;
            total_trials = total_trials + 1;
        end   
        
    end
    
end
cdata_mean = cdata_all/total_trials;
t = linspace(-2.5, 2.5, 513);
imagesc(t, frq, cdata_mean);
axis square; box off;
ylim([20, 200]); xlim([-1, 1]);
caxis([-10, 10]);
set(gca, 'TickDir', 'out'); set(gca, 'FontSize', 14)
set(colorbar, 'TickDir', 'out')
colormap('default')
xline(0, '--k', 'LineWidth', 1)
title('Exemplar PSD: Motor-Responsive Electrode')

%% Auditory
cdata_all = zeros;
total_trials = 0;
% Loop through all days
figure;
for a = 1:n_days
    
    day_struct = c71.day(a);
    n_run = size(day_struct.data_out, 2);    
    
    for b = 1:n_run
        
        temp_data = double(day_struct.data_out{1, b});
        
        % Load indices from specific run
        temp_stim_idxs = day_struct.ECoG.trigger_index{1, b}.stim_idxs;
        temp_cue_idxs = day_struct.ECoG.trigger_index{1, b}.cue_idxs;   
        temp_button_idxs = day_struct.ECoG.trigger_index{1, b}.button_idxs;
        temp_stim_schedule = day_struct.ECoG.stimulus_schedule{1, b};
        
        % Get total number of trials
        temp_n_trials = size(temp_stim_idxs, 2);
        
        % Get trials w/ specific cue/stimulus modalities, 
        aud_stim_trials = find(temp_stim_schedule(:, 5) == 1); 
        
        % Loop across all trials and get time segments for PSD      
        cue_segment = zeros(temp_n_trials, 2);
        bp_segment = zeros(temp_n_trials, 2);
        post_bp_segment = zeros(temp_n_trials, 2);
   
        for c = 1:temp_n_trials
            
            stim_segment(c, 1) = temp_stim_idxs(c);
            stim_segment(c, 2) = temp_stim_idxs(c) + 1200;

            bp_segment(c, 1) = temp_button_idxs(c) - 600;
            bp_segment(c, 2) = temp_button_idxs(c) + 600;
            
            post_bp_segment(c, 1) = temp_button_idxs(c) + 1200;
            post_bp_segment(c, 2) = temp_button_idxs(c) + 2400;  
            
        end
        
        % Loop across all trials to calculate spectrogram

        for e = 1:numel(aud_stim_trials)
            
            % Get trial # to analyze
            f = aud_stim_trials(e, 1);
            
            stim_start_t = stim_segment(f, 1);
            stim_end_t = stim_segment(f, 2);
            bp_start_t = bp_segment(f, 1);
            bp_end_t = bp_segment(f, 2);
            post_bp_start_t = post_bp_segment(f, 1);
            post_bp_end_t = post_bp_segment(f, 2); 
            
            % Get Spectrogram
            h1 = PlotSpectrogram(temp_data(stim_start_t - 2400:stim_start_t + 2399, aud_electrode), fs);
            h = findobj(gca, 'Type', 'Image');
            cdata = h.CData;
            frq = h.YData;
            cdata_all = cdata_all + cdata;
            total_trials = total_trials + 1;
        end   
        
    end
    
end
cdata_mean = cdata_all/total_trials;
t = linspace(-2, 2, 513);
figure;
imagesc(t, frq, cdata_mean);
axis square; box off;
ylim([20, 200]); xlim([-1, 1]);
caxis([-10, 10]);
colormap('default')
xline(0, '--k', 'LineWidth', 1)
title('Exemplar PSD: Auditory Stimulus-Responsive Electrode')


%% Fig. 2B
data_out = c71.day.data_out;

n_run = size(data_out, 2);
electrode_list = [68, 69];
electrode_list_aud = [9, 13];
e_1 = electrode_list(1); e_2 = electrode_list(2);
a_e_1 = electrode_list_aud(1); a_e_2 = electrode_list_aud(2);
f_1 = 70; f_i = 1; f_2 = 170;

% Meta information
ECoG = c71.day.ECoG;
all_stim_schedule = c71.day.ECoG.stimulus_schedule;

% Loop through all runs
z = 1;

% Load data from specific run
temp_data = data_out{1, z};
data_for_cohere = double(temp_data);
fs = ECoG.srate;

[b_hg, a_hg] = butter(2, [70, 170]/600, 'bandpass');
bp_data_1 = filtfilt(b_hg, a_hg, data_for_cohere(:, e_1));
bp_data_2 = filtfilt(b_hg, a_hg, data_for_cohere(:, e_2));
aud_bp_data_1 = filtfilt(b_hg, a_hg, data_for_cohere(:, a_e_1));
aud_bp_data_2 = filtfilt(b_hg, a_hg, data_for_cohere(:, a_e_2));
% Load indices from specific run
temp_stim_idxs = ECoG.trigger_index{1, z}.stim_idxs;
temp_cue_idxs = ECoG.trigger_index{1, z}.cue_idxs;   
temp_button_idxs = ECoG.trigger_index{1, z}.button_idxs;
temp_stim_schedule = all_stim_schedule{1, z};

% Get total number of trials/number of auditory trials in the
% particular run being analyzed
temp_n_trials = size(temp_stim_idxs, 2);

% Find which trials had a tactile component to the cue
tac_trials = find(temp_stim_schedule(:, 3) == 1);
av_trials = find(temp_stim_schedule(:, 3) == 0); % NOTE: these still may have had tactile CUE if "catch" trial

% Preallocating memory
cue_segment = zeros(temp_n_trials, 2);
stim_segment = zeros(temp_n_trials, 2);
bp_segment = zeros(temp_n_trials, 2);
post_bp_segment = zeros(temp_n_trials, 2);

% Loop across all trials to create time segments for coherence
% analyses
for i = 1:temp_n_trials
    cue_segment(i, 1) = temp_cue_idxs(i);
    cue_segment(i, 2) = temp_cue_idxs(i) + 1200;
    
    stim_segment(i, 1) = temp_stim_idxs(i);
    stim_segment(i, 2) = temp_stim_idxs(i) + 1200;

    bp_segment(i, 1) = temp_button_idxs(i) - 600;
    bp_segment(i, 2) = temp_button_idxs(i) + 600;

    post_bp_segment(i, 1) = temp_button_idxs(i) + 1200;
    post_bp_segment(i, 2) = temp_button_idxs(i) + 2400;        
end    

i = 9;

% All motor in this exemplar
figure; 
n_samples = post_bp_segment(i, 2) - cue_segment(i, 1);
x_vec = linspace(0,n_samples/1200,n_samples + 1)';
plot(bp_data_1(cue_segment(i, 1):post_bp_segment(i, 2)), 'b', 'LineWidth', 1);
hold on
plot(bp_data_2(cue_segment(i, 1):post_bp_segment(i, 2)), 'r', 'LineWidth', 1);
xline(bp_segment(i, 1) - cue_segment(i, 1), '--k', 'LineWidth', 2);
xline(bp_segment(i, 2) - cue_segment(i, 1), '--k', 'LineWidth', 2);
xline(cue_segment(i, 1)-cue_segment(i, 1), '--k', 'LineWidth', 1);
xline(cue_segment(i, 2)-cue_segment(i, 1), '--k', 'LineWidth', 1);
xticks(0:1200:7200); xticklabels(0:1:6);
box off; axis square;
set(gca, 'TickDir', 'out'); set(gca, 'FontSize', 14);
ylim([-25, 25])

[pxx_cue1, fxx_cue1] = pwelch(data_for_cohere(cue_segment(i, 1):cue_segment(i, 2), e_1), [], [], [70:1:170], fs);
[pxx_cue2, fxx_cue2] = pwelch(data_for_cohere(cue_segment(i, 1):cue_segment(i, 2), e_2), [], [], [70:1:170], fs);
[pxx_bp1, fxx_bp1] = pwelch(data_for_cohere(bp_segment(i, 1):bp_segment(i, 2), e_1), [], [], [70:1:170], fs);
[pxx_bp2, fxx_bp2] = pwelch(data_for_cohere(bp_segment(i, 1):bp_segment(i, 2), e_2), [], [], [70:1:170], fs);

% Combined Figure
figure;
subplot(221)
plot(bp_data_1(cue_segment(i, 1):cue_segment(i, 2)), 'b', 'LineWidth', 1); hold on; 
ylim([0, 8])
plot(bp_data_2(cue_segment(i, 1):cue_segment(i, 2)), 'r', 'LineWidth', 1);
title('Stimulus, 0.15', 'FontSize', 18); xlim([300, 450]); ylim([-10, 10]);
xticks([300, 450]); xticklabels([0, 0.125])
yticks([-10, -5, 0, 5, 10]);
axis square; box off;
set(gca, 'TickDir', 'out');set(gca, 'FontSize', 16)

subplot(222)
plot(fxx_cue1, log10(pxx_cue1), 'b', 'LineWidth', 2); hold on;
plot(fxx_cue2, log10(pxx_cue2), 'r', 'LineWidth', 2);
ylim([-4, 0]); yticks([-4, -2, 0])
xlim([70, 170]); xticks([70, 120, 170]);
axis square; box off;
set(gca, 'TickDir', 'out');set(gca, 'FontSize', 16)

subplot(223);
plot(bp_data_1(bp_segment(i, 1): bp_segment(i, 2)), 'b', 'LineWidth', 1);
hold on
plot(bp_data_2(bp_segment(i, 1):bp_segment(i, 2)), 'r', 'LineWidth', 1);
title('Button Press, 0.24', 'FontSize', 18); xlim([300, 450]); ylim([-8, 8]);
axis square; box off;
xticks([300, 450]); xticklabels([0, 0.125])
yticks([-8, -4, 0, 4, 8]);
axis square; box off;
set(gca, 'TickDir', 'out');set(gca, 'FontSize', 16)

subplot(224);
plot(fxx_bp1, log10(pxx_bp1), 'b', 'LineWidth', 2); hold on;
plot(fxx_bp2, log10(pxx_bp2), 'r', 'LineWidth', 2);
ylim([-3, 0]); yticks([-4, -2, 0])
xlim([70, 170]); xticks([70, 120, 170]);
axis square; box off;
set(gca, 'TickDir', 'out');set(gca, 'FontSize', 16)

