%%% build an input structure for data-processing in FieldTrip
load CSD_Data_Eyes_Open.mat;
u = new;

[chan, smp, epochs] = size(u);

% remove DC-offset and linear trends
for i = 1:epochs
    chan_mn = mean(u(:,:,i),2); % compute DC-offset
    dat = bsxfun(@minus,u(:,:,i),chan_mn); % remove DC-offset
    dat = detrend(dat,'linear'); % remove linear trend
    store(:,:,i) = dat;
end

data = []; % init data structure
 
data.fsample = 512;
data.label = chan_labels'; % chan labels stored in chan_labels
data.cgf=[]; % null config as data is not preprocessed in FieldTrip
data.trial = cell(1,1); for i = 1:epochs, data.trial{1,i} = store(:,:,i); end % data stored in u
data.time = cell(1,1); time = 0:1/512:1; for i = 1:epochs, data.time{1,i}=time; end % assign time index not NB as frequency domain analysis

%% no parametric spectral analysis

cfg            = [];
cfg.method     = 'mtmfft';
cfg.output     = 'fourier';
cfg.foilim     = [0 45];
%cfg.channel    = {'F3' 'F7' 'P3' 'P7'};
%cfg.channelcmb = {'F3' 'P3'; 'F3' 'P7'};
cfg.tapsmofrq  = 2; % smoothing tapering alt. Hanning
freq = ft_freqanalysis(cfg, data); % frequency and time-frequency analysis
fd = ft_freqdescriptives(cfg, freq); % provides spectrum and std error

%% compute coherence & granger causality

cfg = [];
cfg.method = 'coh';
coherence = ft_connectivityanalysis(cfg, freq); % performs connectivity analysis

cfg = [];
cfg.method = 'granger';
granger_c = ft_connectivityanalysis(cfg, freq); % performs connectivity analysis

%% plot results

cfg = [];
cfg.parameter = 'grangerspctrm';
%cfg.channel = {'F3', 'F7', 'P3', 'P7'};
figure;ft_connectivityplot(cfg, granger_c); % plots connectivity analysis

% figure
% for row=1:3
% for col=1:3
%   subplot(3,3,(row-1)*3+col);
%   plot(granger_c.freq, squeeze(granger_c.grangerspctrm(row,col,:)))
%   ylim([0 1])
% end
% end

cfg.parameter = 'cohspctrm';
figure;ft_connectivityplot(cfg, coherence); % plots connectivity analysis

% figure
% for row=1:3
% for col=1:3
%   subplot(3,3,(row-1)*3+col);
%   plot(coherence.freq, squeeze(coherence.cohspctrm(row,col,:)))
%   ylim([0 1])
% end
% end