clear all dbstop if error % list the subjects: allsubj = {'P1', 'P3', 'P11', 'P12', 'P14', 'P15', 'P17', 'P19', 'P20', 'P21', 'P22', 'P26', 'P28', 'P29', 'P8', 'P9', 'P18', 'P27'}; blockConds = {'H', 'L', 'I'}; % blocked conditions use_sbj = [1:length(allsubj)] sessions = {[2 2 3 ], [2 2 3], [2 2 3], [5 5 4] , [3 3 2], [2 2 3], [2 2 3], [3 3 3], [2 2 3], [3 3 2], [2 2 4], [3 3 2], [2 2 3], [3 3 2], [3 3 2] ,[2 2 3], [3 3 2] ,[2 2 3]}; % list the block numbers for each block type of each subject. In this cell % array, make one row per subject. Hicoh, Locoh, Interleaved blocks = {[1:6] [1:6] [1:12]; [1:6] [1:6] [1:10]; [1:5] [1:5] [1:10]; [1:6] [1:6] [1:12]; [1:5] [1:5] [1:10]; [1:5] [1:5] [1:10]; [1:5] [1:5] [1:10]; [1:5] [1:5] [1:10]; ... [1:5] [1:5] [1:10]; [1:5] [1:5] [1:10]; [1:5] [1:5] [1:10]; [1:5] [1:5] [1:10]; [1:5] [1:5] [1:10]; [1:5] [1:5] [1:10]; [1:6] [1:6] [1:12]; [1:5] [1:5] [1:10]; [1:5] [1:5] [1:10]; [1:5] [1:5] [1:10]}; %1 - start task %4 - Fix on %8 9 - mapping on (cue) % 24 25 - dots on (left action/ right action) % 20 21 - response (left/right) % So list all the relevant triggers that we'll be concerned with (there might be others, like onset of fixation point at beginning of trial etc): relevant_triggers = [1 4 24 25 8 9 20 21]; % define other relevant trigger codes cuetrig = [8 9]; resptrig = [20 21]; stimtrig = [24 25]; datafolder = ''; %%Set path for Behaviour (.mat) data here eegfolder = ''; %%Set path for EEG (.bdf) data here % EEG parameters: fs = 1024; % sample rate nchan=128; % number of EEG channels next=8; % number of external channels (EOG etc) % load in a structure 'chanlocs' containing the x,y,z locations of each of the 128 scalp channels in the cap. load chanlocsBioSemi128; % note this was made by calling>> readlocs('cap128.loc') , and, since the locations for the 8 externals (129:136) are all (1,0,0), getting rid of those by calling chanlocs = chanlocs(1:128) LPF = 1; % 1 = low-pass filter the data, 0=don't. % Low pass filter-- It has a 3dB freq a bit lower than fc, at about 35 % Hz, and it has specially strong attenuation at exactly 50Hz (mains in Ireland) fc=37; % Low Pass Filter cutoff in Hz - this is what it aims for, but have to zoom in to freq to see what the -3dB frequency REALLY is % Get FIR filter weights B for Hamming-windowed sinc filter (see Semmlow 2004 signal proc book): wc = fc/fs*2*pi; B=[]; L=137; % For fs=1024 Hz, we figured out that L=137 and fc=37 Hz gives a 3dB cutoff freq of 34Hz, and great attenuation at 50 Hz for l=1:L n = l-ceil(L/2); if n==0, B(l)=wc/pi; else, B(l) = sin(wc*n)/(pi*n); end end B = B.*hamming(L)'; % figure; % freqz(B,1,10000,fs); detrend_data=1; % detrend the data (1) or not (0)? % define the epoch to extract: epoch_limits_ms = [-150 2450]; ts = round(epoch_limits_ms(1)/1000*fs):round(epoch_limits_ms(2)/1000*fs); % hence get the list of sample points relative to a given event that we need to extract from the continuous EEG t = ts*1000/fs; % hence get the timebase for the epoch (which we can plot average ERPs against) in milliseconds %define the baseline window: BLwin = [-100 0]; %was -100 -25 % Now loop through the subjects and extract the single trial ERPs and the other important informaiton about every trial: for s=use_sbj if strcmp(allsubj{s}, 'P1') blockConds = {'H', 'L', 'HL'}; %Filenames recorded differently for P1 else blockConds = {'H', 'L', 'I'}; end % make empty vectors and matrices for each thing we want to save on every single trial: modir = []; % motion direction, coh = []; % coherence. deadline = []; %Deadline config=[]; % condition, 1= high value on left, 2=high value on right corrLR=[]; % Correct side - which button would have been the correct one to press, 1=left, 2=right (cuedir) respLR=[]; % which button was actually pressed on the trial, 1=left, 2=right RT=[]; %More precise RT in seconds SOA=[]; %Time between cue and stimulus in sample points EEGrespLR=[]; % which button was actually pressed on the trial, 1=left, 2=right EEGRT=[]; % RT for each trial, in sample points ButtRT = []; %RT that was used for calculating awards, less precise, in seconds Reward = []; % Reward given to participant blocknum = [];% also keep track of block number so we can interpolate on per-block basis (cos bad electrodes usually fixed between blocks) cond = []; erp = []; % the EEG epochs. Whereas all of the above are 1-D row vectors, this will have dimensions 136 channels x epoch timepoints x trial - so the 3rd dimension should match the length of all preceding vectors % now list the eeg files to be used for each subject: filenames = []; filecond=[]; f=0; % initialize as empty a list of filenames and corresponding conditions for c=1:length(blockConds) for b=blocks{s,c} f=f+1; filenames{f} = [eegfolder '\' allsubj{s} '_' blockConds{c} num2str(b) '.bdf']; matfilenames{f} = [datafolder '\' allsubj{s} blockConds{c} num2str(sessions{s}(c)) num2str(b) '.mat']; filecond(f) = c; if strcmp(allsubj{s}, 'P1') matfilenames{f} = [datafolder '\' allsubj{s} '_' blockConds{c} num2str(b) '.mat']; elseif strcmp(allsubj{s}, 'P12') filenames{f} = [eegfolder '\' allsubj{s} blockConds{c} num2str(sessions{s}(c)) num2str(b) '.bdf']; elseif strcmp(allsubj{s}, 'P19') || strcmp(allsubj{s}, 'P20') || strcmp(allsubj{s}, 'P21') || strcmp(allsubj{s}, 'P22') || strcmp(allsubj{s}, 'P26') || strcmp(allsubj{s}, 'P28') || strcmp(allsubj{s}, 'P29') || strcmp(allsubj{s}, 'P18') || strcmp(allsubj{s}, 'P27') filenames{f} = [eegfolder '\' allsubj{s} blockConds{c} num2str(b) '.bdf']; end end if(strcmp(allsubj{s}, 'P22') & c<3) %P22 accidently had 2 HL sessions. Including all the data here p22_s3blocks = {[1:5], [1:5]}; for b=p22_s3blocks{c} f=f+1; filenames{f} = [eegfolder '\' allsubj{s} blockConds{c} '3' num2str(b) '.bdf']; matfilenames{f} = [datafolder '\' allsubj{s} blockConds{c} num2str(sessions{s}(c)+1) num2str(b) '.mat']; filecond(f) = c; end end end % Now go through these files, read them in and get their data: for f=1:length(filenames) EEG = pop_biosig(filenames{f}); % read in EEG - this is an EEGLAB function that outputs a structure 'EEG' with lots of fileds, the most important of which is EEG.data - the actual EEG data! if (strcmp(allsubj{s}, 'P21') && f<=10) || (strcmp(allsubj{s}, 'P28') && f>10)||(strcmp(allsubj{s}, 'P22') && ismember(f,[6:10 16:20])) EEG = pop_resample( EEG, fs); EEG.data = EEG.data([1:128 257:264],:); else %if (strcmp(allsubj{s}, 'P14') && f<=10) %Just in case wrong config file was used, let's do this for everyone. EEG = pop_resample( EEG, fs); end % so trigs and stimes are the exact same length; trigs has the trigger codes and stimes the SAMPLE POINTS in the % continuous EEG where those triggers happened. numev = length(EEG.event); % total number of events clear trigs stimes if strcmp(allsubj{s}, 'P8') for i=1:numev if strcmp(EEG.event(i).type(1:3), 'con') trigs(i)= str2num(EEG.event(i).type(11:end)); % trigger codes else trigs(i)= str2num(EEG.event(i).type); disp(trigs(i)) end stimes(i)=round(EEG.event(i).latency); % sample points in continuous data at which the triggers were sent end else for i=1:numev trigs(i)= str2num(EEG.event(i).type(11:end)); % trigger codes stimes(i)=round(EEG.event(i).latency); % sample points in continuous data at which the triggers were sent end end % Shave off beginning and end bits based on triggers, because sometimes experimenter is delayed starting or % stopping the recording, and the EEG might contain effects of the subject stretching out, yawning, dancing a jig goodstuff_start = max(stimes(find(ismember(trigs,relevant_triggers),1)) - fs, 0); % go back 1 second before first trigger that matters goodstuff_end = min(stimes(find(ismember(trigs,relevant_triggers),1,'last')) + 2*fs, size(EEG.data,2)); % go 2 sec beyond the last relevant trigger EEG.data(:,goodstuff_end+1:end)=[]; EEG.data(:,1:goodstuff_start)=[]; stimes = stimes-goodstuff_start; % stimes is adjusted here also of course, so its sample points correspond to the now truncated continuous EEG %%Fix electrode mistake for P19 if strcmp(allsubj{s}, 'P19') & filecond(f) ==3 temp = EEG.data([133,134],:); EEG.data([133,134],:) = EEG.data([135,136],:); EEG.data([135,136],:)=temp; end % detrend? if detrend_data, EEG.data = detrend(EEG.data')'; end % LP Filter if LPF, for q=1:nchan, EEG.data(q,:)=conv(EEG.data(q,:),B,'same'); end; end % it's an FIR filter so we can implement directly by convolution % find indices of triggers corresponding to the events we want to time-lock to (targets) targs = find(trigs==8 | trigs==9); %%get block parameters from mat file blockparams = get_block_params(matfilenames{f}); [Dotdir Config Cuedir ButtLR EventRT RewardRT Rew Coh] = get_trial_data(matfilenames{f}); if length(Dotdir)~=length(targs) if strcmp(allsubj{s}, 'P8') && f==23 Dotdir = Dotdir(2:end); Config= Config(2:end); Cuedir= Cuedir(2:end); ButtLR = ButtLR(2:end); EventRT= EventRT(2:end); RewardRT = RewardRT(2:end); Rew = Rew(2:end); Coh= Coh(2:end); elseif strcmp(allsubj{s}, 'P3') && f==22 targs = targs(3:end); else disp('problem') keyboard end end %%Hpf for emg [Bh, Ah] = butter(4,20*2/1024,'high'); % Now loop through the single trials and grow the trial-descriptor vectors and 'erp' matrix for n=1:length(targs) if stimes(targs(n))+ts(end) > size(EEG.data,2), continue; end % first, check whether the epoch relative to the current event is even contained within the bounds of the EEG data (maybe the recording was stopped mid-trial?) - if not, skip the trial ('continue') blocknum = [blocknum f]; coh = [coh Coh(n)]; % deadline = [deadline blockparams.DL]; cond = [cond filecond(f)]; modir = [modir Dotdir(n)]; config = [config Config(n)]; % 1 = high-on-left, 2=high-on-right corrLR=[corrLR Cuedir(n)]; % Correct side - which button would have been the correct one to press, 1=left, 2=right (cuedir) respLR=[respLR ButtLR(n)]; % which button was actually pressed on the trial, 1=left, 2=right RT=[RT EventRT(n)]; % RT for each trial, in seconds ButtRT = [ButtRT RewardRT(n)]; %RT that was used for calculating awards, less precise Reward = [Reward Rew(n)]; nextstimind = find(ismember(trigs,stimtrig) & stimes>stimes(targs(n)),1); % find the index of the next stimulus in trigs/stimes if ~isempty(nextstimind) SOA = [SOA stimes(nextstimind)-stimes(targs(n))]; end nextrespind = find(ismember(trigs,resptrig) & stimes>stimes(targs(n)),1); % find the index of the next response in trigs/stimes if ~isempty(nextrespind) EEGRT = [EEGRT stimes(nextrespind)-stimes(targs(n))]; EEGrespLR = [EEGrespLR trigs(nextrespind)-19]; % subtract 11 as a quick way to translate response trigger code to 1's and 2's for left/right else % if there WAS no next response, set the response parameters for this trial as 'not a number' EEGRT = [EEGRT nan]; EEGrespLR = [EEGrespLR nan]; end % Now extract the epoch ep = EEG.data(:,stimes(targs(n))+ts); %Baseline correction BLamp = mean(ep(:,find(t>BLwin(1) & t<=BLwin(2))),2); ep = ep - repmat(BLamp,[1,length(t)]); erp = cat(3,erp,ep); % now add the current epoch onto the growing 'erp' matrix by concatenation along the 3rd dimension end end % Now that we have all trials from all blocks, save everything for this subject: save([allsubj{s} '_raw'],'erp','t','modir','coh','deadline','config','corrLR','respLR','RT','EEGRT', 'ButtRT', 'EEGrespLR', 'SOA', 'Reward','blocknum','cond','filenames'); end