UEL Baby Dev Lab

runWindowedAnalyses.m
function runWindowedAnalyses

% Set paths
diskLetter = 'F';
rootPath = strcat(diskLetter, ':\FrontierPaperAnalysis\');
filePaths = {strcat(diskLetter, ':\Frontiers Analysis\Stage_6_preproc_data\')};
% Set target conditions, sbn8 = steady beats normal 8Hz
targetConditions = {'sbn8'; 'sbn2'; 'sbn4'; 'sbn6'; 'sbn10'; 'sbn12'};
% 'sbn2'; 'sbn4'; 'sbn6'; 'sbn10'; 'sbn12';
% Set participants
participants = {'_inf'; '_mum'};
% Set which visit in the longitudinal study
visits = [1,2,3,4];
% Set any special channels that we want to look at. Depreciated
% centralChans = [11,13,19,32,46,47,48,49,56];
centralChans = [31, 32, 33, 38, 48]; % 38, 47, 32, 31
centralOn = 0;
% Set the frequency resolution, upper and lower freuency bounds and target
% frequency
freqAndRes = [0.01, 6, 10, 8];
% Set analysis type
analysisType = 0;

addpath(genpath(rootPath));

%% 8Hz
% Full windows
timeWindow = [1 122880]; titleEx = 'Full window';
runAnalyses(rootPath, filePaths, centralOn, centralChans, targetConditions, participants, timeWindow, titleEx, freqAndRes, diskLetter, analysisType)
% createGrandAverages

% 2 min windows
timeWindow = [1 61440;61441 122880]; titleEx = '2 min windows';
runAnalyses(rootPath, filePaths, centralOn, centralChans, targetConditions, participants, timeWindow, titleEx, freqAndRes, diskLetter, analysisType)

% 1 min windows
timeWindow = [1 30720;30721 61440;61441 92160;92161 122880]; titleEx = '1 min windows';
runAnalyses(rootPath, filePaths, centralOn, centralChans, targetConditions, participants, timeWindow, titleEx, freqAndRes, diskLetter, analysisType)

% 30 sec windows
timeWindow = [1 15360;15361 30720;30721 46080;46081 61440;61441 76800;76801 92160;92161 107520;107521 122880]; titleEx = '30 sec windows';
runAnalyses(rootPath, filePaths, centralOn, centralChans, targetConditions, participants, timeWindow, titleEx, freqAndRes, diskLetter, analysisType)


end

 
runAnalyses.m

function runAnalyses(rootPath, filePaths, centralOn, centralChans, targetConditions, participants, timeWindow, titleEx, freqAndRes, diskLetter, analysisType)

addpath(genpath(rootPath));

if analysisType == 0 | analysisType == 1
    calcWindowedGrandAvgFFT(filePaths, centralOn, centralChans, targetConditions, participants, timeWindow, titleEx, freqAndRes, diskLetter)
end

if analysisType == 0 | analysisType == 2
    calcWindowedGrandAvgPLVAndEntropy(filePaths, centralOn, centralChans, targetConditions, participants, timeWindow, titleEx, freqAndRes, diskLetter)
end

if analysisType == 0 | analysisType == 3
    analyses = {'FFT'; 'PLV'; 'Entropy'};
    %     windowLens = {'Full window'; '2 min windows'; '1 min windows'; '30 sec windows'};
    %     windowSegs = [1,2,4,8];

    if strcmp(titleEx, 'Full window')
        winSeg = 1;
    elseif strcmp(titleEx, '2 min windows')
        winSeg = 2;
    elseif strcmp(titleEx, '1 min windows')
        winSeg = 4;
    elseif strcmp(titleEx, '30 sec windows')
        winSeg = 8;
    end

    for filePath = 1:length(filePaths)
        for analysis = 1:length(analyses)
            rootPath = strcat(filePaths{filePath}, 'Analysed Data\');
            destPath = strcat(rootPath, analyses{analysis}, '\', titleEx, '\');

            cd(destPath);
            for participant = 1:length(participants)
                for targetCondition = 1:length(targetConditions)
                    for w = 1:winSeg
                        searchString = strcat('*',targetConditions{targetCondition}, participants{participant}, '_', num2str(w), '_all elecs_', titleEx, '.mat');
                        files = dir(searchString);

                        for file = 1:length(files)
                            % This checks to see if we have done this before and
                            % will skip this file if we have, if we've not done
                            % even one time segment it will do them all again.
                            foundFiles = 0;
                            for checkTimeWindow = 1:winSeg
                                checkFilename = strcat('SNR_', files(file).name(1:end-4), '.mat');
                                snrDestPath = strcat(rootPath, 'SNR\', analyses{analysis}, '\', titleEx, '\');
                                if exist(strcat(snrDestPath, checkFilename))
                                    foundFiles = foundFiles + 1;
                                end
                            end

                            if foundFiles < size(timeWindow, 1)

                                fprintf(strcat('Processing - ', files(file).name(1:end-4), '\n'))

                                data = load(strcat(files(file).folder, '\', files(file).name));
                                if strcmp(analyses{analysis}, 'FFT')
                                    data = data.fft;
                                    scale = squeeze(data(2,1,:))';
                                    data = squeeze(data(1,:,:));
                                elseif strcmp(analyses{analysis}, 'PLV')
                                    data = data.plvData;
                                    scale = data.scale;
                                    data = data.data;
                                elseif strcmp(analyses{analysis}, 'Entropy')
                                    data = data.entropyData;
                                    scale = data.scale;
                                    data = data.data;
                                end

                                calcSNR(rootPath, analyses{analysis}, titleEx, files(file).name, data, scale, freqAndRes)

                            else
                                fprintf(strcat('Skipping - ', files(file).name(1:end-4), ' - already processed.\n'))
                            end
                        end
                    end
                end
            end
        end
    end
end

end
 
calcWindowedGrandAvgFFT.m

function calcWindowedGrandAvgFFT(filePaths, centralOn, centralChans, targetConditions, participants, timeWindow, titleEx, freqAndRes, diskLetter)

if nargin == 0
    % Used for testing the script without having to run a container
    diskLetter = 'F';
    rootPath = strcat(diskLetter, ':\FrontierPaperAnalysis\');
    filePaths = {strcat(diskLetter, ':\Frontiers Analysis\Stage_6_preproc_data\')};
    targetConditions = {'sbn2'; 'sbn4'; 'sbn6'; 'sbn8'; 'sbn10'; 'sbn12'};
    participants = {'_inf'; '_mum'};
    % centralChans = [11,13,19,32,46,47,48,49,56];
    centralChans = [31, 32, 33, 38, 48]; % 38, 47, 32, 31
    centralOn = 0;
    timeWindow = [1 30720;30721 61440;61441 92160;92161 122880];
    titleEx = '1 min windows';
    freqAndRes = [0.01, 6, 10, 8];
end

% General settings
titleExtra = 'all elecs_';
Fs = 512;

% Specific central electrodes?
if centralOn
    titleExtra = 'central elecs_';
end

for p = 1:length(filePaths)
    outputDir = strcat(filePaths{p}, 'Analysed Data\FFT\', titleEx, '\');
    if ~exist(outputDir); mkdir(outputDir); end

    % Cycle through all conditions and inf/mum to give all condition possibilities
    for t = 1:length(targetConditions)
        for prt = 1:length(participants)
            pptCounter = 0;
            rawFoundParticipants = {}; % Some functions take special characters, some don't so we have both lists.
            foundParticipants = {};
            targetEEG = [];
            
            % Collect all the files in either clean or raw
            cd(filePaths{p})
            filenames = dir('*.mat');
            
            % Set our target condition based on the loops
            target = strcat(targetConditions(t), participants(prt));
            
            % Loop through all files found
            for i = 1:length(filenames)
                % Only move forward if the filename matches the target
                if contains(filenames(i).name, target)

                    % This checks to see if we have done this before and
                    % will skip this file if we have, if we've not done
                    % even one time segment it will do them all again.
                    foundFiles = 0;
                    for checkTimeWindow = 1:size(timeWindow,1)
                        checkFilename = strcat('FFT_', filenames(i).name(1:end-4), '_', num2str(checkTimeWindow), '_', titleExtra, titleEx, '.mat');
                        if exist(strcat(outputDir, checkFilename))
                            foundFiles = foundFiles + 1;
                        end
                    end

                    if foundFiles < size(timeWindow, 1)

                        % Load the EEG, and sort the weirdness I've created
                        EEG = load(strcat(filenames(i).folder, '\', filenames(i).name));
    
                        try
                            EEG = EEG.EEG;
                        catch
                            if contains(filenames(i).name, 'inf')
                                EEG = EEG.inf_EEG;
                            elseif contains(filenames(i).name, 'mum')
                                EEG = EEG.mum_EEG;
                            else
                                disp(strcat('Weird variable names in ', filenames(i).name))
                            end
                        end
    
                        if length(EEG.data) < 122880
                            z = zeros(64, 122880 - length(EEG.data));
                            EEG.data = [EEG.data z];
                        end
    
                        for window = 1:size(timeWindow,1)
                            fname = filenames(i).name(1:end-4);
                            fname = strrep(fname, '_', '\_');
                            fname = strcat(fname, '\_', num2str(window));

                            if size(foundParticipants) == [0 0]
                                foundParticipants = {fname};
                                rawFoundParticipants = {filenames(i).name(1:end-4)};
                            else
                                foundParticipants(end+1) = {fname};
                                rawFoundParticipants(end+1) = {filenames(i).name(1:end-4)};
                            end
    
                            pptCounter = pptCounter + 1;
    
                            if ~centralOn
                                targetEEG(pptCounter,:,:) = EEG.data(1:64, timeWindow(window, 1):timeWindow(window, 2));
                            else
                                tempTargetEEG = [];
                                for roiChan = 1:length(centralChans)
                                    tempTargetEEG = [tempTargetEEG; EEG.data(roiChan,timeWindow(window, 1):timeWindow(window, 2))];
                                end
                                targetEEG(pptCounter,:,:) = tempTargetEEG;
                            end
                        end
                    else
                        fprintf(strcat('Skipping - ', filenames(i).name(1:end-4), ' - already processed.\n'))
                    end
                end
            end

            fileCount = 1;
            for fftPPT = 1:size(targetEEG, 1)
                fft = [];
                for fftElec = 1:size(targetEEG, 2)
                    [fftRes, fftHzScale, dBfft] = myFFT(squeeze(targetEEG(fftPPT, fftElec, :))',Fs,0,50);
                    [minValue, closestIndex] = min(abs(fftHzScale-25)); % We're filtering at 25Hz so no point keeping data past this point.
                    fftRes = fftRes(1:closestIndex);
                    fftHzScale = fftHzScale(1:closestIndex);
                    dBfft = dBfft(1:closestIndex);

                    fft(1, fftElec, :) = fftRes;
                    fft(2, fftElec, :) = fftHzScale;
                    fft(3, fftElec, :) = dBfft;

                end

                filename = strcat('FFT_', rawFoundParticipants{fftPPT}, '_', num2str(fileCount), '_', titleExtra, titleEx, '.mat');
                savename = strcat(outputDir, filename);
                save(savename, 'fft');

                if fileCount == size(timeWindow, 1)
                    fileCount = 1;
                else
                    fileCount = fileCount + 1;
                end
            end
        end
    end
end

end
 
myFFT.m
function [fftFFR, HzScale, dBfft] = myFFT(data,Fs,toPlot,xlimVar,fileName)
%
% Compute and plot FFT
%
% Version 1 - September 2017
%
% Tim Schoof - t.schoof@ucl.ac.uk
% ----------------------------------
if nargin < 5
   fileName = 'file name not specified';
elseif nargin < 4
    xlimVar = 1200;
end
% zero-pad data
 data = [data, zeros(1,length(data))];
% compute FFT
L=length(data);   
NFFT = 2^nextpow2(L);       % Next power of 2 from length of y
Y = fft(data,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1);
fftFFR = 2*abs(Y(1:NFFT/2+1));
% single-sided amplitude spectrum re peak amplitude - not used for plotting
dBref = 1;
dBfft = 20*log10(fftFFR/dBref);
HzScale = [0:(Fs/2)/(length(fftFFR)-1):round(Fs/2)]'; % frequency 'axis'
% plot FFT (if required)
if toPlot == 1
    plot(HzScale,fftFFR,'LineWidth',2)
    xlim([0 xlimVar])
    set(0, 'DefaulttextInterpreter', 'none')
    title(['',fileName,''])
    ylabel('Spectral magnitude')
    xlabel('Frequency (Hz)')
end
 
calcWindowedGrandAvgPLVAndEntropy.m
function calcWindowedGrandAvgPLVAndEntropy(filePaths, centralOn, centralChans, targetConditions, participants, timeWindow, titleEx, freqAndRes, diskLetter)

if nargin == 0
    diskLetter = 'F';
    rootPath = strcat(diskLetter, ':\FrontierPaperAnalysis\');
    filePaths = {strcat(diskLetter, ':\Frontiers Analysis\Stage_6_preproc_data\')};
    targetConditions = {'sbn2'; 'sbn4'; 'sbn6'; 'sbn8'; 'sbn10'; 'sbn12'; 'jbn5p8'; 'jbn10p8'; 'jbn20p8'; 'jmb5p8'; 'jmb10p8'; 'jmb20p8'; 'joob5p8'; 'joob10p8'; 'joob20p8'}; % {'ft1'};
    participants = {'_inf'; '_mum'};
    % centralChans = [11,13,19,32,46,47,48,49,56];
    centralChans = [31, 32, 33, 38, 48]; % 38, 47, 32, 31
    centralOn = 0;
    timeWindow = [1 30720;30721 61440;61441 92160;92161 122880];
    titleEx = '1 min windows';
    freqAndRes = [0.01, 6, 10, 8];
end

% General settings
titleExtra = 'all elecs_';
Fs = 512;
Colours=GenColours;

if centralOn
    titleExtra = 'central elecs_';
end

for p = 1:length(filePaths)
    plvOutputDir = strcat(filePaths{p}, 'Analysed Data\PLV\', titleEx, '\');
    entropyOutputDir = strcat(filePaths{p}, 'Analysed Data\Entropy\', titleEx, '\');

    if ~exist(plvOutputDir); mkdir(plvOutputDir); end
    if ~exist(entropyOutputDir); mkdir(entropyOutputDir); end

    for t = 1:length(targetConditions)
        for prt = 1:length(participants)
            pptCounter = 0;
            rawFoundParticipants = {}; % Some functions take special characters, some don't so we have both lists.
            foundParticipants = {};
            targetEEG = [];
            
            cd(filePaths{p})
            filenames = dir('*.mat');
            target = strcat(targetConditions(t), participants(prt));
            
            for i = 1:length(filenames)
                if contains(filenames(i).name, target)

                    % This checks to see if we have done this before and
                    % will skip this file if we have, if we've not done
                    % even one time segment it will do them all again.
                    foundFiles = 0;
                    for checkTimeWindow = 1:size(timeWindow,1)
                        checkFilename = strcat('PLV_', filenames(i).name(1:end-4), '_', num2str(checkTimeWindow), '_', titleExtra, titleEx, '.mat');
                        if exist(strcat(plvOutputDir, checkFilename))
                            foundFiles = foundFiles + 1;
                        end
                    end

                    if foundFiles < size(timeWindow, 1)
                        EEG = load(strcat(filenames(i).folder, '\', filenames(i).name));
                        
                        try
                            EEG = EEG.EEG3;
                        catch
                            try
                                EEG = EEG.EEG;
                            catch
                                if contains(filenames(i).name, 'inf')
                                    EEG = EEG.inf_EEG;
                                elseif contains(filenames(i).name, 'mum')
                                    EEG = EEG.mum_EEG;
                                else
                                    disp(strcat('Weird variable names in ', filenames(i).name))
                                end
                            end
                        end
    
                        if length(EEG.data) < 122880
                            z = zeros(64, 122880 - length(EEG.data));
                            EEG.data = [EEG.data z];
                        end
    
                        for window = 1:size(timeWindow,1)
                            fname = filenames(i).name(1:end-4);
                            fname = strrep(fname, '_', '\_');
                            fname = strcat(fname, '\_', num2str(window));
        
                            if size(foundParticipants) == [0 0]
                                foundParticipants = {fname};
                                rawFoundParticipants = {filenames(i).name(1:end-4)};
                            else
                                foundParticipants(end+1) = {fname};
                                rawFoundParticipants(end+1) = {filenames(i).name(1:end-4)};
                            end
                            
                            pptCounter = pptCounter + 1;
                            
                            if ~centralOn
                                targetEEG(pptCounter,:,:) = EEG.data(1:64, timeWindow(window, 1):timeWindow(window, 2));
                            else
                                tempTargetEEG = [];
                                for roiChan = 1:length(centralChans)
                                    tempTargetEEG = [tempTargetEEG; EEG.data(roiChan,timeWindow(window, 1):timeWindow(window, 2))];
                                end
                                targetEEG(pptCounter,:,:) = tempTargetEEG;
                            end
                        end
                    else
                        fprintf(strcat('Skipping - ', filenames(i).name(1:end-4), ' - already processed.\n'))
                    end
                end
            end

            try size(targetEEG) ~= [0 0];
                disp(strcat('No data for\t', target{:}))
            catch

                fileCount = 1;
                for tempPPT = 1:size(targetEEG, 1)
                    entropyAcrossChan = [];
                    plvAcrossChan = [];
                    
                    for tempChan = 1:size(targetEEG,2)
                        EEG_sig = targetEEG(tempPPT,tempChan,:);
                        EEG_sig = squeeze(EEG_sig)';
                        listentropy = [];
                        listplv = [];

                        %frequencies to test lowFreq:step:highFreq
                        freqToTest = freqAndRes(2):freqAndRes(1):freqAndRes(3);

                        for k = freqToTest
                            F = k;

                            ts = 0:1/Fs:(length(EEG_sig)/Fs)-1/Fs; % time base
                            sig_STIM = sin(2.0*pi*F*ts);
                            phase_sig1 = angle(hilbert(sig_STIM));
                            phase_sig2 = angle(hilbert(EEG_sig));

                            [~, Ntrials] = size(phase_sig1);
                            % compute PLV
                            e = exp(1i*(phase_sig1 - phase_sig2));
                            plv = abs(sum(e,2)) / Ntrials;

                            % compute Entropy
                            nb_bins = 80;
                            entropy = Bentropy(phase_sig1 - phase_sig2,nb_bins);
                            entropy = (log2(Ntrials)-entropy)/log2(Ntrials);

                            listentropy = [listentropy; entropy];
                            listplv = [listplv; plv];
                        end

                        entropyAcrossChan = [entropyAcrossChan; listentropy'];
                        plvAcrossChan = [plvAcrossChan; listplv'];

                    end
                    entropyData.data = entropyAcrossChan;
                    entropyData.scale = freqToTest;
                    
                    filename = strcat('Entropy_', rawFoundParticipants{tempPPT}, '_', num2str(fileCount), '_', titleExtra, titleEx, '.mat');
                    savename = strcat(entropyOutputDir, filename);
                    save(savename, 'entropyData');

                    plvData.data = plvAcrossChan;
                    plvData.scale = freqToTest;
                    
                    filename = strcat('PLV_', rawFoundParticipants{tempPPT}, '_', num2str(fileCount), '_', titleExtra, titleEx, '.mat');
                    savename = strcat(plvOutputDir, filename);
                    save(savename, 'plvData');

                    if fileCount == size(timeWindow, 1)
                        fileCount = 1;
                    else
                        fileCount = fileCount + 1;
                    end
                end
            end
        end
    end  
end

end

 
Bentropy.m
function [H] = Bentropy(x,nb_bins)
[counts,binCenters] = hist(x,nb_bins);
binWidth = diff(binCenters);
binWidth = [binWidth(end),binWidth];
nz = counts>0;
frequency = counts(nz)/sum(counts(nz));
H = -sum(frequency.*log(frequency./binWidth(nz)));
end
 
calcSNR.m
function calcSNR(rootPath, analysis, win, filename, data, scale, freqAndRes)

% Data is expected to be chans x time/spectra, you can't do multiple ppts
% at the same time, so either one by one or averaged. If you average first
% you are likely to ruin the data, fyi.

if nargin == 0
    type = 0;
    freqAndRes = [0.01 6 10 8];
    fft = load('C:\Users\James Ives\Documents\Frontiers Analysis\Stage_6_preproc_data\Analysed Data\FFT\Full window\FFT_1021_1_sbn8_inf_1_all elecs_Full window.mat');
    scale = squeeze(fft.fft(2,1,:))';
    data = squeeze(fft.fft(1,:,:));
end

plotMe = 0; 

SNROutputDir = strcat(rootPath, 'SNR\', analysis, '\', win, '\');
if ~exist(SNROutputDir); mkdir(SNROutputDir); end

Colours=GenColours;

% This checks to see if we have done this before and
% will skip this file if we have, if we've not done
% even one time segment it will do them all again.
if strcmp(win, 'Full window')
    timeWindow = 1;
elseif strcmp(win, '2 min windows')
    timeWindow = 2;
elseif strcmp(win, '1 min windows')
    timeWindow = 4;
elseif strcmp(win, '30 sec windows')
    timeWindow = 8;
end

foundFiles = 0;
for checkTimeWindow = 1:timeWindow
    checkFilename = strcat('SNR_', filename(1:end-4), '_', num2str(checkTimeWindow), '_all elecs_', win, '.mat');
    if exist(strcat(SNROutputDir, checkFilename))
        foundFiles = foundFiles + 1;
    end
end

% calculate SNR measure
% from sofie's thesis:
% the amplitude spectrum was computed with a high spectral resolution
% (0.017 Hz, 1/ 59.38 s) resulting in a very high signal-to-noise ratio (SNR)
% (Norcia et al., 2015; Regan, 1989).

% SNR spectra were computed for each electrode by dividing the value at
% each frequency bin by the average value of the 20 neighboring frequency
% bins (12 bins on each side, i.e., 24 bins, but excluding the 2 bins
% directly adjacent and the 2 bins with the most extreme values).
if plotMe
    fig(1) = figure('units','normalized','outerposition',[0 0 1 1]);
end
clear SNR_Ret
for chans=1:64
    freqcnt=0;
    forlegend{chans}=num2str(chans);
    
    for freq=freqAndRes(2):freqAndRes(1):freqAndRes(3)
        % Finds the closes index to the target freq (ignore the minValue, unused)
        [minValue,closestIndex] = min(abs(freq-scale));

        % Calcs the bounds of the indices to be used for SNR, trims if necessary
        beforestart=closestIndex-12;
        if beforestart<1
            beforestart=1;
        end

        afterend=closestIndex+12;
        if afterend>length(data)
            afterend=length(data);
        end

        % Finds the index values before and after the target index,
        % excluding the indices adjacent
        beforeend=closestIndex-2;
        afterstart=closestIndex+2;


        ValueAtFreq=data(chans, closestIndex);
        ComparisonVals=horzcat(data(chans, beforestart:beforeend),data(chans, afterstart:afterend))';
        % drop off the value showing the biggest absolute difference from the ValueAtFreq
        [maxValue,closestIndex2] = max(abs(ValueAtFreq-ComparisonVals));
        ComparisonVals(closestIndex2)=NaN;
        % do it again
        [maxValue,closestIndex3] = max(abs(ValueAtFreq-ComparisonVals));
        ComparisonVals(closestIndex3)=NaN;
        freqcnt=freqcnt+1;
        SNR_Ret(freqcnt,1,chans)=freq;
        SNR_Ret(freqcnt,2,chans)=ValueAtFreq/mean(ComparisonVals, 'omitnan');
    end
    squeeze(SNR_Ret(:,1:2,chans));
    if plotMe
        plot(squeeze(SNR_Ret(:,1,chans)'),squeeze(SNR_Ret(:,2,chans)'),'Color',Colours(chans,:))
        hold on
    end
end


if plotMe
    mSNR = mean(SNR_Ret,3);
    plot(squeeze(SNR_Ret(:,1,chans)),mSNR(:,2), 'k', 'LineWidth', 3);
    pfTitle = strrep(filename, '_', '\_'); pfTitle = strcat('SNR_', pfTitle(1:end-4));
    ylabel('SNR'); xlabel(strcat('freq (Hz)')); title(pfTitle);
end

SNR=SNR_Ret;

savename = strcat('SNR_', filename); save(strcat(SNROutputDir, '\', savename), 'SNR');

end
 
randPermutationsByPptAnalysis.m
function randPermuationsByPptAnalysis(freqAndRes)

if nargin == 0
    diskLetter = 'F';
    rootPath = strcat(diskLetter, ':\Frontiers Analysis\Stage_6_preproc_data\Analysed Data\');
    targetConditions = {'sbn8'; }; % {'ft1'}; 'jbn5p8'; 'jbn10p8'; 'jbn20p8'; 'jmb5p8'; 'jmb10p8'; 'jmb20p8'; 'joob5p8'; 'joob10p8'; 'joob20p8'
    participants = {'_inf'; '_mum'};
    visits = [1,2,3,4];
    analyses = {'PLV'; 'Entropy'; 'FFT'};
    titleExs = {'Full window'}; % ; '2 min windows'; '1 min windows'; '30 sec windows'
    winLens = {'*by ppt_Full window.mat'}; % ; '*by ppt_2 min windows.mat'; '*by ppt_1 min windows.mat'; '*by ppt_30 sec windows.mat'
    useAvgs = 0;
    numPerms = 10000;
    freqAndRes = [0.01, 6, 10, 8];
    scale = freqAndRes(2):freqAndRes(1):freqAndRes(3);
end

rng default
plotMe = 0;

randDir = strcat(rootPath, 'Rand Perm data\');
if ~exist(randDir); mkdir(randDir); end
if ~exist(strcat(randDir, 'PPT\')); mkdir(strcat(randDir, 'PPT\')); end
if ~exist(strcat(randDir, 'PNGs\')); mkdir(strcat(randDir, 'PNGs\')); end
if ~exist(strcat(randDir, 'Significance\')); mkdir(strcat(randDir, 'Significance\')); end

for tEx = 1:length(titleExs)
    titleEx = titleExs{tEx};
    fprintf(strcat('Processing - time window - ', titleEx, '\n'))

    if strcmp(titleEx, 'Full window')
        winSeg = 1;
    elseif strcmp(titleEx, '2 min windows')
        winSeg = 2;
    elseif strcmp(titleEx, '1 min windows')
        winSeg = 4;
    elseif strcmp(titleEx, '30 sec windows')
        winSeg = 8;
    end

    for analysis = 1:length(analyses)
        fprintf(strcat('Processing - analyis - ', analyses{analysis}, '\n'))

        if plotMe
            pptTitle = strcat('Random permutation significance analysis ', analyses{analysis}, '_', titleEx, ' data ');
            pptTitle2 = strcat('Avg random permutation significance analysis ', analyses{analysis}, '_', titleEx, ' data ');
            t = now; d = join(split(string(datetime(t,'ConvertFrom','datenum')),':'));
            [ppt] = createPresentation(strcat(randDir, 'PPT\'), strcat(pptTitle, ' - ', d));
            [ppt2] = createPresentation(strcat(randDir, 'PPT\'), strcat(pptTitle2, ' - ', d));
        end

        sigList = {};

        destPath = strcat(rootPath, 'SNR\', analyses{analysis}, '\', titleEx, '\');

        for t = 1:length(targetConditions)
            fprintf(strcat('Processing - condition - ', targetConditions{t}, '\n'))
            for prt = 1:length(participants)
                fprintf(strcat('Processing - participant -', participants{prt}, '\n'))
                for v = 1:length(visits)
                    fprintf(strcat('Processing - visit - ', num2str(visits(v)), '\n'))
                    for w = 1:winSeg
                        fprintf(strcat('Processing - window segment - ', num2str(w), '\n'))

                        searchString = strcat('*', num2str(visits(v)), '_', targetConditions{t}, participants{prt}, '_', num2str(w), '_all elecs_', titleEx, '.mat');
                        pfTitle = strcat('_v', num2str(visits(v)), '_', targetConditions{t}, participants{prt}, '_w', num2str(w), '_all elecs_', titleEx);
                        rawTitle = pfTitle;
                        pfTitle = strrep(pfTitle, '_', '\_');

                        cd(destPath);
                        files = dir(searchString);

                        if length(files) > 1
                            filenames = [];
                            if size(files, 1) > 0
                                % Collect all the data per target
                                % condition/participant type/visit/window
                                % length/type of analysis
                                for file = 1:length(files)
                                    data = load(strcat(files(file).folder, '\', files(file).name));
                                    %                             if strcmp(analyses{analysis}, 'FFT')
                                    %                                 data = data.fft;
                                    %                                 scale = squeeze(data(2,1,:))';
                                    %                                 data = squeeze(data(1,:,:));
                                    %                             elseif strcmp(analyses{analysis}, 'PLV')
                                    %                                 data = data.plvData;
                                    %                                 scale = data.scale;
                                    %                                 data = data.data;
                                    %                             elseif strcmp(analyses{analysis}, 'Entropy')
                                    %                                 data = data.entropyData;
                                    %                                 scale = data.scale;
                                    %                                 data = data.data;
                                    %                             end
                                    data = data.SNR;
                                    scale = squeeze(data(:,1,1))';
                                    data = squeeze(data(:,2,:))';

                                    gData(file,:,:) = data;

                                    filenames = [filenames; {files(file).name(1:end-4)}];

                                end
                            end

                            % Calc num of participants and average over
                            % electrodes
                            nParticipants = size(gData,1);
                            mData = squeeze(mean(gData, 2));

                            % Using the data averaged over electrodes to find
                            % the test datum (at the target frequency) and
                            % create the 10,000 permutations of shuffled data.
                            for p = 1:size(mData,1)

                                testData = squeeze(mData(p, :));
                                test8Hz = testData(201);

                                compSet = [];

                                for perm = 1:numPerms
                                    tempPerm = testData(randperm(length(testData)));

                                    compSet = [compSet tempPerm(201)];
                                end

                                above8Hz = compSet(test8Hz<compSet);
%                                 pc = prctile(compSet, 95)
                                invpc = (100-invprctile(compSet, test8Hz))/100;
%                                 nLess = sum(compSet < test8Hz);
%                                 nEqual = sum(compSet == test8Hz);
%                                 centile = 100 * (nLess + 0.5*nEqual) / length(compSet);

                                %                             returnPermData = {size(aboveTest8Hz,2), test8Hz, compSet};
                                %                             allRandPermData(ppt, elec, :) =  [returnPermData];

                                allCompSets(p, :) = [compSet];
%                                 allTestData(p, :) = [test8Hz, size(aboveTest8Hz, 2), (size(aboveTest8Hz, 2)/numPerms)];
                                allTestData(p, :) = {test8Hz, invpc, above8Hz};

                                if plotMe
                                    fig(1) = figure('units','normalized','outerposition',[0 0 1 1], 'visible', 'off'); %

                                    plot(scale, mData(p,:));
                                    ylabel(strcat('SNR of ', analyses{analysis})); xlabel(strcat('freq (Hz)')); xline(8, '--', 'Stimulation Freq');
                                    s = strrep(files(p).name(1:end-4), '_','\_'); title(s);

                                    [ppt] = addImgToPresentation(strcat(randDir, 'PNGs\'), ppt, strcat(files(p).name(1:end-4)), fig(1));
                                    delete(fig);
                                end

                            end
                            if plotMe
                                fig(2) = figure('units','normalized','outerposition',[0 0 1 1], 'visible', 'off'); % , 'visible', 'off'

                                plot(scale, squeeze(mean(mData,1)), 'k', 'LineWidth', 2);
                                ylabel(strcat('Avg SNR of ', analyses{analysis})); xlabel(strcat('freq (Hz)')); xline(8, '--', 'Stimulation Freq');
                                title(strcat('Avg SNR of ', analyses{analysis}, '\_v', num2str(visits(v)), '\_\', targetConditions{t}, '\', participants{prt}, '\_w', num2str(w), '\_all elecs\_', titleEx));

                                [ppt2] = addImgToPresentation(strcat(randDir, 'PNGs\'), ppt2, strcat('Avg SNR of ', analyses{analysis}, '_v', num2str(visits(v)), '_', targetConditions{t}, participants{prt}, '_w', num2str(w), '_all elecs_', titleEx), fig(2));
                                delete(fig);
                            end

                            save(strcat(randDir, analyses{analysis}, '_', targetConditions{t}, '_', num2str(visits(v)), participants{prt}, '_', num2str(w),  '_all elecs by ppt_', titleEx, '_compSet.mat'), "allCompSets")
                            save(strcat(randDir, analyses{analysis}, '_', targetConditions{t}, '_', num2str(visits(v)), participants{prt}, '_', num2str(w),  '_all elecs by ppt_', titleEx, '_testSet.mat'), "allTestData")

                            percentiles = [];
                            for x = 1:size(allTestData,1)
                                percentiles = [percentiles; allTestData{x,2}];
                            end
                            sig = mean(percentiles);
%                             sig = sum(sum(allTestData(:,2)))/(numPerms * size(mData,1));
                            sigList = [sigList; {strcat(analyses{analysis}, '_v', num2str(visits(v)), '_', targetConditions{t}, participants{prt}, '_w', num2str(w), '_all elecs_', titleEx), sig, filenames, allTestData}];

                            clear allCompSets
                            clear allTestData
                            clear gData
                        else
                            disp('Not enough files, skipping')
                        end
                    end
                end
            end
        end

        save(strcat(randDir, 'Significance\', analyses{analysis}, '_', titleEx, '.mat'), "sigList")

        if plotMe
            close(ppt);
            close(ppt2);
        end
    end
end

end
 
whichFreqOneWayANOVAN.m
function whichFreqOneWayANOVAN()


if nargin == 0
    diskLetter = 'F';
    rootPath = strcat(diskLetter, ':\Frontiers Analysis\Stage_6_preproc_data\Analysed Data\');
    targetConditions = {'sbn2'; 'sbn4'; 'sbn6'; 'sbn8'; 'sbn10'; 'sbn12'}; % {'ft1'}; 'jbn5p8'; 'jbn10p8'; 'jbn20p8'; 'jmb5p8'; 'jmb10p8'; 'jmb20p8'; 'joob5p8'; 'joob10p8'; 'joob20p8'
    participants = {'_inf'; '_mum'};
    visits = [1,2,3,4];
    analyses = {'FFT'; 'PLV'; 'Entropy'};
    titleExs = {'Full window'}; % ; '2 min windows'; '1 min windows'; '30 sec windows'
    winLens = {'*by ppt_Full window.mat'}; % ; '*by ppt_2 min windows.mat'; '*by ppt_1 min windows.mat'; '*by ppt_30 sec windows.mat'
    useAvgs = 0;
    numPerms = 10000;
    freqAndRes = [0.01, 6, 10, 8];
    scale = freqAndRes(2):freqAndRes(1):freqAndRes(3);
end

rng default
% Doesn't work at the moment, we can't get the fig handle our of anova
plotMe = 1;

chanNames = {'1-Fp1' '1-AF7' '1-AF3' '1-F1' '1-F3' '1-F5' '1-F7' '1-FT7' '1-FC5' '1-FC3' '1-FC1' '1-C1' '1-C3' '1-C5' '1-T7' '1-TP7' '1-CP5' '1-CP3' '1-CP1' '1-P1' '1-P3' '1-P5' '1-P7' '1-P9' '1-PO7' '1-PO3' '1-O1' '1-Iz' '1-Oz' '1-POz' '1-Pz' '1-CPz' '1-Fpz' '1-Fp2' '1-AF8' '1-AF4' '1-AFz' '1-Fz' '1-F2' '1-F4' '1-F6' '1-F8' '1-FT8' '1-FC6' '1-FC4' '1-FC2' '1-FCz' '1-Cz' '1-C2' '1-C4' '1-C6' '1-T8' '1-TP8' '1-CP6' '1-CP4' '1-CP2' '1-P2' '1-P4' '1-P6' '1-P8' '1-P10' '1-PO8' '1-PO4' '1-O2'}';

randDir = strcat(rootPath, 'Rand Perm data\');
if ~exist(randDir); mkdir(randDir); end
if ~exist(strcat(randDir, 'PPT\')); mkdir(strcat(randDir, 'PPT\')); end
if ~exist(strcat(randDir, 'PNGs\')); mkdir(strcat(randDir, 'PNGs\')); end
if ~exist(strcat(randDir, 'Significance\')); mkdir(strcat(randDir, 'Significance\')); end



for tEx = 1:length(titleExs)
    titleEx = titleExs{tEx};
    fprintf(strcat('Processing - time window - ', titleEx, '\n'))

    if strcmp(titleEx, 'Full window')
        winSeg = 1;
    elseif strcmp(titleEx, '2 min windows')
        winSeg = 2;
    elseif strcmp(titleEx, '1 min windows')
        winSeg = 4;
    elseif strcmp(titleEx, '30 sec windows')
        winSeg = 8;
    end

    if plotMe
        pptTitle = strcat('ANOVAN Random permutation significance analysis ', titleEx, ' data ');
        t = now; d = join(split(string(datetime(t,'ConvertFrom','datenum')),':'));
        [ppt] = createPresentation(strcat(randDir, 'PPT\'), strcat(pptTitle, ' - ', d));
    end

    for analysis = 1:length(analyses)
        fprintf(strcat('Processing - analyis - ', analyses{analysis}, '\n'))

        sigList = {};
        testData = {};

        destPath = strcat(rootPath, 'SNR\', analyses{analysis}, '\', titleEx, '\');

        for t = 1:length(targetConditions)
            fprintf(strcat('Processing - condition - ', targetConditions{t}, '\n'))

            if strcmp(targetConditions{t}, 'sbn2')
                freqAndRes = [0.01, 1, 4, 2];
            elseif strcmp(targetConditions{t}, 'sbn4')
                freqAndRes = [0.01, 2, 6, 4];
            elseif strcmp(targetConditions{t}, 'sbn6')
                freqAndRes = [0.01, 4, 8, 6];
            elseif strcmp(targetConditions{t}, 'sbn8')
                freqAndRes = [0.01, 6, 10, 8];
            elseif strcmp(targetConditions{t}, 'sbn10')
                freqAndRes = [0.01, 8, 12, 10];
            elseif strcmp(targetConditions{t}, 'sbn12')
                freqAndRes = [0.01, 10, 14, 12];
            end
            scale = freqAndRes(2):freqAndRes(1):freqAndRes(3);

            for prt = 1:length(participants)
                fprintf(strcat('Processing - participant -', participants{prt}, '\n'))
                for v = 1:length(visits)
                    fprintf(strcat('Processing - visit - ', num2str(visits(v)), '\n'))
                    for w = 1:winSeg
                        fprintf(strcat('Processing - window segment - ', num2str(w), '\n'))

                        searchString = strcat('*', num2str(visits(v)), '_', targetConditions{t}, participants{prt}, '_', num2str(w), '_all elecs_', titleEx, '.mat');
                        pfTitle = strcat('_v', num2str(visits(v)), '_', targetConditions{t}, participants{prt}, '_w', num2str(w), '_all elecs_', titleEx);
                        rawTitle = pfTitle;
                        pfTitle = strrep(pfTitle, '_', '\_');

                        cd(destPath);
                        files = dir(searchString);

                        if length(files) > 1
                            filenames = [];
                            if size(files, 1) > 0
                                % Collect all the data per target
                                % condition/participant type/visit/window
                                % length/type of analysis
                                for file = 1:length(files)
                                    data = load(strcat(files(file).folder, '\', files(file).name));
                                    data = data.SNR;
                                    scale = squeeze(data(:,1,1))';
                                    data = squeeze(data(:,2,:))';

                                    gData(file,:,:) = data;

                                    filenames = [filenames; {files(file).name(1:end-4)}];

                                end
                            end

                            % Calc num of participants and average over
                            % electrodes
                            nParticipants = size(gData,1);
%                             gData = squeeze(mean(gData, 1));
                            
                            if strcmp(targetConditions{t}, 'sbn2')
                                temp = gData(:, :, 101); temp = temp(:); %
                                testData(t, prt, v, w, :) = {temp};
                            else
                                temp = gData(:, :, 01); temp = temp(:); %
                                testData(t, prt, v, w, :) = {temp};
                            end
                            clear gData
                        end
                    end
                end
            end
        end

%         input(strcat('Hit enter to run ANOVA on inf data for, ', analyses{analysis}, ':'), 's')
%         infData = padcat([testData{1,1}], [testData{2,1}], [testData{3,1}], [testData{4,1}], [testData{5,1}], [testData{6,1}]);
        allInfData = [];
        allInfConditionNames = [];
        allInfElectrodeNames = [];

        allMumData = [];
        allMumConditionNames = [];
        allMumElectrodeNames = [];

        for row = 1:6

            infDataCol1 = [testData{row,1}];
            mumDataCol1 = [testData{row,2}];
            infDataCol2 = [];
            mumDataCol2 = [];

            for condGroup = 1:size(infDataCol1,1)/64
                allInfElectrodeNames = [allInfElectrodeNames; chanNames];
            end
            for chanGroup = 1:size(infDataCol1)
                allInfConditionNames = [allInfConditionNames; targetConditions(row)];
            end

            for condGroup = 1:size(mumDataCol1,1)/64
                allMumElectrodeNames = [allMumElectrodeNames; chanNames];
            end
            for chanGroup = 1:size(mumDataCol1)
                allMumConditionNames = [allMumConditionNames; targetConditions(row)];
            end


            allInfData = [allInfData; infDataCol1];
            allMumData = [allMumData; mumDataCol1];
%             allGroups = [allGroups; infDataCol2];
        end

        [infP, inftbl, infStats] = anovan(allInfData, {allInfConditionNames, allInfElectrodeNames}, 'display', 'on', 'varNames', {'Conditions', 'Electrodes'});
        if plotMe
            figHandles = findall(groot, 'Type', 'figure');
            fig = figHandles(1);
%             fig.Position = [0 0 1550 775];
            [ppt] = addImgToPresentation(strcat(randDir, 'PNGs\'), ppt, strcat('ANOVA stats ', analyses{analysis}, '_', titleEx, '_inf'), fig);
            delete(fig);
%             fig = figHandles(2);
%             fig.Position = [0 0 1550 775];
%             [ppt] = addImgToPresentation(strcat(randDir, 'PNGs\'), ppt, strcat('ANOVA graph ', analyses{analysis}, '_', titleEx, '_inf'), fig);
%             delete(fig);
        end
%         input(strcat('Hit enter to run multiple comparisons on inf data for, ', analyses{analysis}, ':'), 's')
        [infComparison,infMean,infFigHandle,infGroupNames] = multcompare(infStats);
        if plotMe
            [ppt] = addImgToPresentation(strcat(randDir, 'PNGs\'), ppt, strcat('ANOVAN multiple comparisons ', analyses{analysis}, '_', titleEx, '_inf'), infFigHandle);
            delete(infFigHandle);
        end

%         input(strcat('Hit enter to run ANOVA on mum data for, ', analyses{analysis}, ':'), 's')
        mumData = padcat([testData{1,2}], [testData{2,2}], [testData{3,2}], [testData{4,2}], [testData{5,2}], [testData{6,2}]);
        [mumP, mumtbl, mumStats] = anovan(allMumData, {allMumConditionNames, allMumElectrodeNames}, 'display', 'on', 'varNames', {'Conditions', 'Electrodes'});
        if plotMe
            figHandles = findall(groot, 'Type', 'figure');
            fig = figHandles(1);
%             fig.Position = [0 0 1550 775];
            [ppt] = addImgToPresentation(strcat(randDir, 'PNGs\'), ppt, strcat('ANOVAN ', analyses{analysis}, '_', titleEx, '_mum'), fig);
            delete(fig);
%             fig = figHandles(2);
%             fig.Position = [0 0 1550 775];
%             [ppt] = addImgToPresentation(strcat(randDir, 'PNGs\'), ppt, strcat('ANOVA graph ', analyses{analysis}, '_', titleEx, '_inf'), fig);
%             delete(fig);
        end
%         input(strcat('Hit enter to run multiple comparisons on mum data for, ', analyses{analysis}, ':'), 's')
        [mumComparison,mumMean,mumFigHandle,mumGroupNames] = multcompare(mumStats);
        if plotMe
            [ppt] = addImgToPresentation(strcat(randDir, 'PNGs\'), ppt, strcat('ANOVA multiple comparisons ', analyses{analysis}, '_', titleEx, '_mum'), mumFigHandle);
            delete(mumFigHandle);
        end


    end
if plotMe
    close(ppt);
end
end

end