clear 
close all

% load and visualize headmodel, sourcemodel and electrodes
load ('/Volumes/DATA2/ASSR/templates/bem152.mat')
load ('/Volumes/DATA2/ASSR/templates/lf.mat')
load ('/Volumes/DATA2/ASSR/templates/elec_aligned.mat')

figure;
ft_plot_sens(elec_aligned, 'style', 'or','elecsize',10,'elecshape','point');
ft_plot_mesh(lf.pos(lf.inside,:));
ft_plot_mesh(bem152.bnd(3), 'edgecolor', 'none', 'facecolor', 'r')
ft_plot_mesh(bem152.bnd(2), 'edgecolor', 'none', 'facecolor', 'g')
ft_plot_mesh(bem152.bnd(1), 'edgecolor', 'none', 'facecolor', 'b')
alpha 0.3
view(132, 14)
%%
% Elapsed time is 317.886354 seconds.
subjects = {'45048'	'47723'	'49449'	'48288'	'45670'...
            '49353'	'47522'	'44340'	'45880'	'45611'...
            '46658'	'48993'	'45859'	'45644'	'47790'...
            '45549'	'48411'	'47371'	'48135'	'46637'...
            '45631'	'47498'	'44296'	'49121'	'46658'...
            '49152'	'44212'	'49004'	'45249'};

for i=1:length(subjects)
    datapath = strcat('/Volumes/DATA2/ASSR/data/',subjects{i});
    cd (datapath)
    load data

    %%
    cfg = [];
    cfg.channel = {'all','-E30','-Cz','-E192'};
    data = ft_selectdata(cfg,data);
    %% filter data
    cfg = [];
    cfg.bpfilter = 'yes';
    cfg.bpfilttype = 'firws';
    cfg.bpfreq = [1 40];
    datafilt = ft_preprocessing(cfg,data);
    %% compute covariance
    cfg = [];
    cfg.covariance='yes';
    cfg.vartrllength = 2;
    cfg.covariancewindow = [-.5 .5];
    avgall = ft_timelockanalysis(cfg,datafilt);
    %% perform source reconstruction using the lcmv method
    cfg                    = [];
    cfg.channel            = data.label;
    cfg.elec               = elec_aligned;
    cfg.method             = 'lcmv';
    cfg.grid               = lf;
    cfg.headmodel                = bem152;
    cfg.lcmv.keepfilter    = 'yes';
    cfg.lcmv.fixedori      = 'yes';
    cfg.lcmv.projectnoise  = 'yes';
    cfg.lcmv.weightnorm    = 'nai';
    cfg.lcmv.lambda        = '5%';
    sourceall        = ft_sourceanalysis(cfg, avgall);
    %% average dipole moments within the N1 latency
    ind =find(data.time{1}>=0.17 & data.time{1}<=0.23);
    tmpmom=sourceall.avg.mom(sourceall.inside);
    mom=sourceall.avg.pow(sourceall.inside);
    for ii=1:length(tmpmom)
        mom(ii)=mean(abs(tmpmom{ii}(ind)));
    end;
    %% insert N1 averaged power and save to disc
    source = sourceall;
    source.avg.pow(source.inside)=abs(mom);
    save source source
end;

%% load source reconstructed data
clear
close all
subjects = {'45048'	'47723'	'49449'	'48288'	'45670'...
            '49353'	'47522'	'44340'	'45880'	'45611'...
            '46658'	'48993'	'45859'	'45644'	'47790'...
            '45549'	'48411'	'47371'	'48135'	'46637'...
            '45631'	'47498'	'44296'	'49121'	'46658'...
            '49152'	'44212'	'49004'	'45249'};

for i=1:length(subjects)
    datapath = strcat('/Volumes/DATA2/ASSR/data/',subjects{i});
    cd (datapath)
    load source
    source.cfg = [];
    all{i}=source;
end;
%% grand average 
cfg=[];
cfg.parameter = 'avg.pow';
gaall = ft_sourcegrandaverage(cfg,all{:});
%%
load ('/Volumes/DATA2/ASSR/templates/mri152.mat')
cfg=[];
cfg.parameter='pow';
sourceNAIint=ft_sourceinterpolate(cfg, gaall, mri152);
%%
sourceNAIint.mask =( sourceNAIint.pow > max(sourceNAIint.pow(:))*.95 |...
                     sourceNAIint.pow < min(sourceNAIint.pow(:))*.9);
cfg               = [];
cfg.method        = 'ortho';
cfg.funparameter  = 'pow';
cfg.maskparameter = 'mask';
cfg.funcolormap = 'parula';
ft_sourceplot(cfg,sourceNAIint);