Study_1_preprocess.m 29.9 KB
Newer Older
Amit Goldenberg's avatar
Amit Goldenberg committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438
%% EEG ANALYSES FOR AMIT'S MOTIVATED CONTAGION EXPERIMENT %%
clear all; close all; clc; 
%% SETTINGS %% 

do = struct('import',0 , 'reimport',0, ...                                      % flow control - a structure that monitors our progress
            'ica',0, 'reica',0, ...                                             % rieca is an instruction whether to override ICA same with the others
            'epoch',0, 'reepoch',0, ... 
            'build', 0,  'analyze', 0);

% Andero paths (change)
    % dataDir = 'L:\REGULATION\Amit 2016\';                                       % directory for large datafiles
    % scriptDir = 'S:\Teadus\affect regulation\16 motivated contagion\';          % directory for scripts and other files (mirrored to gitLab)
    % eeglabDir = 'C:\Program Files\EegLab\eeglab\';                              % eeglab program directory
 
% Amit paths
    dataDir = 'G:\My Drive\research\emotional transfer\emotional_conformity_EEG_2017\Study_1_testing conformity\Results\'
    scriptDir = 'G:\My Drive\research\emotional transfer\emotional_conformity_EEG_2017\Study_1_testing conformity\Results\'  
    eeglabDir = 'C:\Program Files\MATLAB\R2017a\toolbox\eeglab14_1_1b\';  
    

% Daniel paths
    %dataDir = 'D:\emotion contagion eeg\amit_conformity_data\';
    %scriptDir = 'D:\emotion contagion eeg\amit_conformity\EEG_analysis_conformity_task\';
    %eeglabDir = 'C:\Program Files\MATLAB\R2016a\toolbox\eeglab_current\eeglab13_4_4b\';

