% This script loads in the gamma band and speech powers for all the eight
% subjects, computes the -log p values for correlating the gamma band power
% and speech power for all the eight subjects, and plot on a generic brain
% model (stored in ALL_brain.mat) using the mapped electrode locations
% (stored in file ALL_electrodes_linear.mat).

clear all; close all; clc;

% Set directory to the present working directory.
dir=pwd;

% Define the sampling rate for the gamma band power and speech power
% signals. The data corresponds to a sampling rate of 100 Hz.
samprate=100;

% Define the time latencies (in seconds) to compute -log p values for and 
% to plot cortical activations for.
t=[-0.44,-0.24,-0.12,0,0.08,0.22,0.44];

% For each time latency, shift the gamma signal by the time lag/lead, and
% then compute the p-value using the function scampedpvalue.m.
pvlue=[];

% List of subjects 
subdir={'A','B','C','D','E','F','G','H'};

for i=1:size(t,2)                                                              
          pval=[];rval=[];
          fprintf(1,['Time:' num2str(t(i)*1000) 'ms \n']);
          data=[];
for sub=1:8
    fprintf(1,['Computing p value for subject' subdir{sub} '\n']);
    gamma=csvread([dir '/Data/Sub' subdir{sub} '_gammabandpower_overt.csv']);
    speech=csvread([dir '/Data/Sub' subdir{sub} '_speechpower.csv']);
    
    % If studying overt and covert speech, uncomment the following lines of
    % code:
    % gamma_covert=csvread([dir '/Data/Sub' num2str(sub) '_gammabandpower_overt.csv']);
    % s=min(size(gamma,1),size(gamma_covert,1));
    % gamma=gamma(1:s,:);gamma_covert=gamma_covert(1:s,:);
    
    %Compute the number of samples that the gamma signal needs to delayed
    %or advanced by, based on the time latency.
    numsamples=round(t(i)*samprate);
    
    %Shift the gamma signal by the number of samples based on a positive,
    %zero, and negative time latency.
    if (t(i)>0)
            gammanew=gamma(numsamples+1:end,:);
            gammanew=cat(1,gammanew,zeros(numsamples,size(gamma,2)));
     end
     if (t(i)==0)
            gammanew=gamma;
     end
     if (t(i)<0)
			numsamples=abs(numsamples);
            gammanew=zeros(size(gamma,1),size(gamma,2));
            gammanew(numsamples+1:end,:)=abs(gamma(1:end-numsamples,:));
     end
     
     % Call the function scrambledpvalue.m that returns the -log p values 
     % from correlating the shifted gamma signal and speech.
     p=scrambledpvalue(gammanew,speech);
     
     % If studying overt and covert speech, uncomment the following line of
     % code:
     % p=scrambledpvalue(gammanew,gamma_covert);;
     pval=cat(1,pval,p);
end
      pvlue=cat(2,pvlue,pval);
end

pvlue(isnan(pvlue))=0;

% mn and mx signify the lower and upper limits, respectively, of the 
% colorbar corresponding to the cortical activation plot.
mn=10;
mx=75;

for i=1:7
    subplot(2,4,i)

    % Call the function tickerplots.m that will plot the -log p values found
    % above on a cortical brain model, by mapping individual subject electrodes
    % to a generic brain model.
    tickerplots(pvlue(:,i),mn,mx);
    title(['Time:' num2str(t(i)*1000) 'ms']);
end
set(gcf,'color',[1 1 1]);

