Alex's Powcorr code

%This script assumes preprocessing has already been done in EEGLab.

%%
close all;
tic

%%
top = 'O:\Labmembers\Alex';    %Base folder pathing
addpath ('C:\Users\xxvi\Documents\EEG Lab');  %Add EEGlab
startup; %Function from the EEGLab folder adds field trip and runs ft_defaults
loadspec=fullfile(top,'Sample Files'); %Concatenate path to data
savespec =fullfile(top,'Sample Files');

%%Read & convert data
subj  = {'0066'};
event = {'1ppi'};
for s = 1:length(subj)

  for e = 1:length(event)

      %Read data: basic structure
      [ALLEEG,EEG,CURRENTSET,ALLCOM] = eeglab;
      %Load set
      EEG = pop_loadset( 'filename', ['NAPLS_01_',subj{s}, '_',event{e},'_mr.set'],'filepath', loadspec);
      %Store in EEG array. Why add 0 if auto indexes into next space?
      [ALLEEG, EEG, CURRENTSET] = eeg_store( ALLEEG, EEG, 0 );
      EEG = eeg_checkset( EEG ); %Checks structure consistency

      %convert to FLDTRP
      E2F = eeglab2fieldtrip( EEG, 'preprocessing') ;
      %Rewrite to savespec directory
      fn1=fullfile(savespec, ['NAPLS_01_',subj{s}, '_',event{e},'_mr.mat']);
      save(fn1,'E2F');

  end
end

%%
%Power analysis via wavelet method, multiplication in frequency domain
for s = 1:length(subj)

  for e = 1:length(event)

    load( fullfile(savespec, ['NAPLS_01_',subj{s}, '_',event{e},'_mr.mat']))

    cfg = [];
        cfg.method = 'wavelet';
        cfg.foi  = 4:0.5:8; %Delta: 1-4, Theta: 4-8, Alpha: 8-12, Beta: 13-30, Gamma: 30-50
        cfg.toi = -.5:.01:2.5;
        %Use thin wavelets to separate time into subwindows for network analysis/time. Try 3,5
        cfg.width = 3./cfg.foi;
        cfg.output = 'powandcsd'; %Return power and cross-spectra
        cfg.channel = 'all';
        cfg.channelcmb ={'all' 'all'}; %Channel pairs
        cfg.trials ='all';
        cfg.keeptrials = 'yes'; %'yes' returns individual trials
        cfg.pad = 'nextpow2'; %Resolution
    Freq = ft_freqanalysis(cfg,E2F);

    %Freq.powspctrm = squeeze(Freq.powspctrm); %Eliminate singleton dimension(why am I getting that anyways?)
    %Freq.crsspctrm = squeeze(Freq.crsspctrm); %JK this doesn't help

    %Change filename for different frequencies
    fn2 = fullfile(savespec, ['NAPLS_01_',subj{s}, '_',event{e},'_mr','theta.mat']);
    save(fn2,'Freq');

  end
end

%%
%Power correlations between EEG channels
for s = 1:length(subj)

  for e = 1:length(event)

     %Change filename for different frequencies
     load (fullfile(savespec, ['NAPLS_01_',subj{s}, '_',event{e},'_mr','theta.mat']));

     cfg = [];
        cfg.method = 'powcorr';
        fpg.output = 'powcorr';
        cfg.channel = 'all'; %All channels
        cfg.channelcmb ={'all' 'all'}; %All channel pairs
     powcorr = ft_connectivityanalysis(cfg,Freq);

     %Change filename for different frequencies
     fn3=fullfile(savespec, ['NAPLS_01_',subj{s}, '_',event{e},'_mr','thetapowcorr.mat']);
     save(fn3, 'powcorr');

  end
end

%%
%Phase-locking value from Freq data
for s = 1:length(subj)

  for e = 1:length(event)

     %Change filename for different frequencies
     load (fullfile(savespec, ['NAPLS_01_',subj{s}, '_',event{e},'_mr','theta.mat']));

     cfg = [];
        cfg.method = 'plv';
        cfg.channel = 'all'; %All channels
        cfg.channelcmb ={'all' 'all'}; %All channel pairs
        plv = ft_connectivityanalysis(cfg,Freq);

     %Change filename for different frequencies
     fn4 = fullfile(savespec, ['NAPLS_01_',subj{s}, '_',event{e},'_mr','thetaPLV.mat']);
     save(fn4, 'plv');

  end
end

%%
%Spectral coherence from Freq data
for s = 1:length(subj)

  for e = 1:length(event)

     %Change filename for different frequencies
     load (fullfile(savespec, ['NAPLS_01_',subj{s}, '_',event{e},'_mr','theta.mat']));

     cfg = [];
        cfg.method = 'coh';
        cfg.channel = 'all'; %All channels
        cfg.channelcmb ={'all' 'all'}; %All channel pairs
     coh = ft_connectivityanalysis(cfg,Freq);

     fn5 = fullfile(savespec, ['NAPLS_01_',subj{s}, '_',event{e},'_mr','thetacoh.mat']);
     save(fn5, 'coh');

  end
end

%%
%Plot power correlations or phase value
for s = 1:length(subj)

  for e = 1:length(event)

     %Change filename for different frequencies
     fullfile(savespec, ['NAPLS_01_',subj{s}, '_',event{e},'_mr','thetapowcorr.mat']);

     cfg = [];
        %cfg.parameter = 'plvspctrm';
        cfg.parameter = 'powcorrspctrm';
        cfg.xlim = 'maxmin';
        cfg.zlim = 'maxmin';
        cfg.channel = 'all';
     figure;
     ft_connectivityplot(cfg,powcorr);

  end
end