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