function [ evoked_boot, std_dev, evoked_orig, epoch ] = bootstrap_epochs( epoch, n_boot, savepath)
%Bootstrap_epochs Bootstraps epoched data
%   Bootstrappes epoched data and returns n_boot evoked data sets plus
%   standard deviation versus time

switch nargin
    case 2
        b_save = 0;
        savepath = '';
    case 3
        b_save = 1;
end

% Average
evoked_orig = cell(1,size(epoch{1},2));
for pos = 1:size(epoch,2)
    evoked_orig{pos} = mean(epoch{pos},1);
end    

% Bootstrap
evoked_boot = cell(1,size(epoch,2));
std_dev = cell(1,size(epoch,2));
nboot = n_boot;
for pos = 1:size(epoch,2)
    [bootstat,bootsam] = bootstrp(nboot,[],epoch{pos});
    for it = 1:nboot
        evoked_boot{pos}(it,:) = mean(epoch{pos}(bootsam(:,it),:));
    end
     std_dev{pos} = std(evoked_boot{pos},0,1);
end

% Opt: Save
if b_save == 1
    save(savepath,'evoked_boot', 'std_dev', 'evoked_orig');
end 

end

