123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159 |
- clear all; cla; clc;
- currentFile = mfilename('fullpath');
- [pathstr,~,~] = fileparts(currentFile);
- cd(fullfile(pathstr,'..'))
- rootpath = pwd;
- pn.data_eeg = fullfile(rootpath, '..', 'eegmp_preproc', 'data', 'outputs', 'eeg');
- pn.data_erp = fullfile(rootpath, 'data', 'erp');
- pn.data_erf = fullfile(rootpath, 'data', 'erf');
- pn.data = fullfile(rootpath, 'data', 'stats'); mkdir(pn.data);
- pn.tools = fullfile(rootpath, 'tools');
- addpath(fullfile(rootpath, '..', 'eegmp_preproc', 'tools', 'fieldtrip')); ft_defaults
- addpath(genpath(fullfile(pn.tools, '[MEG]PLS', 'MEGPLS_PIPELINE_v2.02b')))
- addpath(fullfile(pn.tools, 'BrewerMap'));
- addpath(fullfile(pn.tools, 'shadedErrorBar'));
- pn.figures = fullfile(rootpath, 'figures');
- %% add seed for reproducibility
- rng(0, 'twister');
- %% load erp
- for ind_id = 1:33
- id = sprintf('sub-%03d', ind_id); disp(id)
- load(fullfile(pn.data_erf, [id,'_erf.mat']));
- for ind_option = 1:numel(conds.behavior)
- if ind_id == 1
- erpgroup.behavior.(conds.behavior{ind_option}) = erf.behavior{ind_option};
- erpgroup.behavior.(conds.behavior{ind_option}) = ...
- rmfield(erpgroup.behavior.(conds.behavior{ind_option}), {'powspctrm'});
- erpgroup.behavior.(conds.behavior{ind_option}).dimord = 'sub_chan_freq_time';
- freq = erpgroup.behavior.(conds.behavior{ind_option}).freq;
- end
- idx_freq_t = freq >2 & freq <7;
- thetagroup.behavior.(conds.behavior{ind_option}).avg(ind_id,:,:,:) = squeeze(nanmean(erf.behavior{ind_option}.powspctrm(:,idx_freq_t,:),2));
- idx_freq_a = freq >8 & freq <25;
- alphagroup.behavior.(conds.behavior{ind_option}).avg(ind_id,:,:,:) = squeeze(nanmean(erf.behavior{ind_option}.powspctrm(:,idx_freq_a,:),2));
- end
- end
- time = erpgroup.behavior.hit.time;
- channels = erpgroup.behavior.hit.label;
- %% get max. channels
- load(fullfile(pn.data, 'd2_taskpls_erf.mat'),...
- 'stat', 'result', 'lvdat', 'lv_evt_list', 'num_chans', 'num_freqs', 'num_time')
- plotData.powspctrm = squeeze(nanmean(nanmean(stat.mask(:,stat.freq>8 & stat.freq<25,stat.time>1 & stat.time<2).*...
- stat.prob(:,stat.freq>8 & stat.freq<25,stat.time>1 & stat.time<2),3),2));
- [~, alphaChanMin] = sort(plotData.powspctrm, 'ascend');
- % plotData.powspctrm = squeeze(nanmean(nanmean(stat.mask(:,stat.freq>2 &stat.freq<7,stat.time>1 & stat.time<2).*...
- % stat.prob(:,stat.freq>2 & stat.freq<7,stat.time>1 & stat.time<2),3),2));
- % [~, thetaChan] = sort(plotData.powspctrm, 'descend');
- %
- plotData.powspctrm = squeeze(nanmean(nanmean(stat.mask(:,stat.freq>8 & stat.freq<25,stat.time>0.5 & stat.time<1).*...
- stat.prob(:,stat.freq>8 & stat.freq<25,stat.time>0.5 & stat.time<1),3),2));
- [~, alphaChanMax] = sort(plotData.powspctrm, 'descend');
- %% visualize for frontal theta
- mergeddata = cat(4, thetagroup.behavior.hit.avg, thetagroup.behavior.miss.avg);
- %figure; imagesc(squeeze(nanmean(thetagroup.old.new.avg-thetagroup.old.old.avg,1)))
- idx_chans = [38, 47];% thetaChan(1:3); channels(idx_chans) %[33,38,47];
- condAvg = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,1:2),2),4));
- h = figure('units','centimeters','position',[0 0 10 8]);
- cla; hold on;
- % new value = old value ? subject average + grand average
- curData = squeeze(nanmean(mergeddata(:,idx_chans,:,1),2));
- curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
- standError = nanstd(curData,1)./sqrt(size(curData,1));
- l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
- {'color', 'k','linewidth', 2}, 'patchSaturation', .1);
- curData = squeeze(nanmean(mergeddata(:,idx_chans,:,2),2));
- curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
- standError = nanstd(curData,1)./sqrt(size(curData,1));
- l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
- {'color', 'r','linewidth', 2}, 'patchSaturation', .1);
- xlabel('Time (s) from stim onset')
- legend([l1.mainLine, l2.mainLine],{'hit', 'miss'}, ...
- 'location', 'southeast'); legend('boxoff')
- xlim([-.5 1.50]); %ylim(YLim)
- ylabel({'theta power';'(log10)'});
- xlabel({'Time (s)'});
- set(findall(gcf,'-property','FontSize'),'FontSize',12)
- figureName = ['d2_erf_theta_trace'];
- saveas(h, fullfile(pn.figures, figureName), 'epsc');
- saveas(h, fullfile(pn.figures, figureName), 'png');
- %% visualize for posterior alpha
- mergeddata = cat(4, alphagroup.behavior.hit.avg, alphagroup.behavior.miss.avg);
- idx_chans = alphaChanMin(1:6); channels(idx_chans) %[33,38,47];
- condAvg = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,1:2),2),4));
- h = figure('units','centimeters','position',[0 0 10 8]);
- cla; hold on;
- % new value = old value ? subject average + grand average
- curData = squeeze(nanmean(mergeddata(:,idx_chans,:,1),2));
- curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
- standError = nanstd(curData,1)./sqrt(size(curData,1));
- l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
- {'color', 'k','linewidth', 3}, 'patchSaturation', .1);
- curData = squeeze(nanmean(mergeddata(:,idx_chans,:,2),2));
- curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
- standError = nanstd(curData,1)./sqrt(size(curData,1));
- l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
- {'color', 'r','linewidth', 3}, 'patchSaturation', .1);
- xlabel('Time (s) from stim onset')
- legend([l1.mainLine, l2.mainLine],{'hit', 'miss'}, ...
- 'location', 'north'); legend('boxoff')
- xlim([-.5 1.75]);
- ylabel({'alpha power';'(log10)'});
- xlabel({'Time (s)'});
- set(findall(gcf,'-property','FontSize'),'FontSize',15)
- figureName = ['d2_erf_latealpha_trace'];
- saveas(h, fullfile(pn.figures, figureName), 'epsc');
- saveas(h, fullfile(pn.figures, figureName), 'png');
- %% visualize for central alpha
- mergeddata = cat(4, alphagroup.behavior.hit.avg, alphagroup.behavior.miss.avg);
- idx_chans = alphaChanMax(1:3); channels(idx_chans) %[20,21,12];
- condAvg = squeeze(nanmean(nanmean(mergeddata(:,idx_chans,:,1:2),2),4));
- h = figure('units','centimeters','position',[0 0 10 8]);
- cla; hold on;
- % new value = old value ? subject average + grand average
- curData = squeeze(nanmean(mergeddata(:,idx_chans,:,1),2));
- curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
- standError = nanstd(curData,1)./sqrt(size(curData,1));
- l1 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
- {'color', 'k','linewidth', 3}, 'patchSaturation', .1);
- curData = squeeze(nanmean(mergeddata(:,idx_chans,:,2),2));
- curData = curData-condAvg+repmat(nanmean(condAvg,1),size(condAvg,1),1);
- standError = nanstd(curData,1)./sqrt(size(curData,1));
- l2 = shadedErrorBar(time,nanmean(curData,1),standError, 'lineprops', ...
- {'color', 'r','linewidth', 3}, 'patchSaturation', .1);
- xlabel('Time (s) from stim onset')
- legend([l1.mainLine, l2.mainLine],{'hit', 'miss'}, ...
- 'location', 'northeast'); legend('boxoff')
- xlim([-.5 1.75]);
- ylabel({'alpha power';'(log10)'});
- xlabel({'Time (s)'});
- set(findall(gcf,'-property','FontSize'),'FontSize',15)
- figureName = ['d2_erf_earlyalpha_trace'];
- saveas(h, fullfile(pn.figures, figureName), 'epsc');
- saveas(h, fullfile(pn.figures, figureName), 'png');
|