% step 3_4 whole grid path_template = '/net/server/mnt/Archive/aut_gamma/matlab_analysys/Zvuk/main_data/'; load('/net/server/mnt/programs/olgas/fieldtrip-20160119/template/sourcemodel/standard_sourcemodel3d10mm.mat'); template_grid = sourcemodel; clear sourcemodel; Subj_ID = { % 'K0048' % % 'R0047' % 'R0055' % 'R0056' % 'R0057' % 'R0058' % 'R0059' % 'R0060' % 'R0064' 'R0066' 'R0043' 'R0037' 'RR010' 'RR013' }; %version 1 - all individual %version 1 - all individual for Subj_i =1:length(Subj_ID); % try filename_dataLE = [path_template Subj_ID{Subj_i} '_dataLE'] temp = load(filename_dataLE); data = temp.dataLE4save; % filename_grid = [path_template Subj_ID{Subj_i} '_Tgrid'] % temp = load(filename_grid); % gridT = temp.grid4save; filename_vol = [path_template Subj_ID{Subj_i} '_vol'] temp = load(filename_vol); vol = temp.vol4save; % % % % filename_grid = [path_template Subj_ID{Subj_i} '_grid'] % % temp = load(filename_grid); % % grid = temp.grid4save; % % filename_mri = [path_template Subj_ID{Subj_i} '_mri'] temp = load(filename_mri); mri = temp.mri4save; % % % % %preprocessing cfg = []; cfg.bpfilter = 'yes'; % apply lowpass filter cfg.bpfreq = [35 45]; data_LE_bpf3545 = ft_preprocessing(cfg, data); %% 3) compute the beamformer filter for this region cfg = []; cfg.covariance = 'yes'; cfg.channel = 'MEG';% check separately fo grad and mag cfg.vartrllength = 2; cfg.covariancewindow = 'all'; timelock = ft_timelockanalysis(cfg, data_LE_bpf3545); % create the subject specific grid, using the template grid that has just been created cfg = []; cfg.grid.warpmni = 'yes'; cfg.grid.template = template_grid; cfg.grid.nonlinear = 'yes'; % use non-linear normalization cfg.mri = mri;%TD11_mri{1}; %cfg.grid.unit ='mm'; %gridWtemplate_auto = ft_prepare_sourcemodel(cfg); gridT = ft_prepare_sourcemodel(cfg); filename_grid = [path_template Subj_ID{Subj_i} '_Tgrid'] grid4save = gridT; save(filename_grid,'grid4save','-v7.3'); % vol_cm = ft_convert_units(vol, 'cm'); % grid_test = gridWtemplate_auto; % % make a figure of the single subject headmodel, and grid positions % figure; hold on; % ft_plot_vol(vol_cm, 'edgecolor', 'none', 'facealpha', 0.4); % ft_plot_mesh(grid_test.pos(grid_test.inside,:)); cfg = []; cfg.method = 'lcmv'; %cfg.hdmfile = 'SubjectCMC.hdm'; cfg.headmodel = vol;%_nii_neuromag; %or project grad to common headmodel? cfg.grid = gridT; cfg.keepfilter = 'yes'; cfg.lcmv.fixedori='yes'; %use optimal orientation source = ft_sourceanalysis(cfg, timelock); filename_source = [path_template Subj_ID{Subj_i} '_sourceLE_Tgrid'] source4save = source; save(filename_source,'source4save','-v7.3'); % 4) reconstruct time course of the activity in the region of interest for vox_indx = 1:length(source.avg.filter) if ~isempty(source.avg.filter{vox_indx,1}) % construct the 3-D virtual channel at the location of interest beamformer_freqDISC{vox_indx} = source.avg.filter(vox_indx); end; end; chansel = ft_channelselection('MEG', data_LE_bpf3545.label); % find the names chansel = match_str(data_LE_bpf3545.label, chansel); % find the indices sourcedata = []; % sourcedata.label = {'x', 'y', 'z'}; % sourcedata.label = {'opt'}; % sourcedata.label = {'opt'}; sourcedata.time = data_LE_bpf3545.time; for i=1:length(data_LE_bpf3545.trial) for vox_indx = 1:length(source.avg.filter) if ~isempty(source.avg.filter{vox_indx,1}) sourcedata.trial{i,vox_indx} = beamformer_freqDISC{vox_indx}{1,1} * data_LE_bpf3545.trial{i}(chansel,:); end; end; end; % %view results % cfg = []; % cfg.viewmode = 'butterfly';%'vertical'; % you can also specify 'butterfly' % ft_databrowser(cfg, timelock_source); % % cfg = []; % cfg.viewmode = 'vertical'; % you can also specify 'butterfly' % cfg.channel = 'MEG1312'; % ft_databrowser(cfg, data_LE_bpf3545); % 5) computer complex FFT, exctract ITC as at sensor-level source_temp = sourcedata; for vox_indx = 1:length(source.avg.filter) if ~isempty(source.avg.filter{vox_indx,1}) source_temp.trial = sourcedata.trial(:, vox_indx); source_temp.label = {'temp'}; cfg = []; cfg.method = 'wavelet'; cfg.output = 'fourier'; cfg.keeptrials = 'yes';% or 'no', return individual trials or average %(default = 'no') cfg.tapsmofrq = 4; %tapping smoothing parameter cfg.t_ftimwin = 0.3; cfg.foi = [37:3:43]; cfg.toi = [-0.5:0.025:0.8]; freq_LE_grid{vox_indx} = ft_freqanalysis(cfg, source_temp); F = freq_LE_grid{vox_indx}.fourierspctrm; freq_LE_grid{vox_indx}.totalpow_mean = squeeze(mean(abs(F),1)); freq_LE_grid{vox_indx}.itpc = squeeze(abs(sum(F,1) ./ sum(abs(F),1))); nTrials = size(F,1); r = randperm(nTrials); sel50 = r(1:50); F50 = freq_LE_grid{vox_indx}.fourierspctrm([sel50],:,:,:); freq_LE_grid{vox_indx}.itpc_50tr = squeeze(abs(sum(F50,1) ./ sum(abs(F50),1))); %timelock cfg = []; timelock_source{vox_indx} = ft_timelockanalysis(cfg, source_temp); cfg = []; cfg.method = 'wavelet'; cfg.output = 'fourier'; cfg.keeptrials = 'yes'; %or 'no', return individual trials or average %(default = 'no') cfg.tapsmofrq = 4; %tapping smoothing parameter cfg.t_ftimwin = 0.3; cfg.foi = [1:2:100]; cfg.toi = [-0.5:0.025:0.8]; cfg.foi = [37:3:43]; cfg.toi = [-0.5:0.025:0.8]; freq_LEar_evoked = ft_freqanalysis(cfg,timelock_source{vox_indx}); F = freq_LEar_evoked.fourierspctrm; freq_LE_grid{vox_indx}.evokedpow = squeeze(sum(abs(F),1)); end; end; filename_freqLEgrid = [path_template Subj_ID{Subj_i} '_freqLEgrid_Tgrid'] for vxl = 1:length (freq_LE_grid) freqLEgrid{vxl} = freq_LE_grid{vxl}; end; save(filename_freqLEgrid,'freqLEgrid','-v7.3'); filename_ERPsourceLE = [path_template Subj_ID{Subj_i} '_ERPsourceLE_Tgrid'] for vxl = 1:length (freq_LE_grid) ERPsourceLE{vxl} = timelock_source{vxl}; end; save(filename_ERPsourceLE,'ERPsourceLE','-v7.3'); %catch %end; end; % %grid % for i = 1:size(source2_K0012.avg.filter,1) % temp1(i,2)= isempty(source2_K0012.avg.filter{i}); % end; % % % source2.pos = template_grid.pos; % source2_K0012.pos = template_grid.pos; % % cfg = []; % source_EvPow_gravr = ft_sourcegrandaverage(cfg, source2_K0012, source2); % % % source_mom = source; % time = source.time; % time_ms = [0.2 0.5]; % t1 = find(time>time_ms(1)); % t2 = find(time>time_ms(2)); % % for vox_indx = 1:length(source.avg.filter) % if ~isempty(source.avg.filter{vox_indx,1}) % %assign itpc +avg over time % source_mom.avg.pow(vox_indx) = squeeze(mean(source.avg.mom{vox_indx}(t1(1):t2(1)))); % end; % end; % %test4plot % cfg = []; % cfg.method = 'slice'%'surface';%%'ortho'%'slice';%'ortho'; % % cfg.atlas = '/net/server/mnt/programs/olgas/fieldtrip-20160119/template/atlas/aal/ROI_MNI_V4.nii'; % cfg.funparameter = 'avg.pow';%'itpc'; % % cfg.maskparameter = cfg.funparameter; % %cfg.funcolorlim = [0.0 40.2]; % % cfg.opacitylim = [0.0 1.2]; % % cfg.opacitymap = 'rampup'; % % cfg.title = 'ITCP'; % ft_sourceplot(cfg, source_mom); %