% RELATIVE PATHS (don't change)
    addpath(scriptDir); addpath(eeglabDir);  % add the second part if it doesn't wrok            % add scripts and eeglab to MATLAB path for this session
    impdir = [dataDir, 'raw\'];                                                 % raw EEG files 
    indir = [dataDir,'import\']; mkdir(indir);                                  % EEG files converted to EEGLAB format
    icadir = [dataDir,'output\ica\']; mkdir(icadir);                            % EEG files with ICA solutions
    epochdir = [dataDir,'output\epoch_f01\']; mkdir(epochdir);                      % fully preprocessed EEG files

% DESIGN
    d.photoType  = {'emo' 'ntr'};
    d.phase  = {'indiv' 'social'};
    d.condition = {'low' 'high' 'same' 'none'};
    
eeglab;                                                                     % open eeglab

%% IMPORT %% takes files from the raw and move them ot the import folder in eeg format
while do.import==1;                                                         % flow control question 
    % Importing the EEG data, attaching events from excel file; relabeling the channael names and rereferencing. 
    pop_editoptions('option_storedisk', 0);                                 % set eeglab to hold actual data from several datasets in MATLAB (as opposed to loading them from disk only when needed which is useful at later stages ) 
    implist = dir([impdir, '*.raw']);                                       % make a fresh list of the contents of the raw files directory (input of this loop)
    %input to the loop
    %dir will give you the content of the folder. Structure. 
    inlist = dir([indir, '*.set']);                                         % make a fresh list of the contents of the eeglab-converted files (output of this loop)
    %output to the loop. we will end up with one .set file. 
    
    [a, a, trialinfo] = xlsread([scriptDir 'trialinfo.xlsx']);              % import the event file - this should be the big consolidated file   
    %gives you three output 1. all numbers. 2. all string. 3. cell which
    %holds both types of data. 
    trialinfo = trialinfo([1 1+find([trialinfo{2:end,strcmp('screen', trialinfo(1,:))}]==1)],:); % remove participants we cannot use
    % indexing trialinfo with a subset of trialinfo. 2:end - rows I want to
    % take. strcom - look for string call screen in the first row. 
    subList = unique([trialinfo{2:end,strcmp('student_id', trialinfo(1,:))}]); % build a cell with unique subject IDs
    subList = setdiff(subList,[6065275 6002886 5934980]);                                     % leave out subject with trigger problems
    
    % unique just filters our repeats 
    for s = subList;                                                        % loop through subject IDs
        if ~do.reimport; if any(strcmp([num2str(s),'.set'], {inlist(:).name})); continue; end; end  % if a given subject has already been imported, skip if required in do
        % if ~do.reimpot is =! 1 then  if - student idea and convert it to
        % .set and is included in the inlist name then keep on going and go
        % to the next test. 
        ALLEEG = []; EEG = [];                                              % erase anything in the eeglab. 
        %claers the eeg lab in case it had previous information. 
        % IMPORT
        impFiles = regexp({implist.name}, num2str(s));                      % find all datafiles that contain the selected subject...
        for i = 1:length(impFiles); if ~isempty(impFiles{i});               % ... and use the first to...
            EEG = pop_readsegegi([impdir implist(i).name]);                 % ... import them all . raw files to eeg files, where the magic happens. 
            [ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, CURRENTSET,'setname',num2str(s),'overwrite','on');  % names the dataset
            break
        end; end
    
        if isempty(EEG); continue; end                                      % skip the following if there's no EEG for the selected subject
        %checks if there is not eeg file it will move on to the next
        %subject. 
        
        % DOWN-SAMPLE 
        EEG = pop_resample( EEG, 500);                                      % downsample to 500 Hz
        [ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET); 

        % RE-REFERENCE
        EEG = pop_reref(EEG, [23 24] ,'refloc',...
        struct('labels',{'Cz'},'Y',{0},'X',{5.4492e-16},'Z',{8.8992},'sph_theta',{0},'sph_phi',{90},'sph_radius',{8.8992},'theta',{0},'radius',{0},'type',{'REF'},'ref',{''},'urchan',{[]},'datachan',{0}));
        
        % IMPORT CHANNEL LOCATIONS
        EEG=pop_chanedit(EEG, 'load',{[scriptDir 'locfile.ced'] 'filetype' 'autodetect'});
        [ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET); 
         
        % ADD EVENT INFORMATION
        EEG.subject = num2str(s);                                           % subject name 
        subinfo = trialinfo(1+find([trialinfo{2:end,strcmp('student_id',trialinfo(1,:))}]==s),:); % make a headerless matrix of trialinfo from the selected subject
        % subinfo is just a matrix for one participant
        count.one = 0; count.two = 0;                                       % initialize trial counters for the first and second block
        %creating two counters to keep track of the trial in each 
        for e = 1:length(EEG.event)                                         % loop through all EEG events
            if strcmp('PIC1', EEG.event(e).type)                            % if the event is picture onset in the first phase
                count.one = count.one + 1;                                  % move the trial counter for that phase
                row = find([subinfo{:,strcmp('order1',trialinfo(1,:))}]==count.one); % find the row in trialinfo corresponding to the given trial
                if isempty(row); continue; end                              % if the trial number is missing from behavioural file, skip what follows
                row = row(1);                                               % if there are more than one trial numbers in behavioural file, take the first
                for f = {'student_id', 'date', 'id', 'condition', 'photo', 'photoType'};  % loop through block-independent columns describing the stimulus...
                    [EEG.event(e:e+2).(f{:})] = deal(subinfo{row,strcmp(f{:},trialinfo(1,:))}); % ... and write their values to the event structure
                %the event structure included 
                end                    
                [EEG.event(e:e+2).phase] = deal(1);                         % generate a new field denoting the phase of the experiment
                for f = {'rate1','rt1', 'picot1', 'rateot1', 'order1'};     % loop through block-dependent columns describing the specific trial ...
                    [EEG.event(e:e+2).(f{:}(1:end-1))] = deal(subinfo{row,strcmp(f{:},trialinfo(1,:))}); % and write their values to event structure (to fields with block-independent names)
                end                    
            elseif strcmp('PIC3', EEG.event(e).type)                        % do the same for trials in the second block
                count.two = count.two + 1;
                row = find([subinfo{:,strcmp('order2',trialinfo(1,:))}]==count.two);
                if isempty(row); continue; end                              % if the trial number is missing from behavioural file, continue
                row = row(1);                                               % if there are more than one trial numbers in behavioural file, take the first
                for f = {'student_id', 'date', 'id', 'condition', 'photo', 'photoType', 'grate'}; 
                    [EEG.event(e:e+2).(f{:})] = deal(subinfo{row,strcmp(f{:},trialinfo(1,:))});
                end                    
                [EEG.event(e:e+2).phase] = deal(2);                                      
                for f = {'rate2','rt2', 'picot2', 'rateot2', 'order2', 'groupot2'};
                    [EEG.event(e:e+2).(f{:}(1:end-1))] = deal(subinfo{row,strcmp(f{:},trialinfo(1,:))});
                end                    
            end
            if any(strcmp({'PIC1' 'PIC3'}, EEG.event(e).type))  
                [EEG.event(e:e+2).photoType] = deal(d.photoType{EEG.event(e).photoType});
                [EEG.event(e:e+2).phase] = deal(d.phase{EEG.event(e).phase});
                if any(ismember(EEG.event(e).condition,[1:4]))
                    [EEG.event(e:e+2).condition] = deal(d.condition{EEG.event(e).condition});
                end
            end
        end        
        [ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);       % save changes to EEGLAB structure within Matlab
        EEG = pop_saveset(EEG, 'filename', EEG.subject, 'filepath', indir, 'savemode', 'twofiles'); % save changes to the import directory on disk
    end
    do.import = 0;                                                          % mark this phase as done
end

%% RUN ICA %% crates a foled in the output folder 
while do.ica == 1 
    % same is CPa but without the assumption that the variables are orthogenal 
    % class of algo that seperates the variances into max independent sources 
    inlist = dir([indir , '*.set']); icalist = dir([icadir, '*.set']); %create a list of all the .set files in thie folder                                              % update file lists
    for s = 1:length(inlist);   %we are looping through the list           
        
        if ~do.reica; if any(strmatch(inlist(s).name, {icalist(:).name})); continue; end; end
        % if don't want to overwrite (this is a just a loop to get you
        % started from the middle in cases your computer crashed) and if
        % the name alread exists then skip. 
        ALLEEG = []; EEG = []; CURRENTSET = 1; % we clean the structure
        EEG = pop_loadset('filename', inlist(s).name, 'filepath', indir);  % we load the eeg lab .set file
        [ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);% save the change to eeglab (i guess). 
        
        [ALLEEG EEG CURRENTSET] = pop_copyset(ALLEEG, 1, 2);   % this line makes a copy of the dataset for the high pass.                                                      % copy a separaate dataset for ICA training
        [ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET); 
        % WE are conducting high pass but the high pass destroys the
        % strength of the amplitue so we are just doing that in order to
        % get the ICA (which resuires high pass filtering)
        EEG = pop_eegfiltnew(EEG, [], 1, 3300, 1, [], 0);   % half-amplitude at 0.75 Hz, transition band edge at 1 Hz                                                   % filter the training data at 1 Hz high
        [ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET); 

        EEG = pop_epoch(EEG, {'PIC1', 'PIC3'}, [-1 3], 'epochinfo', 'yes');   %we're definine the segments which will be used for the ICA                                           % extract training epochs 
        [ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET); 

        [EEG ALLEEG(1).rejhist.preICA] = pop_rejspec(EEG, 1, [1:EEG.nbchan], -35, 35, 15, 30, 0, 0);   %cleaning the data of muscle noise        % remove segments with high muscle noise (15-30 Hz)
        [ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);% whenever we see a lot of power for specific point.

        EEG = pop_runica(EEG, 'extended', 1, 'interupt', 'on');     %actually running the ICA                                                % find an ICA solution
        [ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET); 
        
        [ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, CURRENTSET, 'retrieve', 1);   % write the solution of the ICA on the data.                            % write the ICA solution back to original data
        EEG = pop_editset(EEG, 'icachansind', 'ALLEEG(2).icachansind', 'icaweights', 'ALLEEG(2).icaweights', 'icasphere', 'ALLEEG(2).icasphere'); 
        [ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);     

        EEG = pop_saveset(EEG, 'filename', EEG.filename, 'filepath', icadir, 'savemode', 'twofiles');
    end
    do.ica = 0;
end

%% EPOCH AND CLEAN %% segments and throws out bad components (eye movement)
while do.epoch == 1;
    
    icalist = dir([icadir,'*.set']); epochlist = dir([epochdir,'*.set']); % refresh different folder content lists
              
    inlist = dir([indir,'*.set']);
    
    for s = 1:length(icalist); % we are looking through the ica directory 
        
        if ~do.reepoch; if any(strmatch(icalist(s).name, {epochlist(:).name})); continue; end; end    % overwrite check

        ALLEEG = []; EEG = []; CURRENTSET = 1;                                                      % clean workspace and load dataset from ICA folder
        
        EEG = pop_loadset('filename', icalist(s).name, 'filepath', icadir); 
        [ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);

        
        if isempty(find(EEG.reject.gcompreject));                                   % if the file does not have markings for ocular ICs ...
            
            EEG = pop_epoch(EEG, {'CROS'}, [0 4], 'epochinfo', 'yes');              % epoch just for visualisation 
            EEG = pop_rmbase(EEG, [800 1000]);                                                    
            [ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);                           

            pop_eegplot(EEG, 1, 1, 1); pop_eegplot(EEG, 0, 1, 1);                   % open scroll plots (channels, components)
            pop_selectcomps(EEG, [1:length(EEG.reject.gcompreject)]);               % open component plot

            uiwait; close(gcf); close(gcf);                                         % wait until the plots are closed

            gcompreject = EEG.reject.gcompreject;                                   % save the decisions to a variable

            ALLEEG = []; EEG = [];                                                  % delete the segmented file

            EEG = pop_loadset('filename', icalist(s).name, 'filepath', icadir);     % open a fresh copy (without segments)
            [ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);

            EEG.reject.gcompreject = gcompreject;                                   % write the decisions to the newly opened file

            EEG = pop_saveset(EEG, 'filename', EEG.filename, 'filepath', icadir, 'savemode', 'twofiles'); % save the update ICA file back to disk
    
        end                            
        % we may want to adopt an automatized algorithm instead of this thing
        
        EEG = pop_eegfiltnew(EEG, 1, 40, 1650, 0, [], 0);                        % filtering 1htz lower cutof and 90upper cutof (for erp 30-40)
        [ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET); 
        
        EEG = pop_epoch(EEG, {'PIC1' 'PIC3'}, [-0.2 1.5], 'epochinfo', 'yes');  % extract epochs 1 second before the picture and 3 after the picture. We may want to shortening
        EEG = pop_rmbase(EEG, [-200 0]);                                    % remove baseline - we establish 200 ms before the picture as the baselie. 
        [ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, CURRENTSET, 'overwrite', 'on');    % save 
        
        EEG = pop_subcomp(EEG, find(EEG.reject.gcompreject), 0);            % remove the IC components marked as artifacts - we are statistically smoothing whatever the ICA found. 
        [ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET); 
        
        % detect artifacts while ignoring a bunch of channels
        EEG = rejSegChan(EEG, {'VEOG1','VEOG2','HEOG1','HEOG2','O9', 'O10', 'F9', 'F10', 'Fn'}, 100);      % call an algorithm for rejecting artifactual epochs and channels (see help rejSegChan) - allows a deviation of 100 micro volts before rejecting
        %however, if there is one channel that doens't work we don't want
        %to remove all the data. We therefore check if this is a crucial
        %channel and keep the data. 
        %all of this filtering is occuring within the eeg strcuture.       
        [ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET); 
                        
        EEG = pop_saveset(EEG, 'filename', EEG.subject, 'filepath', epochdir, 'savemode', 'twofiles');     
    end
    do.epoch = 0;
end       

%% BUILD ANALYSIS STRUCTURE %% looks at the epoch folder (inside the output) and creates a 
while do.build == 1;    
    
    epochlist = dir([epochdir,'*.set']); 
    
    % RESEARCHER INPUT
    d.photoType  = {'emo' 'ntr'};
    d.phase  = {'indiv' 'social'};
    d.condition = {'low' 'high' 'same' 'none'};
    minEpoch = 1;                                                           % minimum number of epochs that needs to be retained per design ceafter artifact rejection before the subject is excluded
    roi = {};                                                               %  {'Pz' 'Cz' 'P3' 'P4'} define region(s) of interest, if applicable
    
    % INITIALIZE VARIABLES
    p = struct('times', {[]}, 'chanlocs', readlocs([scriptDir,'locfile.ced']), 'data', {{}}, 'subs', {{}});  % p is a custom structure we'll build to be used later as input to visualizaton functions
    s = 0; bkg = {};                                                        % s counts subjects; bkg is the meta-data matrix
    for fld = fields(d)'; p.labels.(fld{:}) = d.(fld{:}); end               % take factor and level labels from the d structure and add to the p 

    if ~isempty(roi)                                                        % if regions of interest have been defined, add a dummy channel location for each ROI
        chf = fields(p.chanlocs)';
        for r = 1:length(roi)
            [roiLoc, roiLoc] = intersect({p.chanlocs.labels},roi{r});
            p.chanlocs(end+1).labels = [p.chanlocs(roiLoc).labels];
            for f = 2:length(chf)
                p.chanlocs(end).(chf{f}) = mean([p.chanlocs(roiLoc).(chf{f})]);
            end
        end
    end    
        
    for i = 1:length(epochlist);                                                    % loop through files
        
        EEG = []; ALLEEG = []; data = []; n = [];                                   % initiate/clear various variables                                                      
        
        EEG = pop_loadset('filename', epochlist(i).name, 'filepath', epochdir);     % load file
        [ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, CURRENTSET);   
        
        s = s + 1;                                                                  % advance the subjects counter
        p.subs{s} = EEG.subject;                                                    % record subject code to p structure
        
        p.times = EEG.times;                                                        % add time vector for plotting
        
        bkg = getBkg(EEG, d, bkg, minEpoch, readlocs([scriptDir,'locfile.ced']));   % call a function to build meta-data matrix
        
        if bkg{end,end} == 0                                                        % if any of the design cells has less than required trials... 
            s = s - 1;                                                              % ... drop the subject from the p structure
            continue
        end
        
        epoch = epochList(EEG);                                                     % get information about each epoch in correct format
        p.epochFields = fields(epoch)';
        
        EEG = pop_interp(EEG, readlocs([scriptDir,'locfile.ced']), 'spherical');    % interpolate missing channels - if you threw out a channel the system will try to fill in the data
        [ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);        

        if ~isempty(roi)                                                            % if region of interest is specified, create a new channels by averaging the constituent channels
            for r = 1:length(roi)
                [roiLoc, roiLoc] = intersect({EEG.chanlocs.labels},roi{r});
                EEG.data(end+1,:,:) = mean(EEG.data(roiLoc,:,:),1);
            end    
        end
        
        % AVERAGE ERPs FOR EACH DESIGN CELL 
        for photoType = 1:2                                                         % loop through levels of photoType factor
            t.a = strcmp(d.photoType{photoType},{epoch.photoType});                 % find all trials for a given level (all level labels were converted to strings by epochList, so we match using strcmp)
            for phase = 1:2                                                         % loop through 2. factor levels
                t.b = strcmp(d.phase{phase},{epoch.phase});
                for condition = 1:4                                 
                    t.c = strcmp(d.condition{condition},{epoch.condition});
                    inds = find(t.a.*t.b.*t.c);                                     % find trials belonging to one design cell...
                    p.data{photoType,phase,condition}(:,:,s) = mean(EEG.data(:,:,inds),3);  % average trials in the selected cell
                    p.raw.dat{photoType,phase,condition,s} = EEG.data(:,:,inds);  % trial data
                    p.raw.inf{photoType,phase,condition,s} = squeeze(struct2cell(epoch(inds)))';
                end
            end
        end
    end
    
    save([epochdir date '-data.mat'], 'p');                                         % save the p structure to disk
    xlswrite([epochdir date '-export.xlsx'], bkg, 'bkg');                           % save the meta-data file
    
    do.build = 0;
end       

%% BUILD ANALYSIS STRUCTURE %%
while do.build == 1;    
    
    epochlist = dir([epochdir,'*.set']); 
    
    % RESEARCHER INPUT
    d.photoType  = {'emo' 'ntr'};
    d.phase  = {'indiv' 'social'};
    d.condition = {'low' 'high' 'same' 'none'};
    minEpoch = 1;                                                           % minimum number of epochs that needs to be retained per design ceafter artifact rejection before the subject is excluded

    ersp = struct('times', {[]}, 'chanlocs', readlocs([scriptDir,'locfile.ced']), 'data', {{}}, 'subs', {{}});  % p is a custom structure we'll build to be used later as input to visualizaton functions    
	ersp.frqs = {'theta' [3 7]; 'alpha' [8 12]; 'beta' [15 25]};
    ersp.params = {'cycles',[1 4],'freqs', [3 30], 'baseline',[-1250 -250],'trialbase','full','basenorm','off','freqscale','linear','scale','abs','timesout',length(epochsegment{2}(1):0.0039063:epochsegment{2}(2))/8,'padratio',1,'plotphase','off','plotersp','on','plotitc','on','verbose','off'};
	
    % INITIALIZE VARIABLES
    s = 0; bkg = {};                                                        		% s counts subjects; bkg is the meta-data matrix
    for fld = fields(d)'; p.labels.(fld{:}) = d.(fld{:}); end               		% take factor and level labels from the d structure and add to the p 
        
    for i = 1:length(epochlist);                                                    % loop through files
        
        EEG = []; ALLEEG = []; data = []; n = [];                                   % initiate/clear various variables                                                      
        
        EEG = pop_loadset('filename', epochlist(i).name, 'filepath', epochdir);     % load file
        [ALLEEG EEG CURRENTSET] = pop_newset(ALLEEG, EEG, CURRENTSET);   
        
        s = s + 1;                                                                  % advance the subjects counter
               
        bkg = getBkg(EEG, d, bkg, minEpoch, readlocs([scriptDir,'locfile.ced']));   % call a function to build meta-data matrix
        
        if bkg{end,end} == 0                                                        % if any of the design cells has less than required trials... 
            s = s - 1;                                                              % ... drop the subject from the p structure
            continue
        end
        
        epoch = epochList(EEG);                                                     % get information about each epoch in correct format
        p.epochFields = fields(epoch)';
        
        EEG = pop_interp(EEG, readlocs([scriptDir,'locfile.ced']), 'spherical');    % interpolate missing channels - if you threw out a channel the system will try to fill in the data
        [ALLEEG EEG CURRENTSET] = eeg_store(ALLEEG, EEG, CURRENTSET);        
        
		% COMPUTE ERSP
		for c = 1:length(EEG.chanlocs)-4
			figure; [a,a,a,times,freqs,a,a,tfdata] = newtimef(EEG.data(c,:,:), EEG.pnts, [EEG.xmin EEG.xmax]*1000, EEG.srate, ersp.params{:});
			power = tfdata.*conj(tfdata);                                           % convert to power values
			power = power./repmat(mean(power,2),[1 length(times) 1 1]);             % full-trial normalization; gain model divide by each trial mean value (based on Grandchamp & Delorme, 2011)
			baseline = ersp.params{strmatch('baseline', ersp.params(1:2:end))*2}; 
			[a lti] = min(abs(times-baseline(1))); [a uti] = min(abs(times-baseline(2))); 
			power = 10*log10(power./repmat(mean(mean(power(:,lti:uti,:,:),2),3),[1 length(times) size(power,3) 1]));% actual baseline correction on mean level + transform to log
			% power = power-repmat(mean(power(lfi:ufi,lbi:ubi,:),2),[1 length(times) 1]);                                 % actual baseline correction on trial 
			close(gcf);
	
			% AVERAGE ERSP FOR EACH DESIGN CELL 
			for photoType = 1:2                                                         % loop through levels of photoType factor
				t.a = strcmp(d.photoType{photoType},{epoch.photoType});                 % find all trials for a given level (all level labels were converted to strings by epochList, so we match using strcmp)
				for phase = 1:2                                                         % loop through 2. factor levels
					t.b = strcmp(d.phase{phase},{epoch.phase});
					for condition = 1:4                                 
						t.c = strcmp(d.condition{condition},{epoch.condition});
						inds = find(t.a.*t.b.*t.c);                                     % find trials belonging to one design cell...
						for f = ersp.frqs(:,1)'                                                                                 % keskmistab kohe sagedused kokku
							fval = ersp.frqs{strmatch(f{:},ersp.frqs(:,1),'exact'),2};
							if fval(1) < freqs(1); fval(1) = freqs(1); end; if fval(2) > freqs(end); fval(2) = freqs(end); end; 
							[a lfi] = min(abs(freqs-fval(1))); [a ufi] = min(abs(freqs-fval(2))); 
							ersp.(f{:}).times = times;
							ersp.(f{:}).subs{s} = EEG.subject;                                                    % record subject code to p structure
						end					
						ersp.(f{:}).data{photoType,phase,condition}(c,:,s) = squeeze(mean(mean(power(lfi:ufi,:,inds),1),3));  % average trials in the selected cell
					end
				end
			end
		end
    end
    
    save([epochdir date '-ERSP.mat'], 'p');                                         % save the p structure to disk
    
    do.buildERSP = 0;
end     

%% VISUALIZE AND EXPORT %% 
while do.analyze == 1;
        
    buildlist = dir([epochdir,'*.mat']);                                    % find .mat files from epochdir
    load([epochdir buildlist(end).name]);                                   % open the last one

    % there are 120 emo and 60 ntr pic, therefore the maximus should be 30
    % and 15, so something is off in the 

    % VIUSAL STATISTICS %
    %plotERPeffect(data-structure, factor-to-plot, cells-to-include, p-level, correctoon, time-to-plot, baseline, plot type)
    % see help plotERPeffect

    % Envelope plot
    %plotERPeffect(p,1,{1:2,1:2,1:4}, 1, 'none', [-200 1500], [-200 0], 'envp');        % plot grand average waveforms from each channel on a single plot

    % ERP plots
    %plotERPeffect(p,1,{1:2,1,[1:3]}, 0.05, 'none', [-200 1500], [-200 0], 'erp');        % main effect to picture type across conditions. 
    %plotERPeffect(p,2,{1,1:2,[1:3]}, 0.05, 'none', [-200 1500], [-200 0], 'erp');        % main effect on indiv vs. social 
    plotERPeffect(p,3,{1,2,[1:3]}, 0.05, 'none', [-200 1500], [-200 0], 'erp');      % e.g. condition main effect from the second phase with cluster-based statistics 

    % SCALP PLOTS
    %plotERPeffect(p,1,{1:2,2,[1:3]}, 0.99, 'no', [270 470], [-200 0], 'scalpDiff');     % plot the scalp distribution of the difference between the last and first level of the selected dimension 

    %plotERPeffect(p,1,{1:2,1,1:3}, 1.0, 'no', [170 260], [-200 0], 'scalpComp');      % plot the scalp distribution of the grand average  
    
    % EXPORT (see help exportERPeffect)
    %exportERPeffect([epochdir date  '-export.xlsx'], p,{1:2,1:2,1:4},...
    %    {'P3', {'Pz'}, [300 500];...
    %     'LPP', {'Pz'}, [500 800];...
    %     'SLW' , {'Pz'}, [800 1100]}, [-200 0], 'ch', NaN, 'no' ,'yes');                                 

    do.analyze = 0;
end