OSCR

40 Hz audiovisual stimulation improves sustained attention and related brain oscillations.

Code ↔ Paper

13 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 13 matches · 4 of them tie a paragraph to a whole file, not to given lines: weak matches, whose lines are not tinted
  1. [1] § Results › 40 Hz flicker is correlated with decreases in delta activity during a vigilance task ↔ EEG Analysis/S4_ALLSubjects_NoCut_NoNotch_2_100hz.m, lines 4300–4333 · score 0.87 · 13–30 Hz, 30–37 Hz, 8–13 Hz, 39–41 Hz, power spectral densities, 1–4 Hz
  2. [2] § Results › 40 Hz flicker increases low-alpha functional connectivity which is correlated with better behavior performance ↔ EEG Analysis/S4_ALLSubjects_NoCut_NoNotch_2_100hz.m, lines 2058–2100 · score 0.78 · 10–13 Hz, top quartile, upper alpha, lower alpha, 8–10 Hz, channel pairs
  3. [3] § Results › 40 Hz flicker increases low-alpha functional connectivity which is correlated with better behavior performance ↔ EEG Analysis/S4_ALLSubjects_NoCut_NoNotch_2_100hz.m, lines 2058–2100 · score 0.75 · 10–13 Hz, top quartile, upper alpha, lower alpha, 8–10 Hz, channel pair
  4. [4] § Results › 40 Hz flicker is correlated with decreases in delta activity during a vigilance task ↔ EEG Analysis/S4_ALLSubjects_NoCut_NoNotch_2_100hz.m, lines 3178–3220 · score 0.74 · 30–37 Hz, 8–13 Hz, 39–41 Hz, power spectral density, 4–8 Hz, PSD
  5. [5] § Methods › EEG data analyses ↔ EEG Analysis/S3_CompletePreprocessing.m, lines 14–46 · score 0.70 · high pass filter, EEGLAB, 4 seconds, preprocessed, ICA, epochs
  6. [6] § Methods › EEG data analyses ↔ EEG Analysis/Functions/GenMatCode-main/Numeric/SignalProcessing/power/welchSpecLuTrial.m, the whole file · a weak match · score 0.67 · Signal Processing, power spectral density, Welch, overlapping, windows
  7. [7] § Methods › EEG data analyses ↔ EEG Analysis/Functions/GenMatCode-main/Numeric/SignalProcessing/GT_welchPsdTrial.m, the whole file · a weak match · score 0.67 · Signal Processing, power spectral density, Welch, overlapping, windows
  8. [8] § Methods › Statistical approach ↔ EEG Analysis/Functions/GenMatCode-main/Statistics/ANOVAandMixEffect/myfriedman.m, lines 1–83 · score 0.65 · Kruskal Wallis, Post hoc, ANOVA, mixed, alpha
  9. [9] § Methods › EEG data analyses ↔ EEG Analysis/S0_ConvertBDFtoSet_and_Preprocess_for_Syncing.m, lines 14–44 · score 0.61 · high pass filter, 4 seconds, preprocessed, ICA, epochs, noise
  10. [10] § Methods › EEG data analyses ↔ EEG Analysis/Functions/GenMatCode-main/Numeric/SignalProcessing/wpli_TrialIndex.m, the whole file · a weak match · score 0.60 · weighted phase lag, volume conduction, WPLI, connectivity, signal
  11. [11] § Results › 40 Hz flicker increases low-alpha functional connectivity which is correlated with better behavior performance ↔ EEG Analysis/S4_ALLSubjects_NoCut_NoNotch_2_100hz.m, lines 1745–1808 · score 0.56 · High functional connectivity, alpha band, channel pairs, permutation, frequency band, WPLI
  12. [12] § Results › 40 Hz flicker improved accuracy and reaction time in a vigilance task ↔ EEG Analysis/Functions/GenMatCode-main/Statistics/ANOVAandMixEffect/myfriedman.m, lines 1–83 · score 0.56 · Kruskal Wallis, post hoc, rank, sum
  13. [13] § Results › 40 Hz flicker increases low-alpha functional connectivity which is correlated with better behavior performance ↔ EEG Analysis/Functions/GenMatCode-main/Numeric/SignalProcessing/wpli_TrialIndex.m, the whole file · a weak match · score 0.55 · Weighted Phase Lag, volume conduction, WPLI, connectivity

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

MATLAB · 4,365 lines · 195 KB · no license · 5 matches

  1. %% Step 4 - Does PSD. (derived from Step2_WPLI by Lu Z)TrialType
  2. % Must add GenMatCode-main and all subfolders to path before running
  3. S4Start = tic;
  4. currDate = strrep(datestr(datetime), ':', '_');
  5. currDate = datestr(datetime, 'yy-mm-dd_HHMMSSFFF');
  6. scriptName = mfilename;
  7. %% Load All EEGs
  8. % LoadPath='Y:\singer\LuZhang\Project6-EEG\Results\Step0-PreparingData\';
  9. cd("2_CheckSync_Outputs\")
  10. %% Get newest file
  11. newestFolder = getNewestFolder();
  12. %cd(newestFolder);
  13. %newestEEGData = getNewestFile(); % Automatically gets the newest .mat file in folder
  14. %newestEEGData=
  15. %% Select file to load
  16. inputDatafile = '01-Apr-2024 13_38_57EEGDataStep2CheckSyncOutput.mat' %#ok<NOPTS> % can replace variable with hardcoded EEG file
  17. tic
  18. load(inputDatafile)
  19. %load('01-Apr-2024 13_38_57EEGDataStep2CheckSyncOutput.mat')
  20. %% Create folder
  21. disp(['EEG data loaded:', inputDatafile])
  22. toc
  23. % load('02-Feb-2024 17_46_26EEGData4EyeBlinkRemoval.mat')
  24. cd ..
  25. %cd ..
  26. SaveFolder=['4_COH\' currDate '_' scriptName];
  27. mkdir(SaveFolder);
  28. save([SaveFolder '\Step4DataInputSummary.mat'],'inputDatafile','currDate','scriptName')
  29. %% Create output text
  30. %Generate output message
  31. SummaryTextOutput='';
  32. % if isempty(SummaryTextOutput)
  33. % SummaryTextOutput = 'None';
  34. % end
  35. SummaryTextOutput = sprintf('Step4 started on %s\nScript used: %s \n', currDate, scriptName);
  36. SummaryTextOutput = sprintf('%s \nInput dataset: %s \n', SummaryTextOutput, inputDatafile);
  37. % SummaryTextOutput = sprintf('%s \n%d out of %d files synced:\n%s', SummaryTextOutput, iCsv, nCsv, subjectsList);
  38. % endDateAndTime = char(datetime('now','TimeZone','local','Format','d-MMM-y HH:mm:ss'));
  39. % endDuration = toc(tStart);
  40. % SummaryTextOutput = sprintf('%s\nStep1Syncing completed at: %s\nSyncing duration: %g seconds.', SummaryTextOutput, endDateAndTime, endDuration);
  41. % Add toc
  42. %Display the output message in the command window
  43. disp(SummaryTextOutput);
  44. %Create a .txt file with the output message
  45. filePath = fullfile(SaveFolder, [scriptName, currDate '_Summary.txt']);
  46. fileID = fopen(filePath, 'w');
  47. if fileID == -1
  48. error('Failed to open or create the file: %s', filePath);
  49. else
  50. fprintf(fileID, '%s', SummaryTextOutput);
  51. fclose(fileID);
  52. fprintf('File created and written: %s\n', filePath);
  53. end
  54. %% Set Channels, Locations for Heatmap
  55. % (Necessary if already done in Step 0?)
  56. channels = {'Fp1', 'AF3', 'F7', 'F3', 'FC1', 'FC5', 'T7', 'C3', 'CP1', 'CP5', 'P7', 'P3', 'Pz', 'PO3', 'O1', 'Oz', 'O2', 'PO4', 'P4', 'P8', 'CP6', 'CP2', 'C4', 'T8', 'FC6', 'FC2', 'F4', 'F8', 'AF4', 'Fp2', 'Fz', 'Cz', 'EXG1', 'EXG2', 'EXG3', 'EXG4', 'EXG5', 'EXG6', 'EXG7', 'EXG8'};
  57. EEGchInd=1:32;
  58. EEGch=channels(EEGchInd);
  59. ChNTotal=length(channels);
  60. NeedFields={'labels','theta','radius'};
  61. tempN=fieldnames(EEGList{1}.chanlocs);
  62. tempN1=tempN;
  63. tempN2=tempN;
  64. tempN1([1 2 11 12])=[];
  65. NeedI=setdiff(1:length(tempN),[1 2 12]);
  66. tempN2([2 12])=[];
  67. ChanPos=FieldName2Struct(tempN1);
  68. TempTable=cell2table(squeeze(struct2cell(EEGList{1}.chanlocs)));
  69. ChInd=table2array(cell2table(table2array(TempTable(11,:))));
  70. IndEx=find(ChInd>=33);
  71. ChIndWritten=[];
  72. DataTemp=zeros(length(tempN1),length(EEGchInd))+nan;
  73. for iFile=1:length(EEGList)
  74. TempTable=cell2table(squeeze(struct2cell(EEGList{iFile}.chanlocs)));
  75. ChI=table2array(cell2table(table2array(TempTable(11,:))));
  76. IndEx=find(ChI>=33);
  77. Invalid=find(ChI>=33);
  78. ChIndV=ChI;
  79. ChIndV(Invalid)=[];
  80. TempTable(:,Invalid)=[];
  81. InforAdd=TempTable;
  82. InforAdd([1 2 11 12],:)=[];
  83. [ChAdd,I1]=setdiff(ChIndV,ChIndWritten);
  84. temp=table2array(cell2table(table2array(InforAdd)));
  85. DataTemp(:,ChAdd)=temp(:,I1);
  86. ChIndWritten=find(isnan(DataTemp(1,:))==0);
  87. if length(ChIndWritten)==length(EEGchInd)
  88. break
  89. end
  90. end
  91. for ifield=1:length(tempN1)
  92. ChanPos=setfield(ChanPos,{1,1},tempN1{ifield},DataTemp(ifield,:));
  93. end
  94. ChanPos.labels=EEGch;
  95. ChanEEGLab=rmfield(ChanPos,'labels');
  96. tempName=fieldnames(ChanEEGLab);
  97. for iCh=1:length(EEGchInd)
  98. for iN=1:length(tempName)
  99. dataTemp=getfield(ChanPos,{1},tempName{iN});
  100. ChanEEGLab=setfield(ChanEEGLab,{1,iCh},tempName{iN},dataTemp(iCh));
  101. end
  102. ChanEEGLab(iCh).labels=EEGch{iCh};
  103. ChanEEGLab(iCh).urchan=iCh;
  104. end
  105. %% Assign Condition Groups *** Run to set GroupName
  106. FlickerSubj{4}=[20:27 30:33 44:51 53 55:58 62:64 66:69 71:73 75 77:81 83:87]; % All control subjects (Random and Light) no cuts
  107. GroupName{4}='BothControls'; % 'BothControls' = Random and Light together. previously 'Random' or 'Light'
  108. FlickerSubj{5}=[20:27 30:33 44:51 53 55:58 62:63 66 68 69 71:73 75 77:81 83:87]; % All control subjects minus 64 and 67 (for RT)
  109. GroupName{5}='BothControlsRT'; % 'BothControls' = Random and Light together. previously 'Random' or 'Light'
  110. FlickerSubj{1}=[10:19 34:42 52 54 59:61 65 70 74 76 82]; % SubjID of 40Hz group.
  111. GroupName{1}='40Hz';
  112. FlickerSubj{2} = [62:64 66:69 71:73 75 77:81 83:87]; % All SubjID for Light Group (no cuts)
  113. GroupName{2}='Light';
  114. FlickerSubj{3}=[20:27 30:33 44:51 53 55:58]; % SubjID of Random flicker group
  115. GroupName{3} = 'Random';
  116. FlickerSubj{6} = [62:63 66 68 69 71:73 75 77:81 83:87]; % All SubjID for Light Group for RT (removed s064 and s067 due to incorrect average RTs
  117. GroupName{6}='LightRT'; % No 64 67 for WPLI vs RT
  118. % Updated Groups (includes 2023 EEGs)
  119. % FlickerSubj{1}=[20:28 30:33 45 48:50 56:58]; % (20 total) SubjID of Random flicker group (Accuracy cuts)
  120. % FlickerSubj{2}=[10:19 39 40 52 54 59 60 74 76 82]; % (20 total) % SubjID of 40Hz group. Newest 40 Hz EEGs added (2023) (Accuracy Cut
  121. % FlickerSubj{3}=[]; % check if ID is still matching
  122. %% Get Accuracies for all subjects
  123. Acc=zeros(length(FileStruct),1)+nan;
  124. SubjsAvgRT = zeros(length(FileStruct),1)+nan;
  125. for iFile=1:length(FileStruct)
  126. if ~isempty(FileStruct(iFile).Subj)
  127. %% May need to uncomment out the below
  128. SubjID(iFile) = str2num(FileStruct(iFile).Subj(end-1:end));
  129. Acc(iFile)=length(FileStruct(iFile).hits)/(length(FileStruct(iFile).hits)+length(FileStruct(iFile).misses));
  130. %% Calculate Avg RT per subject using data found in FileStruct->dotsynch
  131. isHit = (FileStruct(iFile).dotsynch(:,3)==1); % get logical index for all trials that are Hits (misses and premature hits will = 0)
  132. colorchangeTimes = FileStruct(iFile).dotsynch(isHit,1); % get time of color change for Hit trials only
  133. subjRTtimes = FileStruct(iFile).dotsynch(isHit,2); % get RT time for Hit trials (this should ignore premature hits)
  134. subjRTduration = subjRTtimes - colorchangeTimes; % subtract RT time from color change time to get duration of RT (AKA the reaction time)
  135. subjRTdurInSecs = subjRTduration/512; % sample rate is generally 512 samples per second
  136. SubjsAvgRT(iFile) = mean(subjRTdurInSecs);
  137. end
  138. end
  139. %% SubjG assignment for WPLI calculation
  140. SubjG{1}=[];
  141. SubjG{2}=[];
  142. SubjG{3}=[];
  143. SubjG{4}=[];
  144. SubjG{5}=[];
  145. SubjG{6}=[];
  146. [~,SubjG{1},~]=intersect(SubjID,FlickerSubj{1}); % 40Hz
  147. [~,SubjG{2},~]=intersect(SubjID,FlickerSubj{2}); % Light
  148. [~,SubjG{3},~]=intersect(SubjID,FlickerSubj{3}); % Random
  149. [~,SubjG{4},~]=intersect(SubjID,FlickerSubj{4}); % Both Controls
  150. [~,SubjG{5},~]=intersect(SubjID,FlickerSubj{5}); % BothControlsRT
  151. [~,SubjG{6},~]=intersect(SubjID,FlickerSubj{6}); % LightRT
  152. %% Set Trial Groups (Hit, Miss, Etc.)
  153. ChCount=[];
  154. ChIndList={};
  155. TrialType{1}=1; %%%%%Hit trial
  156. TrialType{2}=0; %%%%%Miss trial
  157. TrialType{3}=[0 1]; %%%%%All trial
  158. % TrialType{4}=-1; %%%%%Premature trial
  159. TrialTypeName{1}='Hit'; %%%%%Hit trial
  160. TrialTypeName{2}='Miss'; %%%%%Miss trial
  161. TrialTypeName{3}='HitAndMiss'; %%%%%All trial
  162. % TrialTypeName{4}='Premature'; %%%%%Premature trial
  163. %% PSD parameters
  164. psdParameter.Fs=512;
  165. psdParameter.window=1024; % can increase window to 1024 from 512 8/18/24 - less smoth, but higher res
  166. psdParameter.noverlap = psdParameter.window/2; % can change overlap to 256 or half window size
  167. psdParameter.nfft=1024; % decrease nfft from 1024 to 512 - less smooth but decrease resolution
  168. nFre=psdParameter.nfft/2+1; % modified pwelch see Lu genmat code on git
  169. %% Epoch, Samp Rate, Make Save Trial Folder
  170. % (When is data epoched? What if data is already epoched from preprocessing)
  171. DataTimeRange=[-4 1]; %%%4 seconds before color-change and 1s after.
  172. AnaRange=[-4 0]; %%4s before color-change
  173. % AnaRange=[0 1]; %%1s after color change
  174. SampRate=512;
  175. SampI=(AnaRange-DataTimeRange(1))*SampRate;
  176. SampI=SampI(1)+1:SampI(2);
  177. % clear CohGroup
  178. parfor iFile=1:length(EEGList)
  179. ChTempN(iFile)=size(EEGList{iFile}.AllChData,1);
  180. end
  181. SavePath = ['4_COH\' currDate '_' scriptName '\'];
  182. if ~exist(SavePath, 'dir')
  183. mkdir(SavePath)
  184. end
  185. %% Get coherence for all pairs of channels PER SUBJECT!
  186. % This takes a looong time! (Data from here goes in TrialCrossSpec)
  187. % SaveTrialSubj='Y:\singer\LuZhang\Project6-EEG\Results\Step2-COH\TrialCrossSpec\';
  188. SaveTrialSubj=['4_COH\' currDate '_' scriptName '\TrialCrossSpec\'];
  189. mkdir(SaveTrialSubj);
  190. parpool(12)
  191. tic
  192. for iFile=1:length(EEGList)
  193. SaveTemp=[SaveTrialSubj EEGList{iFile}.filename(1:4) '\'];
  194. mkdir(SaveTemp);
  195. if ~isempty(EEGList{iFile})
  196. ChTempN=size(EEGList{iFile}.AllChData,1);
  197. TrialTypeTemp=EEGList{iFile}.TrialType;
  198. TrialI=[];
  199. %% Calculate for pair of channels
  200. tic
  201. for iCh=1:ChTempN
  202. for jCh=iCh:ChTempN
  203. clear TempTrial1 tempSig1 TempTrial2 tempSig2 TrialSpec
  204. tempSig1=squeeze(EEGList{iFile}.AllChData(iCh,SampI,:));
  205. tempSig2=squeeze(EEGList{iFile}.AllChData(jCh,SampI,:));
  206. if sum(sum(isnan(tempSig1)))>1||sum(sum(isnan(tempSig2)))>1
  207. continue
  208. end
  209. % tic
  210. for iTrial=1:length(TrialTypeTemp) % loop by trial
  211. if length(TrialTypeTemp)>1
  212. TempTrial1(iTrial).Data=tempSig1(:,iTrial);
  213. TempTrial2(iTrial).Data=tempSig2(:,iTrial);
  214. else
  215. TempTrial1(iTrial).Data=tempSig1;
  216. TempTrial2(iTrial).Data=tempSig2;
  217. end
  218. TempTrial1(iTrial).Time=([1:length(TempTrial1(iTrial).Data)]-1)/512;
  219. TempTrial2(iTrial).Time=([1:length(TempTrial2(iTrial).Data)]-1)/512;
  220. end
  221. % clear TrialSpec
  222. %%%Old version to calculate CrossSpec,tested equal to
  223. %%%new version
  224. % [TrialSpec1.Sxy,TrialSpec1.Sxx,TrialSpec1.Syy,TrialSpec1.w,TrialSpec1.options,ValidIndex]=coh_TrialData(TempTrial1,TempTrial2,psdParameter);
  225. %%%Old version to calculate CrossSpec,tested equal to
  226. % %%%new version
  227. % psdParameter.noverlap=500;
  228. % psdParameter.nfft=512;
  229. % psdParameter.window=512;
  230. % nFre=psdParameter.nfft/2+1;
  231. [TrialSpec.Sxy,TrialSpec.Sxx,TrialSpec.Syy,TrialSpec.w,TrialSpec.options,ValidIndex]=crossspec_EqualTriL(TempTrial1,TempTrial2,psdParameter); %% find in genmat code
  232. save([SaveTemp 'Ch' num2str(iCh) 'Ch' num2str(jCh) '.mat'],'TrialSpec','ValidIndex','psdParameter');
  233. % a=crossspec_Trial(TrialSpec);
  234. % figure;
  235. % plot(a.Fre,abs(((a.wpli))))
  236. % figure;
  237. % plot(a.Fre,abs(mean((a.wpli(1:30,:)))))
  238. % figure;
  239. % plot(a.Fre,abs(mean((a.wpli(1:30,:)))))
  240. % figure;
  241. % plot(TempTrial1(4).Data);hold on;plot(TempTrial2(4).Data,'r.')
  242. for iTrialType=1:length(TrialType) % group: hit, miss, etc.
  243. TrialI=[];
  244. parfor j=1:length(TrialType{iTrialType})
  245. TrialI=union(TrialI,find(TrialTypeTemp==TrialType{iTrialType}(j)));
  246. end
  247. TrialI=intersect(TrialI,ValidIndex);
  248. if ~isempty(TrialI)
  249. % CohGroup{iGroup,iFile}{iCh,jCh}=Coh_TrialIndex(TrialSpec1,TrialI);
  250. CohGroup{iTrialType,iFile}{iCh,jCh}=crossspec_TrialIndex(TrialSpec,TrialI);
  251. %% %confirmed Old and New version of Cross-Spectrum results in same coherence results.
  252. % D1=Coh_TrialIndex(TrialSpec1,TrialI);
  253. % D2=crossspec_TrialIndex(TrialSpec,TrialI);
  254. % figure;
  255. % plot(D1.Fre,(D1.Cxy));hold on;
  256. % plot(D2.Fre,(D2.Cxy),'r.');hold on;c
  257. %%%confirmed Old and New version of Cross-Spectrum results in same coherence results.
  258. end
  259. end
  260. % toc
  261. %% figure;
  262. % Temp=CohGroup{1,iFile}{iCh,jCh};
  263. % subplot(2,1,1)
  264. % plot(Temp.Fre,(Temp.Cxy));
  265. % subplot(2,1,2)
  266. %
  267. % plot(Temp.Fre,log(abs(Temp.Pxx)));
  268. % hold on;
  269. % plot(Temp.Fre,log(abs(Temp.Pyy)));
  270. end
  271. end
  272. toc
  273. end
  274. end
  275. toc
  276. % Check if everything above works ***
  277. %% Save workspace to be used for Step 5: WPLITrial Group. This step takes a long time, creates 20gb file!
  278. COHSaveFileName =['COHdata_forWPLITrialType_' currDate '_' scriptName '.mat'];
  279. COH_Save_Path = [SavePath COHSaveFileName];
  280. save(COH_Save_Path,'-v7.3') % Check this
  281. % load(COH_Save_Path) % loading should not be necessary as all variables in workspace should be in that save file
  282. %% For visualization - Set-up - Must run before plotting anything below
  283. % (of what? COH,WPLI and PSD?)
  284. load('chanPosColin27');
  285. %
  286. Fre=CohGroup{1,1}{1,2}.Fre; % may need to import CohGroup from a COHdata file
  287. FBand=[2 100]; % consider changing [1 100] to [2 100] due to normalization
  288. FreInd=find(Fre>=FBand(1)&Fre<=FBand(2));
  289. Fplot=Fre(FreInd);
  290. PSDall=zeros(length(FileStruct),length(FreInd),ChNTotal,length(TrialType))+nan;
  291. COHall=zeros(length(FileStruct),length(FreInd),ChNTotal,ChNTotal-1,length(TrialType))+nan;
  292. WPLIall=zeros(length(FileStruct),length(FreInd),ChNTotal,ChNTotal-1,length(TrialType))+nan;
  293. %
  294. for iTrialType=1:length(TrialType)
  295. for iFile=1:length(FileStruct)
  296. for iCh=1:size(CohGroup{iTrialType,iFile},1)
  297. for jCh=iCh+1:size(CohGroup{iTrialType,iFile},2)
  298. if isempty(CohGroup{iTrialType,iFile}{iCh,jCh})
  299. continue;
  300. end
  301. if iCh==1
  302. PSDall(iFile,:,iCh,iTrialType)=CohGroup{iTrialType,iFile}{iCh,jCh}.Pxx(FreInd);
  303. PSDall(iFile,:,jCh,iTrialType)=CohGroup{iTrialType,iFile}{iCh,jCh}.Pyy(FreInd);
  304. end
  305. COHall(iFile,:,iCh,jCh,iTrialType)=CohGroup{iTrialType,iFile}{iCh,jCh}.Cxy(FreInd);
  306. WPLIall(iFile,:,iCh,jCh,iTrialType)=CohGroup{iTrialType,iFile}{iCh,jCh}.wpli(FreInd);
  307. end
  308. end
  309. end
  310. end
  311. %% Peak Alpha WPLI Distribution Histogram
  312. % Get all alpha values
  313. alphaLowerLimitFreqHz = 8;
  314. alphaUpperLimitFreqHz = 13;
  315. % Find indices of values in Fre between 8 and 13 (inclusive)
  316. allAlphaFreIndices = find(Fre >= alphaLowerLimitFreqHz & Fre <= alphaUpperLimitFreqHz);
  317. HitMissTrialType = 3;
  318. WPLIallAlpha = squeeze(WPLIall(:,allAlphaFreIndices,1:32,1:32,HitMissTrialType)); % size(WPLIall) ans = 67 197 40 39 3- 67subs x allFreqs x iCh x jCh x TrialType
  319. % Find the peak value within alpha (dimension 2)
  320. [peakAlphaValues, peakAlphaIndices] = max(WPLIallAlpha, [], 2);
  321. % Reshape the result to 3D
  322. PeakWPLIallAlpha = squeeze(peakAlphaValues);
  323. %% Plot distribution of PeakWPLIallAlpha
  324. % Flatten PeakWPLIallAlpha to 1D for distribution analysis
  325. PeakWPLIallAlphaFlat = PeakWPLIallAlpha(:);
  326. % Remove NaN values
  327. PeakWPLIallAlphaFlat = PeakWPLIallAlphaFlat(~isnan(PeakWPLIallAlphaFlat));
  328. % Calculate the total number of data points (channel pairs)
  329. nDataPoints = numel(PeakWPLIallAlphaFlat);
  330. % Compute top percentiles
  331. top25Percent = prctile(PeakWPLIallAlphaFlat, 75);
  332. top10Percent = prctile(PeakWPLIallAlphaFlat, 90);
  333. top5Percent = prctile(PeakWPLIallAlphaFlat, 95);
  334. top1Percent = prctile(PeakWPLIallAlphaFlat, 99);
  335. top0_1Percent = prctile(PeakWPLIallAlphaFlat, 99.9);
  336. top0_01Percent = prctile(PeakWPLIallAlphaFlat, 99.99);
  337. % Count data points greater than each percentile
  338. countAbove25Percent = sum(PeakWPLIallAlphaFlat > top25Percent);
  339. countAbove10Percent = sum(PeakWPLIallAlphaFlat > top10Percent);
  340. countAbove5Percent = sum(PeakWPLIallAlphaFlat > top5Percent);
  341. countAbove1Percent = sum(PeakWPLIallAlphaFlat > top1Percent);
  342. countAbove0_1Percent = sum(PeakWPLIallAlphaFlat > top0_1Percent);
  343. countAbove0_01Percent = sum(PeakWPLIallAlphaFlat > top0_01Percent);
  344. % Display the results
  345. fprintf('Top 25%% WPLI value: %.4f, Data points above: %d\n', top25Percent, countAbove25Percent);
  346. fprintf('Top 10%% WPLI value: %.4f, Data points above: %d\n', top10Percent, countAbove10Percent);
  347. fprintf('Top 5%% WPLI value: %.4f, Data points above: %d\n', top5Percent, countAbove5Percent);
  348. fprintf('Top 1%% WPLI value: %.4f, Data points above: %d\n', top1Percent, countAbove1Percent);
  349. fprintf('Top 0.1%% WPLI value: %.4f, Data points above: %d\n', top0_1Percent, countAbove0_1Percent);
  350. fprintf('Top 0.01%% WPLI value: %.4f, Data points above: %d\n', top0_01Percent, countAbove0_01Percent);
  351. % Plot the histogram of PeakWPLIallAlpha
  352. figure;
  353. histogram(PeakWPLIallAlphaFlat, 'Normalization', 'probability', 'BinWidth', 0.02);
  354. hold on;
  355. % Add vertical lines for top percentiles
  356. xline(top25Percent, '--k', 'Top 25%', 'LineWidth', 1.5);
  357. xline(top10Percent, '--c', 'Top 10%', 'LineWidth', 1.5);
  358. xline(top5Percent, '--r', 'Top 5%', 'LineWidth', 1.5);
  359. xline(top1Percent, '--g', 'Top 1%', 'LineWidth', 1.5);
  360. xline(top0_1Percent, '--b', 'Top 0.1%', 'LineWidth', 1.5);
  361. xline(top0_01Percent, '--m', 'Top 0.01%', 'LineWidth', 1.5);
  362. % Label the axes
  363. xlabel('WPLI');
  364. ylabel('Fraction');
  365. % Add a title including the total number of data points
  366. title(['Distribution of Peak Alpha WPLI (Total data points: ', num2str(nDataPoints), ')']);
  367. % Improve plot appearance
  368. grid on;
  369. % % Display the size of the resulting array
  370. % disp('Size of the resulting 3D array:');
  371. % disp(size(PeakWPLIallAlpha));
  372. %% Plot the histogram of PeakWPLIallAlpha (SuppFig4A MS version)
  373. % May need to run the previous section first! MKA 2025-03-18
  374. figure;
  375. histogram(PeakWPLIallAlphaFlat, 'Normalization', 'probability', 'BinWidth', 0.02);
  376. hold on;
  377. % Add vertical lines for top percentiles
  378. xline(top25Percent, '--k', 'Top 25%', 'LineWidth', 1.5);
  379. % Label the axes
  380. xlabel('WPLI');
  381. ylabel('Fraction');
  382. % Add a title including the total number of data points
  383. title(['Distribution of Peak Alpha WPLI (Total data points: ', num2str(nDataPoints), ')']);
  384. % Improve plot appearance
  385. grid on;
  386. saveas(gcf, fullfile('Fig4Panels', 'SuppFig4A_PeakWPLI_Alpha.svg'));
  387. %% Peak Alpha WPLI 40Hz vs Light ALL CHANNELS (not used) - search "stats preceding fig4D" for signif channels only
  388. % Find the peak WPLI value within alpha (dimension 2)
  389. [peakAlphaValues, peakAlphaIndices] = max(WPLIallAlpha, [], 2, "includemissing");
  390. alphaFreqHz = 8:.5:13;
  391. % Initialize a copy of peakAlphaIndices
  392. peakAlphaIndicesNaN = peakAlphaIndices;
  393. % If all alpha WPLI values are NaN, then set the max index to NaN
  394. peakAlphaIndicesNaN(all(isnan(WPLIallAlpha),2)) = NaN;
  395. % Reshape the result to 3D
  396. PeakAlphaIndicesNaN3D = squeeze(peakAlphaIndicesNaN);
  397. % % Check with single participant
  398. % singleparticiantAllChPairPeakAlpha = squeeze(PeakAlphaIndicesNaN3D(1,:,:))
  399. %
  400. % % Convert Indices to Correct Corresponding Frequnecy (Hz)
  401. % % Find valid indices (values between 1 and 11)
  402. % validIdx = singleparticiantAllChPairPeakAlpha >= 1 & singleparticiantAllChPairPeakAlpha <= 11;
  403. %
  404. % % Initialize the output array with NaN, preserving original shape
  405. % singleallchpfreqs = NaN(size(singleparticiantAllChPairPeakAlpha));
  406. %
  407. % % Perform mapping only for valid indices
  408. % singleallchpfreqs(validIdx) = alphaFreqHz(singleparticiantAllChPairPeakAlpha(validIdx));
  409. % Convert Indices (3D) to Correct Corresponding Frequncy (Hz)
  410. validIdx = PeakAlphaIndicesNaN3D >= 1 & PeakAlphaIndicesNaN3D <= 11;
  411. % Initialize the output array with NaN, preserving original shape
  412. PeakAlphaFreqs = NaN(size(PeakAlphaIndicesNaN3D));
  413. % Perform mapping only for valid indices
  414. PeakAlphaFreqs(validIdx) = alphaFreqHz(PeakAlphaIndicesNaN3D(validIdx));
  415. % singleallchpfreqs = alphaFreqHz(singleparticiantAllChPairPeakAlpha)
  416. GroupSubjs_40Hz = 1; % 40Hz group
  417. GroupSubjs_Light = 6; % LightRT group
  418. % % Extract Peak Alpha WPLI into groups: 40 Hz & Light
  419. % % PeakWPLIallAlpha: 67 subs x 32ch x 32ch
  420. % PeakAlphaWPLI_40HzGroup = PeakWPLIallAlpha(SubjG{GroupSubjs_40Hz},:,:);
  421. % PeakAlphaWPLI_LightGroup = PeakWPLIallAlpha(SubjG{GroupSubjs_Light},:,:);
  422. % % single particiapn
  423. % singledudeWPLIallpair = squeeze(PeakWPLIallAlpha(1,:,:))
  424. % Extract Peak Alpha WPLI freq into groups: 40 Hz & Light
  425. % PeakWPLIallAlpha: 67 subs x 32ch x 32ch
  426. PeakAlphaFreqs_40HzGroup = PeakAlphaFreqs(SubjG{GroupSubjs_40Hz},:,:);
  427. PeakAlphaFreqs_LightGroup = PeakAlphaFreqs(SubjG{GroupSubjs_Light},:,:);
  428. % Get average peak alpha WPLI frequncy for each group for all channel pairs
  429. MeanPeakAlphaFreq_40Hz = squeeze(mean(PeakAlphaFreqs_40HzGroup, 1));
  430. MeanPeakAlphaFreq_Light = squeeze(mean(PeakAlphaFreqs_LightGroup, 1));
  431. % Display average peak alpha frequency tables (32x32)
  432. ChannelLabels = {'Fp1', 'AF3', 'F7', 'F3', 'FC1', 'FC5', 'T7', 'C3', 'CP1', 'CP5', ...
  433. 'P7', 'P3', 'Pz', 'PO3', 'O1', 'Oz', 'O2', 'PO4', 'P4', 'P8', ...
  434. 'CP6', 'CP2', 'C4', 'T8', 'FC6', 'FC2', 'F4', 'F8', 'AF4', 'Fp2', 'Fz', 'Cz'};
  435. fprintf('\nMean Peak Alpha Frequency (40Hz Group):\n');
  436. disp(array2table(MeanPeakAlphaFreq_40Hz, 'VariableNames', ChannelLabels, 'RowNames', ChannelLabels));
  437. fprintf('\nMean Peak Alpha Frequency (Light Group):\n');
  438. disp(array2table(MeanPeakAlphaFreq_Light, 'VariableNames', ChannelLabels, 'RowNames', ChannelLabels));
  439. % Flatten both 32 x 32 arrays into two 1D-arrays (492 channel pairs each)
  440. UpperTriIdx = find(triu(ones(32, 32), 1)); % Indices of upper triangular elements
  441. FlattenedAlphaFreq_40Hz = MeanPeakAlphaFreq_40Hz(UpperTriIdx);
  442. FlattenedAlphaFreq_Light = MeanPeakAlphaFreq_Light(UpperTriIdx);
  443. % Compute average peak alpha frequency across all channel pairs for each group
  444. AvgPeakAlphaFreq_40Hz = mean(FlattenedAlphaFreq_40Hz);
  445. AvgPeakAlphaFreq_Light = mean(FlattenedAlphaFreq_Light);
  446. fprintf('\nAverage Peak Alpha Frequency Across All Channel Pairs:\n');
  447. fprintf('40Hz Group: %.4f Hz\n', AvgPeakAlphaFreq_40Hz);
  448. fprintf('Light Group: %.4f Hz\n', AvgPeakAlphaFreq_Light);
  449. % Test for normality (to decide if t-test or ranksum)
  450. [H_40Hz, p_40Hz] = kstest(FlattenedAlphaFreq_40Hz);
  451. [H_Light, p_Light] = kstest(FlattenedAlphaFreq_Light);
  452. fprintf('\nNormality Test Results:\n');
  453. fprintf('40Hz Group: H = %d, p = %.16f\n', H_40Hz, p_40Hz);
  454. fprintf('Light Group: H = %d, p = %.16f\n', H_Light, p_Light);
  455. % Decide on statistical test
  456. if H_40Hz == 0 && H_Light == 0
  457. % Normally distributed: Use independent t-test
  458. [h_ttest, p_ttest] = ttest2(FlattenedAlphaFreq_40Hz, FlattenedAlphaFreq_Light);
  459. test_used = 't-test';
  460. p_value = p_ttest;
  461. else
  462. % Non-normally distributed: Use Wilcoxon rank-sum test
  463. [p_ranksum, h_ranksum] = ranksum(FlattenedAlphaFreq_40Hz, FlattenedAlphaFreq_Light);
  464. test_used = 'Wilcoxon rank-sum test';
  465. p_value = p_ranksum;
  466. end
  467. % Display results
  468. fprintf('Statistical Test Used: %s\n', test_used);
  469. fprintf('p-value: %.16f\n', p_value);
  470. % Perform one-sided Wilcoxon rank-sum test
  471. [p_ranksum_right, h_ranksum_right] = ranksum(FlattenedAlphaFreq_40Hz, FlattenedAlphaFreq_Light, 'tail', 'right');
  472. [p_ranksum_left, h_ranksum_left] = ranksum(FlattenedAlphaFreq_40Hz, FlattenedAlphaFreq_Light, 'tail', 'left');
  473. % Display results
  474. fprintf('\nOne-Sided Wilcoxon Rank-Sum Test Results:\n');
  475. fprintf('H0: 40Hz <= Light | p-value (right-tailed, 40Hz > Light): %.16f\n', p_ranksum_right);
  476. fprintf('H0: 40Hz >= Light | p-value (left-tailed, 40Hz < Light): %.16f\n', p_ranksum_left);
  477. % Create Violin plots showing distribution
  478. % Create figure
  479. figure;
  480. hold on;
  481. % Combine data for violin plot
  482. groupLabels = [repmat({'40Hz'}, length(FlattenedAlphaFreq_40Hz), 1); ...
  483. repmat({'Light'}, length(FlattenedAlphaFreq_Light), 1)];
  484. data = [FlattenedAlphaFreq_40Hz; FlattenedAlphaFreq_Light];
  485. % Create violin plot
  486. violinplot(data, groupLabels);
  487. % Format plot
  488. title('Violin Plot of Peak Alpha Frequency');
  489. ylabel('Peak Alpha Frequency (Hz)');
  490. xlabel('Group');
  491. ylim([8 13]); % Set y-axis range from 8 Hz to 13 Hz
  492. grid on;
  493. hold off;
  494. %% PSD related parameters
  495. LogPSDraw=log(abs(PSDall));
  496. NoiseInd=find(Fplot>=58&Fplot<=62);
  497. NormBandI=setdiff(1:length(Fplot),NoiseInd);
  498. PSDall=PSDall./repmat(nansum(PSDall(:,NormBandI,:,:),2),1,length(Fplot),1,1);
  499. LogPSD=log(abs(PSDall));
  500. PlotColor2=[1 0 0;0 0 1];
  501. % ParamPSD.ANOVAstats='Anova';
  502. ParamPSD.PlotType=3;
  503. ParamPSD.SigPlot='Anova';
  504. ParamPSD.SigPlot='Ttest';
  505. ParamPSD.CorrName='fdr'; %%%methold for multi-compairson
  506. ParamPSD.Q=0.1;
  507. ParamPSD.Ytick=[0 0.002 0.004];
  508. ParamPSD.LegendShow=0;
  509. ParamPSD.Legend=[];
  510. ParamPSD.TimeRepeatAnova=1;
  511. ParamPSD.GroupRepeatAnova=0;
  512. ParamPSD.RepeatAnova=0;
  513. ParamPSD.TimeCol=Fplot;
  514. ParamPSD.Paired=1;
  515. ParamPSD.BinName='Fre';
  516. ParamPSD.Bin=Fplot;
  517. ParamPSD.TimeComparison=0;
  518. ParamPSD.statisP=1; % 1 to do stats and plot. 0 will do stats, but not plot, will be faster. Uses R. If error, set as 0
  519. ParamPSD.Ytick=[-8:4:0];
  520. ParamPSD.Crit_p=0.05;
  521. % One color for ea of the 3 groups. Blue for 40, Gold/yellow/orange for Light, Red for Random
  522. FlickerColor=[31 125 184; 219 129 50; 150 27 27]/255;
  523. FlickerColor=[31 125 184; 150 27 27; 219 129 50]/255; %40, Random, Light
  524. % FlickerColor=[0.5 0.5 0.5;0.9 0.1 0.3];
  525. load('Functions\GenMatCode-main\Plotfun\Color\colorMapPN.mat')
  526. load('Functions\GenMatCode-main\Plotfun\Color\colorMapPNraw.mat')
  527. %% WPLI - Weight Phase Lag Index - Parameters
  528. ParamWPLI=ParamPSD;
  529. ParamWPLI.Ytick= [0:0.1:0.2]; %#ok<NBRAK2>
  530. SubSaveWPLI=[SavePath 'WPLI\'];
  531. ParamWPLI.SigPlot='Anova';
  532. mkdir(SubSaveWPLI)
  533. ParamWPLI.statisP=1;
  534. P.xLeft=0.01; %%%%%%Left Margin
  535. P.xRight=0.01; %%%%%%Right Margin
  536. P.yTop=0.01; %%%%%%Top Margin
  537. P.yBottom=0.01; %%%%%%Bottom Margin
  538. P.xInt=0.005; %%%%%%Width-interval between subplots
  539. P.yInt=0.005; %%%%%%Height-interval between subplots
  540. %% WPLI - Weight Phase Lag Index - Calculation *** Fig4b - this takes a long time
  541. WPLIStimGroupIndices = [1 3 6]; % 1=40, 3=Random, 6=LightRT
  542. SubjGWPLI = SubjG(WPLIStimGroupIndices); % Which three groups to include
  543. GroupNameWPLI = GroupName(WPLIStimGroupIndices);
  544. todayDate = datestr(now, 'yymmdd');
  545. for iTrialType=3%1:length(TrialType)
  546. SaveTemp=[SubSaveWPLI todayDate '\' TrialTypeName{iTrialType} '\'];
  547. mkdir(SaveTemp)
  548. SubSaveFig=[SaveTemp 'Chan\'];
  549. mkdir(SubSaveFig)
  550. % CH-Ch WPLI plot.tif figure;
  551. iPlot=0;
  552. alphaPeakAmplitudeList = zeros(length(EEGchInd),length(EEGchInd),length(SubjGWPLI));
  553. alphaPeakFrequencyList = zeros(length(EEGchInd),length(EEGchInd),length(SubjGWPLI));
  554. % alphaPeakAmplitudeListEmpty = double.empty(length(EEGchInd),length(EEGchInd),length(SubjG),0)
  555. for iCh=1:length(EEGchInd)
  556. for jCh=iCh+1:length(EEGchInd)
  557. clear DataPlot
  558. for iStimGroup=1:length(SubjGWPLI)
  559. DataPlot{iStimGroup}= squeeze(WPLIall(SubjGWPLI{iStimGroup},:,EEGchInd(iCh),EEGchInd(jCh),iTrialType));
  560. Invalid=isnan(DataPlot{iStimGroup}(:,1));
  561. DataPlot{iStimGroup}(Invalid,:)=[];
  562. end
  563. if isempty(DataPlot{1})||isempty(DataPlot{2})
  564. continue;
  565. end
  566. iPlot=iPlot+1;
  567. subplotLU(length(EEGchInd),length(EEGchInd),iCh,jCh,P);
  568. ParamWPLI.PathSave=[SaveTemp 'Light40Rand' EEGch{iCh} '-' EEGch{jCh}];
  569. % [~,COHComStatis{iCom,iCh,jCh}]=RateHist_GroupPlot(Fplot,DataPlot,FlickerColor,ParamCOH);
  570. tic
  571. RateHist_GroupPlot(Fplot,DataPlot,FlickerColor,ParamWPLI);
  572. % toc
  573. text(50,ParamWPLI.Ytick(end),[EEGch{iCh} '-' EEGch{jCh}]);
  574. set(gca,'xlim',FBand,'xtick',[0:20:120],'xticklabel',[],'yticklabel',[]);
  575. end
  576. end
  577. %% Permutation to set WPLI threshold ***fig4bc
  578. % Create "significant" WPLI threhold curve using permutation of current channel pair data.
  579. % -MKA 2024-12-11
  580. %Set up save folder
  581. saveDate = datestr(datetime, 'yy-mm-dd_HHMMSSFFF');
  582. SaveTemp=[SaveTrialSubj 'WPLIPermutation_' saveDate '\'];
  583. mkdir(SaveTemp);
  584. % Define parameters
  585. nPermutations = 10000; % Number of permutations, 10k or 1mil
  586. nSubjs = length(SubjID); % total number of subjects in analysis
  587. % Define the frequency range of interest: 2-55Hz and every half frequency inbetween
  588. % frequencies = [2:0.5:55]; % start at 1 or 2 hz? Cut off before 60 Hz
  589. % Pre-allocate storage for maximum WPLI curves across frequencies
  590. % WPLIperms = zeros(nPermutations, length(frequencies));
  591. % WPLIallPerm=zeros(length(FileStruct),length(FreInd),ChNTotal,ChNTotal-1,length(TrialType))+nan;
  592. %% Preallocate CohGroupPerm as a cell array
  593. CohGroupPerm = cell(3, nPermutations);
  594. % Begin permutations
  595. tStartPerm = tic;
  596. for iPermutation = 1:nPermutations
  597. % Step 1: Randomly select two subjects and a channel pair
  598. randSubjs = randperm(nSubjs, 2); %randomly choose a random Subject A and Subject B from all three groups.
  599. SubAInd = randSubjs(1);
  600. SubAData = EEGList{1,SubAInd};
  601. SubBInd = randSubjs(2);
  602. SubBData = EEGList{1,SubBInd};
  603. randChs = randperm(32,2); %randomly choose a chan Chi from subject i,  Chj for subject j
  604. Chi = randChs(1);
  605. Chj = randChs(2);
  606. % Step 2: Find the lower number of trials between the two selected subjects
  607. nTrials_SubA = SubAData.trials;
  608. nTrials_SubB = SubBData.trials;
  609. minTrials = min(nTrials_SubA, nTrials_SubB);
  610. minTrialsIndex = 1:minTrials;
  611. % Step 3: Get `minTrials` from each subject
  612. clear TempTrial1 DataA TempTrial2 DataB TrialSpec
  613. DataA=squeeze(SubAData.AllChData(Chi,SampI,minTrialsIndex)); % analog to "tempSigA"
  614. DataB=squeeze(SubBData.AllChData(Chj,SampI,minTrialsIndex)); % analog to "tempSigB"
  615. %% Step 4: Calculate WPLI for selected trials
  616. % WPLIperms(iPermuation, :) = calculateWPLI(DataA, DataB);
  617. % Preallocate TempTrial1 and TempTrial2 as structure arrays "FOR SPEED"
  618. TempTrial1(minTrials).Data = []; % Preallocate Data field
  619. TempTrial1(minTrials).Time = []; % Preallocate Time field
  620. TempTrial2(minTrials).Data = []; % Preallocate Data field
  621. TempTrial2(minTrials).Time = []; % Preallocate Time field
  622. for iTrial=1:minTrials % loop by trial
  623. if minTrials>1
  624. TempTrial1(iTrial).Data=DataA(:,iTrial);
  625. TempTrial2(iTrial).Data=DataB(:,iTrial);
  626. else
  627. TempTrial1(iTrial).Data=DataA;
  628. TempTrial2(iTrial).Data=DataB;
  629. end
  630. TempTrial1(iTrial).Time=([1:length(TempTrial1(iTrial).Data)]-1)/512;
  631. TempTrial2(iTrial).Time=([1:length(TempTrial2(iTrial).Data)]-1)/512;
  632. end
  633. [TrialSpec.Sxy,TrialSpec.Sxx,TrialSpec.Syy,TrialSpec.w,TrialSpec.options,ValidIndex]=crossspec_EqualTriL(TempTrial1,TempTrial2,psdParameter); %% find in genmat code
  634. % save([SaveTemp 'Ch' num2str(Chi) 'Ch' num2str(Chj) '.mat'],'TrialSpec','ValidIndex','psdParameter');
  635. jTrialType = iTrialType;
  636. for jTrialType=3%1:length(TrialType) % group: hit, miss, etc.
  637. %iTrialType hardcoded to 3; 3=Hit&Miss (all trials)
  638. TrialI=[];
  639. parfor j=1:length(TrialType{jTrialType}) % starting pool takes time; does this need to be parfor?
  640. TrialI=union(TrialI,find(TrialTypeTemp==TrialType{jTrialType}(j)));
  641. end
  642. TrialI=intersect(TrialI,ValidIndex);
  643. if ~isempty(TrialI)
  644. % CohGroup{iGroup,iFile}{iCh,jCh}=Coh_TrialIndex(TrialSpec1,TrialI);
  645. CohGroupPerm{iPermutation}=crossspec_TrialIndex(TrialSpec,TrialI);
  646. % CohGroup -> {3x67 cell} -> {3 trialTypes x 67 subjects}
  647. end
  648. end
  649. end
  650. elapsedTime = toc(tStartPerm);
  651. disp(['Elapsed time for script: ', num2str(elapsedTime), ' seconds']);
  652. %Save Permuation variable
  653. filename = [nPermutations 'CohGroupPerm_', datestr(now, 'yyyy-mm-dd'), '.mat'];
  654. save(filename, 'CohGroupPerm');
  655. %% Get WPLI from coherence
  656. Fre=CohGroupPerm{1}.Fre; % may need to import CohGroup from a COHdata file
  657. FBand=[2 100]; % consider changing [1 100] to [2 100] due to normalization
  658. FreInd=find(Fre>=FBand(1)&Fre<=FBand(2));
  659. Fplot=Fre(FreInd);
  660. PSDallPerm=zeros(length(FileStruct),length(FreInd),ChNTotal,length(TrialType))+nan;
  661. COHallPerm=zeros(length(FileStruct),length(FreInd),ChNTotal,ChNTotal-1,length(TrialType))+nan;
  662. WPLIallPerm=zeros(length(FileStruct),length(FreInd),ChNTotal,ChNTotal-1,length(TrialType))+nan;
  663. for iPermutation = 1:nPermutations
  664. if isempty(CohGroupPerm{iPermutation})
  665. continue;
  666. end
  667. if iCh==1
  668. PSDallPerm(iPermutation,:,iCh)=CohGroupPerm{iPermutation}.Pxx(FreInd);
  669. PSDallPerm(iPermutation,:,jCh)=CohGroupPerm{iPermutation}.Pyy(FreInd);
  670. end
  671. COHallPerm(iPermutation,:,iCh,jCh)=CohGroupPerm{iPermutation}.Cxy(FreInd);
  672. WPLIallPerm(iPermutation,:,iCh,jCh)=CohGroupPerm{iPermutation}.wpli(FreInd);
  673. end
  674. end
  675. %% Step 5: Determine the significance threshold for each frequency
  676. % EXTRACT relevant data for wpli threshold
  677. HitMissDataset = CohGroupPerm(3,:);
  678. % Initialize an output cell array of the same size.
  679. extractedStructs = cell(1, numel(HitMissDataset));
  680. % Extract the (31, 32) cell data from each cell in CohGroupPerm.
  681. extractedStructs = cellfun(@(x) x{31, 32}, HitMissDataset, 'UniformOutput', false);
  682. % Convert the extracted 'wpli' data to a 10,000x513 double matrix.
  683. wpliMatrix = cell2mat(cellfun(@(x) x.wpli, extractedStructs, 'UniformOutput', false)');
  684. wpliFreqs = extractedStructs{1,1}.Fre;
  685. % FIND WPLI threshold (value_of_topk)
  686. % Sort each column of wpliMatrix in descending order
  687. sortedMatrix = sort(wpliMatrix, 1, 'descend');
  688. % Extract the top_kth greatest value from each column
  689. p_value_threshold = 0.01; % p-val Significance threshold analogus to alpha value of 0.01
  690. top_0_01 = ceil(nPermutations * p_value_threshold); % Top K for p-value cutoff
  691. value_of_top0_01 = sortedMatrix(top_0_01, :);
  692. p_value_threshold2 = 0.001; % p-val Significance threshold analogus to alpha value of 0.01
  693. top_0_001 = ceil(nPermutations * p_value_threshold2); % Top K for p-value cutoff
  694. value_of_top0_001 = sortedMatrix(top_0_001, :);
  695. p_value_threshold3 = 0.0001; % p-val Significance threshold analogus to alpha value of 0.01
  696. top_0_0001 = ceil(nPermutations * p_value_threshold3); % Top K for p-value cutoff
  697. value_of_top0_0001 = sortedMatrix(top_0_0001, :);
  698. %% Step 6: Plot permutation WPLI threshold curve
  699. % Plot wpliFreqs (x-axis) against value100 (y-axis)
  700. figure;
  701. hold on;
  702. upperFreqLimit = 55; % Upper frequency limit of plot in (hz). Otherwise plot will go to >250 Hz
  703. upperFreqIndex = find(wpliFreqs >= upperFreqLimit, 1, 'first');
  704. % Plot the first line
  705. plot(wpliFreqs, value_of_top0_01, 'LineWidth', 2, 'DisplayName', ...
  706. ['p = ' num2str(p_value_threshold) ', Top K = ' num2str(top_0_01)]);
  707. % Plot the second line
  708. plot(wpliFreqs, value_of_top0_001, 'LineWidth', 2, 'DisplayName', ...
  709. ['p = ' num2str(p_value_threshold2) ', Top K = ' num2str(top_0_001)]);
  710. % Plot the third line
  711. plot(wpliFreqs, value_of_top0_0001, 'LineWidth', 2, 'DisplayName', ...
  712. ['p = ' num2str(p_value_threshold3) ', Top K = ' num2str(top_0_0001)]);
  713. % Add a legend
  714. legend('show', 'Location', 'best');
  715. % Add axis labels and title
  716. xlabel('Frequency (Hz)', 'FontSize', 12);
  717. ylabel('WPLI', 'FontSize', 12);
  718. title([num2str(size(wpliMatrix, 1)) ' Permutations. Threshold Curves for Different p-Values'], 'FontSize', 14);
  719. % Set y-axis to start at zero
  720. xlim([0 55])
  721. % ylim([0, max(max(value_of_top0_01(1:upperFreqIndex)), max(value_of_top0_001(1:upperFreqIndex)), max(value_of_top0_0001(1:upperFreqIndex)))]);
  722. ylim([0, max([max(value_of_top0_01(1:upperFreqIndex)), ...
  723. max(value_of_top0_001(1:upperFreqIndex)), ...
  724. max(value_of_top0_0001(1:upperFreqIndex))])]);
  725. % plot(wpliFreqs(1:index), value_of_topk(1:index), 'LineWidth', 2);
  726. % plot(wpliFreqs(1:index), value_of_topk2(1:index), 'LineWidth', 2);
  727. %
  728. % legend()
  729. %
  730. % % Add axis labels and title
  731. % xlabel('Frequency (Hz)', 'FontSize', 12);
  732. % ylabel('WPLI', 'FontSize', 12);
  733. % title([num2str(nPermutations) ' permutations. Top ' num2str(top_k) ' WPLI value. pvalue threshold of ' num2str(p_value_threshold)], 'FontSize', 14);
  734. %
  735. % % Set y-axis to start at zero
  736. % ylim([0, max(value_of_topk(1:index))]);
  737. grid on; % improve readability of plot
  738. hold off;
  739. % wpli_threshold_curve = zeros(32, 32, length(frequencies));
  740. % for freq_idx = 1:length(frequencies)
  741. % % Sort WPLI values for the current frequency across all permutations
  742. % sorted_values = sort(perm_idx_WPLICurve(:, freq_idx), 'descend');
  743. %
  744. % % Find the WPLI value corresponding to the top `top_k` value
  745. % wpli_threshold_curve(iCh, jCh, freq_idx) = sorted_values(top_k);
  746. % end
  747. %% Real > Permutated WPLI
  748. % get real WPLI - from
  749. load('10kCohGroupPerm_2024-12-17.mat')
  750. upperFreqLimit = 55; % Upper frequency limit of plot in (hz). Otherwise plot will go to >250 Hz
  751. upperFreqIndex = find(wpliFreqs >= upperFreqLimit, 1, 'first');
  752. lowerFreqLimit = 2;
  753. lowerFreqIndex = find(wpliFreqs >= lowerFreqLimit, 1, 'first');
  754. averagePermutatedWPLI_0to55_top0_01 = mean(value_of_top0_01(1:upperFreqIndex), 'omitnan');
  755. averagePermutatedWPLI_0to55_top0_001 = mean(value_of_top0_001(1:upperFreqIndex), 'omitnan');
  756. averagePermutatedWPLI_0to55_top0_0001 = mean(value_of_top0_0001(1:upperFreqIndex), 'omitnan');
  757. averagePermutatedWPLI_2to55_top0_01 = mean(value_of_top0_01(lowerFreqIndex:upperFreqIndex), 'omitnan');
  758. averagePermutatedWPLI_2to55_top0_001 = mean(value_of_top0_001(lowerFreqIndex:upperFreqIndex), 'omitnan');
  759. averagePermutatedWPLI_2to55_top0_0001 = mean(value_of_top0_0001(lowerFreqIndex:upperFreqIndex), 'omitnan');
  760. permutatedWPLI = value_of_top0_001;
  761. %% Initialize the average WPLI storage
  762. AvgRealWPLI = cell(length(EEGchInd), length(EEGchInd)); % Cell array for all channel pairs
  763. % AvgRealWPLI is a 32x32 cell containing 3x197 doubles (average WPLI per stim group (3) per freq (197)
  764. % Loop through all channel pairs 32x32
  765. for iCh = 1:length(EEGchInd)
  766. for jCh = iCh+1:length(EEGchInd) % Avoid duplicates and diagonal (upper triangle)
  767. % Initialize a 2D array to store averages for this channel pair
  768. AvgRealWPLI{iCh, jCh} = zeros(length(SubjGWPLI), size(WPLIall, 2)); % Rows: StimGroups, Cols: Frequencies
  769. % Loop through stimulus groups
  770. for iStimGroup = 1:length(SubjGWPLI)
  771. % Extract WPLI data for this stimulus group and channel pair
  772. DataPlot = squeeze(WPLIall(SubjGWPLI{iStimGroup}, :, EEGchInd(iCh), EEGchInd(jCh), iTrialType));
  773. % Remove invalid rows containing NaN
  774. Invalid = isnan(DataPlot(:, 1));
  775. DataPlot(Invalid, :) = [];
  776. % Compute the average across all valid rows for this stim group
  777. if ~isempty(DataPlot) % Ensure there is data after removing NaN
  778. AvgRealWPLI{iCh, jCh}(iStimGroup, :) = mean(DataPlot, 1); % Average across rows
  779. end
  780. end
  781. end
  782. end
  783. % outcome:
  784. % AvgRealWPLI is a 32x32 cell containing 3x197 doubles (average WPLI per stim group (3) per freq (197)
  785. %% Peak WPLI Extraction
  786. %%% set up
  787. BOIAlphas=[1 4 8 8 10 13 30 39.5 43;4 8 13 10 13 30 37 41.5 100];
  788. BOI2HzAlphas=[2 4 8 8 10 13 30 39.5 43 55;4 8 13 10 13 30 37 41.5 55 100];
  789. BOI2Hz=[2 4 8 13 30 39.5 43 55;4 8 13 30 37 41.5 55 100]; % No LOWER or UPPER ALPHA
  790. BandName={'Delta','Theta','Alpha','LowAlpha','HighAlpha','Beta','Gamma-1','Gamma-E', 'Gamma-55','Gamma-2'};
  791. BandHzNameAlphas={'1-4 Hz','4-8 Hz','8-13Hz','8-10Hz','10-13Hz','13-30 Hz','30-37 Hz','39-41 Hz','43-55 Hz', '55-100 Hz'};
  792. BandHzName2HzAlphas={'2-4 Hz','4-8 Hz','8-13Hz','8-10Hz','10-13Hz','13-30 Hz','30-37 Hz','39-41 Hz','43-55 Hz', '55-100 Hz'};
  793. BandHzName2HzAlphas={'2-4 Hz','4-8 Hz','8-13Hz','13-30 Hz','30-37 Hz','39-41 Hz','43-55 Hz', '55-100 Hz'};
  794. lowerboundHz = 2;
  795. upperboundHZ = 39; % can change to 55 Hz
  796. %%%
  797. % $% Filter frequency bands up to upperbound Hz
  798. BOI_filtered = BOI2Hz(:, BOI2Hz(2, :) <= upperboundHZ);
  799. BandName_filtered = BandName(BOI2Hz(2, :) <= upperboundHZ);
  800. BandHzName_filtered = BandHzName2Hz(BOI2Hz(2, :) <= upperboundHZ);
  801. %
  802. % Preallocate results (cell array for flexibility)
  803. numBands = size(BOI_filtered, 2);
  804. [numChannels, ~] = size(AvgRealWPLI);
  805. peakWPLI = cell(numChannels, numChannels);
  806. % %% TESTING ONLY - Create BOIsFreqIndices - comment when fixed
  807. % % Example: Fre is the frequency vector (197x1 double), BOI_filtered contains the bands up to 55 Hz
  808. % BOIsFreqIndices = cell(1, numBands);
  809. % BOIsFreqValues = cell(1, numBands); % To store frequencies for each band
  810. %
  811. % % Calculate the indices and corresponding frequencies for each band
  812. % for b = 1:numBands
  813. % BOIsFreqIndices{b} = find(Fplot >= BOI_filtered(1, b) & Fplot <= BOI_filtered(2, b));
  814. % BOIsFreqValues{b} = Fplot(BOIsFreqIndices{b}); % Map indices to frequencies
  815. % end
  816. %
  817. % % Display the indices and their corresponding frequencies for each band
  818. % for b = 1:numBands
  819. % fprintf('Band: %s (%s)\n', BandName_filtered{b}, BandHzName_filtered{b});
  820. % fprintf('Indices: %s\n', mat2str(BOIsFreqIndices{b}));
  821. % fprintf('Frequencies (Hz): %s\n', mat2str(BOIsFreqValues{b}));
  822. % fprintf('\n');
  823. % end
  824. %%% Loop over channel pairs
  825. for iCh = 1:length(EEGchInd)
  826. for jCh = 1:length(EEGchInd)
  827. if isempty(AvgRealWPLI{iCh, jCh})
  828. continue; % Skip empty cells
  829. end
  830. data = AvgRealWPLI{iCh, jCh}; % data: 3x197 double (stimulation groups x frequencies)
  831. % Preallocate storage for this channel pair
  832. peakWPLI{iCh, jCh} = zeros(size(data, 1), numBands); % stim groups x frequency bands
  833. % Loop over bands of interest
  834. for iBand = 1:numBands
  835. % Find indices corresponding to the frequency band
  836. freqIndices = Fplot >= BOI_filtered(1, iBand) & Fplot <= BOI_filtered(2, iBand);
  837. % Get peak WPLI for each stimulation group
  838. for group = 1:size(data, 1)
  839. peakWPLI{iCh, jCh}(group, iBand) = max(data(group, freqIndices));
  840. % peakWPLI: 32x32 cell of 3x5 double (3 stim groups, 5 BOIs)
  841. end
  842. end
  843. end
  844. end
  845. % Result is stored in peakWPLI{i, j}(group, iBand), where:
  846. % i, j = channel indices
  847. % group = stimulation group
  848. % iBand = frequency band index
  849. % %% TESTING ONLY - Max WPLI Frequency per band
  850. % % Example variables
  851. % maxWPLIFreqHz = cell(length(EEGchInd), numChannels); % To store the frequencies of max WPLI
  852. %
  853. % % Loop through all channel pairs
  854. % for i = 1:numChannels
  855. % for j = 1:numChannels
  856. % if isempty(AvgRealWPLI{i, j})
  857. % continue; % Skip empty cells
  858. % end
  859. %
  860. % data = AvgRealWPLI{i, j}; % 3x197 double (stimulation groups x frequencies)
  861. % maxWPLIFreqHz{i, j} = zeros(size(data, 1), numBands); % Preallocate for groups x bands
  862. %
  863. % % Loop through stimulation groups
  864. % for group = 1:size(data, 1)
  865. % % Loop through each band
  866. % for iBand = 1:numBands
  867. % freqIndices = BOIsFreqIndices{iBand}; % Get indices for the band
  868. % [~, maxIndex] = max(data(group, freqIndices)); % Find index of max value
  869. % maxWPLIFreqHz{i, j}(group, iBand) = Fplot(freqIndices(maxIndex)); % Map to frequency in Hz
  870. % end
  871. % end
  872. % end
  873. % end
  874. %
  875. % % Display example output for a specific channel pair
  876. % exampleChannelPair = [1, 2]; % Change to any valid pair
  877. % if ~isempty(maxWPLIFreqHz{exampleChannelPair(1), exampleChannelPair(2)})
  878. % fprintf('Max WPLI frequencies for channel pair (%d, %d):\n', exampleChannelPair(1), exampleChannelPair(2));
  879. % for group = 1:numGroups
  880. % fprintf(' %s:\n', GroupNameWPLI{group}); % Use the group name
  881. % for iBand = 1:numBands
  882. % fprintf(' Band: %s (%s) - Max WPLI at %.2f Hz\n', ...
  883. % BandName_filtered{iBand}, BandHzName_filtered{iBand}, ...
  884. % maxWPLIFreqHz{exampleChannelPair(1), exampleChannelPair(2)}(group, iBand));
  885. % end
  886. % end
  887. % end
  888. %%% Plot fraction of channel pairs with WPLI greater than p=0.01 permutated WPLI
  889. % Initialize fraction of channel pairs exceeding the threshold
  890. numGroups = 3; % Number of stimulation groups
  891. numBands = size(BOI_filtered, 2); % Number of frequency bands
  892. fractionExceed = zeros(numGroups, numBands);
  893. averagePermutatedWPLI_2to55_top0_0001 = 0.0894;
  894. % averagePermutatedWPLI_2to55_top0_001 = 0.04
  895. % 0.1202; % top quartile threshold
  896. averagePermutatedWPLIvalue = averagePermutatedWPLI_2to55_top0_0001; % top quartile threshold % averagePermutatedWPLI_2to55_top0_001;
  897. % Count total non-empty cells
  898. numChannels = size(AvgRealWPLI, 1);
  899. totalPairs = 0;
  900. for i = 1:numChannels
  901. for j = 1:numChannels
  902. if ~isempty(AvgRealWPLI{i, j})
  903. totalPairs = totalPairs + 1;
  904. end
  905. end
  906. end
  907. % Loop through stimulation groups and frequency bands
  908. for iGroup = 1:numGroups
  909. for iBand = 1:numBands
  910. exceedCount = 0;
  911. % Loop over all channel pairs
  912. for i = 1:numChannels
  913. for j = 1:numChannels
  914. if isempty(AvgRealWPLI{i, j})
  915. continue; % Skip empty cells
  916. end
  917. % Check if the value for the group and band exceeds the threshold
  918. if peakWPLI{i, j}(iGroup, iBand) > averagePermutatedWPLIvalue %
  919. exceedCount = exceedCount + 1;
  920. end
  921. end
  922. end
  923. % Calculate the fraction for this group and band
  924. fractionExceed(iGroup, iBand) = exceedCount / totalPairs;
  925. end
  926. end
  927. %% Peak Alpha WPLI 40Hz vs Light SIGNIFICANT CH ONLY ( stats preceding fig4D stats)
  928. %%% Groups
  929. GroupSubjs_40Hz = 1; % 40Hz group
  930. GroupSubjs_Light = 6; % LightRT group
  931. %%% Get all alpha values
  932. alphaLowerLimitFreqHz = 8;
  933. alphaUpperLimitFreqHz = 13;
  934. % Find indices of values in Fre between 8 and 13 (inclusive)
  935. allAlphaFreIndices = find(Fre >= alphaLowerLimitFreqHz & Fre <= alphaUpperLimitFreqHz);
  936. HitMissTrialType = 3;
  937. WPLIallAlpha = squeeze(WPLIall(:,allAlphaFreIndices,1:32,1:32,HitMissTrialType)); % size(WPLIall) ans = 67 197 40 39 3- 67subs x allFreqs x iCh x jCh x TrialType
  938. %WPLIallAlpha: (67subs, x 11 freqs x 32 x 32chs)
  939. %%% Find the peak WPLI value within alpha (dimension 2) and the corresponding Index (correspond to
  940. % Freq Hz) of the peak alpha WPLI value
  941. [peakAlphaValues4D, peakAlphaIndices4D] = max(WPLIallAlpha, [], 2);
  942. peakAlphaValues = squeeze(peakAlphaValues4D);
  943. peakAlphaIndices = squeeze(peakAlphaIndices4D);
  944. %%%%%%%%%%%%%% pre-process peak alpha WLPI values
  945. % peakAlphaValues = squeeze(peakAlphaValues); % 4d to 3d
  946. % test with single
  947. peakAlphaValueSingleSub = squeeze(peakAlphaValues(1,:,:));
  948. % Extract Peak Alpha WPLI values for each sub into groups: 40 Hz & Light
  949. PeakAlphaWPLIs_40HzGroup = peakAlphaValues(SubjG{GroupSubjs_40Hz},:,:);
  950. PeakAlphaWPLIs_LightGroup = peakAlphaValues(SubjG{GroupSubjs_Light},:,:);
  951. % Get average peak alpha WPLI for each group for all channel pairs
  952. MeanAlphaPeakWPLI_40Hz = squeeze(mean(PeakAlphaWPLIs_40HzGroup, 1));
  953. MeanAlphaPeakWPLI_Light = squeeze(mean(PeakAlphaWPLIs_LightGroup, 1));
  954. %%% Get alpha peak WPLI frequency from peakWPLI_Freq (see: %% Plot bar graph of # channels exceeding
  955. %%% p=0.0001 WPLI threshold (New Fig4C MS) MKA 2025-02-06)
  956. % Preallocate with NaN
  957. peakAlphaFreq_40HzGroup = NaN(numChannels, numChannels);
  958. peakAlphaFreq_LightGroup = NaN(numChannels, numChannels);
  959. % Apply cellfun with error handling
  960. validCells = ~cellfun(@isempty, peakWPLI_Freq); % Logical mask for non-empty cells
  961. peakAlphaFreq_40HzGroup(validCells) = cellfun(@(x) x(1,3), peakWPLI_Freq(validCells));
  962. peakAlphaFreq_LightGroup(validCells) = cellfun(@(x) x(1,3), peakWPLI_Freq(validCells));
  963. % %%%%%%%%%%%% Get frequency value of the peak alpha WPLI
  964. % % Initialize a copy of peakAlphaIndices
  965. % peakAlphaIndicesNaN = peakAlphaIndices;
  966. %
  967. % % test with single participant
  968. % singleSubPeakAlphaIndexAllChPairs = squeeze(peakAlphaIndices(1,:,:)); % should contain values between 1 and 11
  969. %
  970. % % Convert Indices to Correct Corresponding Frequnecy (Hz)
  971. % % Find valid indices (values between 1 and 11)
  972. % validIdx = singleSubPeakAlphaIndexAllChPairs >= 1 & singleSubPeakAlphaIndexAllChPairs <= 11;
  973. %
  974. % % Initialize the output array with NaN, preserving original shape
  975. % singleallchpfreqs = NaN(size(singleSubPeakAlphaIndexAllChPairs));
  976. %
  977. % % Perform mapping only for valid indices
  978. % singleallchpfreqs(validIdx) = alphaFreqHz(singleSubPeakAlphaIndexAllChPairs(validIdx));
  979. %
  980. % % If all alpha WPLI values are NaN, then set the max index to NaN
  981. % peakAlphaIndicesNaN(all(isnan(WPLIallAlpha),2)) = NaN;
  982. %
  983. % % % Reshape the result to 3D: 67x32x32 double
  984. % % PeakAlphaIndicesNaN3D = squeeze(peakAlphaIndicesNaN);
  985. %
  986. % % Convert Indices (3D) to Correct Corresponding Frequency (Hz)
  987. % validIdx = peakAlphaIndicesNaN >= 1 & peakAlphaIndicesNaN <= 11;
  988. %
  989. % % Initialize the output array with NaN, preserving original shape
  990. % PeakAlphaFreqs = NaN(size(peakAlphaIndicesNaN));
  991. %
  992. % % Perform mapping only for valid indices
  993. % PeakAlphaFreqs(validIdx) = alphaFreqHz(peakAlphaIndicesNaN(validIdx));
  994. %
  995. %
  996. % % Extract Peak Alpha WPLI frequency into groups: 40 Hz & Light
  997. % % PeakWPLIallAlpha: 67 subs x 32ch x 32ch
  998. % PeakAlphaFreqs_40HzGroup = PeakAlphaFreqs(SubjG{GroupSubjs_40Hz},:,:);
  999. % PeakAlphaFreqs_LightGroup = PeakAlphaFreqs(SubjG{GroupSubjs_Light},:,:);
  1000. %
  1001. % % Get average peak alpha WPLI frequncy for each group for all channel pairs
  1002. % MeanPeakAlphaFreq_40Hz = squeeze(mean(PeakAlphaFreqs_40HzGroup, 1));
  1003. % MeanPeakAlphaFreq_Light = squeeze(mean(PeakAlphaFreqs_LightGroup, 1));
  1004. %%% Find channel pairs that exceed WPLI threshold
  1005. % WPLI threshold:
  1006. averagePermutatedWPLI_2to55_top0_0001 = 0.0894;
  1007. averagePermutatedWPLIvalue = averagePermutatedWPLI_2to55_top0_0001;
  1008. % display the averagePermutatedWPLIvalue
  1009. fprintf("averagePermutatedWPLI_2to55_top0_0001: %.4f\n", averagePermutatedWPLIvalue);
  1010. % create logical mask of all channel pairs that have a peak band WPLI that exceeds the WPLI threshold (e.g. 0.0894) for 40 Hz and Light
  1011. % Groups. peakWPLI: 32x32 cell of 3x5 doubles
  1012. % SignificantWPLI_40Hz = MeanAlphaPeakWPLI_40Hz > averagePermutatedWPLIvalue;
  1013. % SignificantWPLI_Light = MeanAlphaPeakWPLI_Light > averagePermutatedWPLIvalue;
  1014. SignificantWPLI_40Hz = exceedMask(:,:,1,3);
  1015. SignificantWPLI_Light = exceedMask(:,:,3,3);
  1016. % count the total number of channel pairs that exceed the WPLI threshold, display the result
  1017. numSignificantPairs_40Hz = sum(sum(SignificantWPLI_40Hz));
  1018. numSignificantPairs_Light = sum(sum(SignificantWPLI_Light));
  1019. fprintf('\nNumber of Significant Channel Pairs (exceeding WPLI threshold) in 40Hz Group: %d\n', numSignificantPairs_40Hz);
  1020. fprintf('Number of Significant Channel Pairs (exceeding WPLI threshold) in Light Group: %d\n', numSignificantPairs_Light);
  1021. %%% use logical mask to extract peak alpha frequencies of signficant WPLI (exceeding threshold) channel pairs for 40 Hz and Light from
  1022. % , display the result
  1023. % SignificantAlphaFreqs_40Hz = MeanPeakAlphaFreq_40Hz;
  1024. SignificantAlphaFreqs_40Hz = peakAlphaFreq_40HzGroup .* SignificantWPLI_40Hz;
  1025. SignificantAlphaFreqs_40Hz(~SignificantWPLI_40Hz) = NaN;
  1026. SignificantAlphaFreqs_Light = peakAlphaFreq_LightGroup .* SignificantWPLI_Light;
  1027. SignificantAlphaFreqs_Light(~SignificantWPLI_Light) = NaN;
  1028. % Compute average peak alpha frequency across all channel pairs that exceed WPLI thresold for each
  1029. % group, display the result
  1030. MeanSignificantAlphaFreq_40Hz = mean(SignificantAlphaFreqs_40Hz(SignificantAlphaFreqs_40Hz > 0));
  1031. MeanSignificantAlphaFreq_Light = mean(SignificantAlphaFreqs_Light(SignificantAlphaFreqs_Light > 0));
  1032. fprintf('\nMean Peak Alpha Frequency (Significant 40Hz Group): %.2f Hz\n', MeanSignificantAlphaFreq_40Hz);
  1033. fprintf('Mean Peak Alpha Frequency (Significant Light Group): %.2f Hz\n', MeanSignificantAlphaFreq_Light);
  1034. % Compute median values
  1035. MedianSignificantAlphaFreq_40Hz = median(SignificantAlphaFreqs_40Hz(SignificantAlphaFreqs_40Hz > 0));
  1036. MedianSignificantAlphaFreq_Light = median(SignificantAlphaFreqs_Light(SignificantAlphaFreqs_Light > 0));
  1037. fprintf('\nMedian Peak Alpha Frequency (Significant 40Hz Group): %.2f Hz\n', MedianSignificantAlphaFreq_40Hz);
  1038. fprintf('Median Peak Alpha Frequency (Significant Light Group): %.2f Hz\n', MedianSignificantAlphaFreq_Light);
  1039. % Compute 25th and 75th percentiles
  1040. Q1_40Hz = prctile(SignificantAlphaFreqs_40Hz(SignificantAlphaFreqs_40Hz > 0), 25);
  1041. Q3_40Hz = prctile(SignificantAlphaFreqs_40Hz(SignificantAlphaFreqs_40Hz > 0), 75);
  1042. Q1_Light = prctile(SignificantAlphaFreqs_Light(SignificantAlphaFreqs_Light > 0), 25);
  1043. Q3_Light = prctile(SignificantAlphaFreqs_Light(SignificantAlphaFreqs_Light > 0), 75);
  1044. fprintf('\n25th Percentile (Significant 40Hz Group): %.2f Hz\n', Q1_40Hz);
  1045. fprintf('75th Percentile (Significant 40Hz Group): %.2f Hz\n', Q3_40Hz);
  1046. fprintf('25th Percentile (Significant Light Group): %.2f Hz\n', Q1_Light);
  1047. fprintf('75th Percentile (Significant Light Group): %.2f Hz\n', Q3_Light);
  1048. % Compute STE (Standard Error of the Mean)
  1049. std_40Hz = nanstd(SignificantAlphaFreqs_40Hz(:)); % Standard deviation
  1050. std_Light = nanstd(SignificantAlphaFreqs_Light(:));
  1051. N_40Hz = sum(~isnan(SignificantAlphaFreqs_40Hz(:))); % Sample size
  1052. N_Light = sum(~isnan(SignificantAlphaFreqs_Light(:)));
  1053. STE_40Hz = std_40Hz / sqrt(N_40Hz);
  1054. STE_Light = std_Light / sqrt(N_Light);
  1055. fprintf('\nStandard Error of the Mean (STE) - 40Hz Group: %.4f Hz\n', STE_40Hz);
  1056. fprintf('Standard Error of the Mean (STE) - Light Group: %.4f Hz\n', STE_Light);
  1057. % Test the distribution of alpha peak frequencies for normality, display results
  1058. [h_40Hz, p_40Hz] = kstest(SignificantAlphaFreqs_40Hz(:));
  1059. [h_Light, p_Light] = kstest(SignificantAlphaFreqs_Light(:));
  1060. fprintf('\nKolmogorov–Smirnov Normality Test (40Hz Group): p = %.8f\n', p_40Hz);
  1061. fprintf('Kolmogorov–Smirnov Normality Test (Light Group): p = %.8f\n', p_Light);
  1062. % Create violin plot of 40 Hz and Light alpha peak frequency distributions
  1063. figure;
  1064. v = violinplot([SignificantAlphaFreqs_40Hz(:), SignificantAlphaFreqs_Light(:)], {'40Hz', 'Light'});
  1065. ylabel('Peak Alpha Frequency (Hz)');
  1066. title(sprintf('Violin Plot of Peak Alpha Frequency Distributions\nWPLI threshold: %.4f', averagePermutatedWPLIvalue));
  1067. % Compute means and number of non-NaN data points
  1068. mean_40Hz = nanmean(SignificantAlphaFreqs_40Hz(:));
  1069. mean_Light = nanmean(SignificantAlphaFreqs_Light(:));
  1070. N_40Hz = sum(~isnan(SignificantAlphaFreqs_40Hz(:)));
  1071. N_Light = sum(~isnan(SignificantAlphaFreqs_Light(:)));
  1072. % Annotate mean values on the plot
  1073. hold on;
  1074. plot(1, mean_40Hz, 'kd', 'MarkerFaceColor', 'k', 'MarkerSize', 8); % Mean for 40Hz
  1075. plot(2, mean_Light, 'kd', 'MarkerFaceColor', 'k', 'MarkerSize', 8); % Mean for Light
  1076. % Annotate median values on the plot
  1077. plot(1, MedianSignificantAlphaFreq_40Hz, 'bs', 'MarkerFaceColor', 'b', 'MarkerSize', 8); % Median for 40Hz
  1078. plot(2, MedianSignificantAlphaFreq_Light, 'bs', 'MarkerFaceColor', 'b', 'MarkerSize', 8); % Median for Light
  1079. % Annotate 25th and 75th percentiles on the plot
  1080. plot(1, Q1_40Hz, 'm^', 'MarkerFaceColor', 'm', 'MarkerSize', 6); % 25th Percentile 40Hz
  1081. plot(1, Q3_40Hz, 'm^', 'MarkerFaceColor', 'm', 'MarkerSize', 6); % 75th Percentile 40Hz
  1082. plot(2, Q1_Light, 'm^', 'MarkerFaceColor', 'm', 'MarkerSize', 6); % 25th Percentile Light
  1083. plot(2, Q3_Light, 'm^', 'MarkerFaceColor', 'm', 'MarkerSize', 6); % 75th Percentile Light
  1084. % Display text for mean, median, and percentiles
  1085. text(1, mean_40Hz + 0.2, sprintf('Mean: %.2f Hz', mean_40Hz), 'HorizontalAlignment', 'center');
  1086. text(2, mean_Light + 0.2, sprintf('Mean: %.2f Hz', mean_Light), 'HorizontalAlignment', 'center');
  1087. text(1, MedianSignificantAlphaFreq_40Hz - 0.2, sprintf('Median: %.2f Hz', MedianSignificantAlphaFreq_40Hz), 'HorizontalAlignment', 'center', 'Color', 'b');
  1088. text(2, MedianSignificantAlphaFreq_Light - 0.2, sprintf('Median: %.2f Hz', MedianSignificantAlphaFreq_Light), 'HorizontalAlignment', 'center', 'Color', 'b');
  1089. text(1, Q1_40Hz - 0.3, sprintf('Q1: %.2f Hz', Q1_40Hz), 'HorizontalAlignment', 'center', 'Color', 'm');
  1090. text(1, Q3_40Hz + 0.3, sprintf('Q3: %.2f Hz', Q3_40Hz), 'HorizontalAlignment', 'center', 'Color', 'm');
  1091. text(2, Q1_Light - 0.3, sprintf('Q1: %.2f Hz', Q1_Light), 'HorizontalAlignment', 'center', 'Color', 'm');
  1092. text(2, Q3_Light + 0.3, sprintf('Q3: %.2f Hz', Q3_Light), 'HorizontalAlignment', 'center', 'Color', 'm');
  1093. % Display sample sizes
  1094. text(1, min(ylim) + 1, sprintf('N = %d', N_40Hz), 'HorizontalAlignment', 'center');
  1095. text(2, min(ylim) + 1, sprintf('N = %d', N_Light), 'HorizontalAlignment', 'center');
  1096. hold off;
  1097. % Rank sum test: testing if Light has significantly greater alpha peak frequency than 40 Hz
  1098. [p_ranksum, h_ranksum] = ranksum(SignificantAlphaFreqs_40Hz(:), SignificantAlphaFreqs_Light(:), 'tail', 'left');
  1099. fprintf('\nRank Sum Test p-value: %.16f\n', p_ranksum);
  1100. % Display peak alpha frequency of channel pairs exceeding WPLI threshold in tables (32x32)
  1101. ChannelLabels = {'Fp1', 'AF3', 'F7', 'F3', 'FC1', 'FC5', 'T7', 'C3', 'CP1', 'CP5', ...
  1102. 'P7', 'P3', 'Pz', 'PO3', 'O1', 'Oz', 'O2', 'PO4', 'P4', 'P8', ...
  1103. 'CP6', 'CP2', 'C4', 'T8', 'FC6', 'FC2', 'F4', 'F8', 'AF4', 'Fp2', 'Fz', 'Cz'};
  1104. fprintf('\nMean Peak Alpha Frequency (Significant channel pairs only - in 40Hz Group):\n');
  1105. disp(array2table(SignificantAlphaFreqs_40Hz, 'VariableNames', ChannelLabels, 'RowNames', ChannelLabels));
  1106. fprintf('\nMean Peak Alpha Frequency (Significant channel pairs only - in Light Group):\n');
  1107. disp(array2table(SignificantAlphaFreqs_Light, 'VariableNames', ChannelLabels, 'RowNames', ChannelLabels));
  1108. fprintf('\nMean Peak Alpha WPLI (All channel pairs in 40Hz Group):\n');
  1109. disp(array2table(MeanAlphaPeakWPLI_40Hz, 'VariableNames', ChannelLabels, 'RowNames', ChannelLabels));
  1110. fprintf('\nMean Peak Alpha WPLI (All channel pairs in Light Group):\n');
  1111. disp(array2table(MeanAlphaPeakWPLI_Light, 'VariableNames', ChannelLabels, 'RowNames', ChannelLabels));
  1112. %% Peak Alpha WPLI 40Hz vs Light SIGNIFICANT CH ONLY (abandoned approached)
  1113. % Initialize logical masks for valid channel pairs
  1114. numChannels = 32;
  1115. validChannels_40Hz = false(numChannels, numChannels);
  1116. validChannels_Light = false(numChannels, numChannels);
  1117. group40HzIndex = 1;
  1118. groupLightIndex = 3;
  1119. % Identify valid channels where WPLI exceeds threshold in each group
  1120. for i = 1:numChannels
  1121. for j = 1:numChannels
  1122. if isempty(peakWPLI{i, j})
  1123. continue; % Skip empty cells
  1124. end
  1125. % Check if peak WPLI is greater than threshold for 40Hz and Light groups
  1126. if peakWPLI{i, j}(group40HzIndex, 1) > averagePermutatedWPLIvalue
  1127. validChannels_40Hz(i, j) = true;
  1128. end
  1129. if peakWPLI{i, j}(groupLightIndex, 1) > averagePermutatedWPLIvalue
  1130. validChannels_Light(i, j) = true;
  1131. end
  1132. end
  1133. end
  1134. % Find common channel pairs where both groups exceed the threshold
  1135. validChannels = validChannels_40Hz & validChannels_Light;
  1136. % Collect flattened peak alpha frequencies for valid channels only
  1137. selectedFreqs_40Hz = [];
  1138. selectedFreqs_Light = [];
  1139. for i = 1:numChannels
  1140. for j = 1:numChannels
  1141. if validChannels(i, j)
  1142. selectedFreqs_40Hz = [selectedFreqs_40Hz; PeakAlphaFreqs_40HzGroup(:, i, j)];
  1143. selectedFreqs_Light = [selectedFreqs_Light; PeakAlphaFreqs_LightGroup(:, i, j)];
  1144. end
  1145. end
  1146. end
  1147. %% Plot the results (no counts)
  1148. figure;
  1149. for iGroup = 1:numGroups
  1150. subplot(1, numGroups, iGroup);
  1151. bar(fractionExceed(iGroup, :));
  1152. title(GroupNameWPLI{iGroup}); % Use group name as title
  1153. xlabel('Frequency Band');
  1154. ylabel('Fraction of Channel Pairs');
  1155. % Customize x-axis labels with both BandName and BandHzName
  1156. xticks(1:numBands);
  1157. xticklabels(arrayfun(@(iBand) [BandName_filtered{iBand}, ' (', BandHzName_filtered{iBand}, ')'], ...
  1158. 1:numBands, 'UniformOutput', false));
  1159. xtickangle(45);
  1160. ylim([0 1]); % Fractions are between 0 and 1
  1161. end
  1162. % Add super title
  1163. sgtitle('Fraction of Channel Pairs Exceeding average permutated p=0.01 WPLI value');
  1164. %% Same as above but with counts above each bar
  1165. figure;
  1166. for iGroup = 1:numGroups
  1167. subplot(1, numGroups, iGroup);
  1168. barHandle = bar(fractionExceed(iGroup, :));
  1169. title(GroupNameWPLI{iGroup}); % Use group name as title
  1170. xlabel('Frequency Band');
  1171. ylabel('Fraction of Channel Pairs');
  1172. % Customize x-axis labels with both BandName and BandHzName
  1173. xticks(1:numBands);
  1174. xticklabels(arrayfun(@(iBand) [BandName_filtered{iBand}, ' (', BandHzName_filtered{iBand}, ')'], ...
  1175. 1:numBands, 'UniformOutput', false));
  1176. xtickangle(45);
  1177. ylim([0 1]); % Fractions are between 0 and 1
  1178. %%% Add exceed count above each bar
  1179. barHeights = fractionExceed(group, :);
  1180. for iBand = 1:numBands
  1181. % Position the text slightly above the bar height
  1182. text(iBand, barHeights(iBand) + 0.02, num2str(round(barHeights(iBand) * totalPairs)), ...
  1183. 'HorizontalAlignment', 'center', 'FontSize', 10);
  1184. end
  1185. end
  1186. % Add super title
  1187. sgtitle('Fraction of Channel Pairs Exceeding average permutated p=0.0001 WPLI value');
  1188. %% Bar graph with counts below
  1189. figure;
  1190. for group = 1:numGroups
  1191. subplot(1, numGroups, group);
  1192. barHandle = bar(fractionExceed(group, :));
  1193. title(GroupNameWPLI{group}); % Use group name as title
  1194. xlabel('Frequency Band');
  1195. ylabel('Fraction of Channel Pairs');
  1196. % Customize x-axis labels with both BandName and BandHzName
  1197. xticks(1:numBands);
  1198. xticklabels(arrayfun(@(iBand) [BandName_filtered{iBand}, ' (', BandHzName_filtered{iBand}, ')'], ...
  1199. 1:numBands, 'UniformOutput', false));
  1200. xtickangle(45);
  1201. ylim([0 1]); % Fractions are between 0 and 1
  1202. % Add exceed count below the top of each bar
  1203. barHeights = fractionExceed(group, :);
  1204. for iBand = 1:numBands
  1205. % Position the text slightly below the bar height
  1206. text(iBand, barHeights(iBand) - 0.02, num2str(round(barHeights(iBand) * totalPairs)), ...
  1207. 'HorizontalAlignment', 'center', 'FontSize', 10, 'VerticalAlignment', 'top');
  1208. end
  1209. end
  1210. % Add super title
  1211. sgtitle('Fraction of Channel Pairs Exceeding average permutated p=0.01 WPLI value');
  1212. %% Plot the results 2025-01-08 MKA fig4D supp? - top quartile
  1213. % BOI=[1 4 8 13 30 39.5;4 8 13 30 37 41.5];
  1214. % BOI2HzFiveBands=[2 4 8 13 39.5;4 8 13 30 41.5];
  1215. % BandNameFiveBands={'Delta','Theta','Alpha','Beta','Gamma-E'};
  1216. % BandHzName={'1-4 Hz','4-8 Hz','8-13Hz','13-30 Hz','39-41 Hz'};
  1217. % BandHzName2HzFiveBands={'2-4 Hz','4-8 Hz','8-13Hz','13-30 Hz','39-41 Hz'};
  1218. BOI2Hz=[2 4 8 8 10 13 30 39.5 43 55;4 8 13 10 13 30 37 41.5 55 100];
  1219. BandName={'Delta','Theta','Alpha','LowAlpha','HighAlpha','Beta','Gamma-1','Gamma-E', 'Gamma-55','Gamma-2'};
  1220. BandHzName2Hz={'2-4 Hz','4-8 Hz','8-13Hz','8-10Hz','10-13Hz','13-30 Hz','30-37 Hz','39-41 Hz','43-55 Hz', '55-100 Hz'};
  1221. numBands = size(BOI2Hz, 2);
  1222. %%% Get WPLI for each stim group for each BOI and ch Pair
  1223. %%%% Loop thru all ch pairs
  1224. for i = 1:numChannels
  1225. for j = 1:numChannels
  1226. if isempty(AvgRealWPLI{i, j})
  1227. continue; % Skip empty cells
  1228. end
  1229. data = AvgRealWPLI{i, j}; % 3x197 double (stimulation groups x frequencies)
  1230. % Preallocate storage for this channel pair
  1231. peakWPLI{i, j} = zeros(size(data, 1), numBands); % stim groups x frequency bands
  1232. % Loop over bands of interest
  1233. for iBand = 1:numBands
  1234. % Find indices corresponding to the frequency band
  1235. freqIndices = Fplot >= BOI2Hz(1, iBand) & Fplot <= BOI2Hz(2, iBand);
  1236. % Get peak WPLI for each stimulation group
  1237. for group = 1:size(data, 1)
  1238. peakWPLI{i, j}(group, iBand) = max(data(group, freqIndices));
  1239. end
  1240. end
  1241. end
  1242. end
  1243. % Result is stored in peakWPLI{i, j}(group, iBand), where:
  1244. % i, j = channel indices
  1245. % group = stimulation group
  1246. % iBand = frequency band index
  1247. %%% Plot fraction of channel pairs with WPLI greater than p=0.01 permutated WPLI
  1248. % Example variables (replace these with actual data)
  1249. % averagePermutatedWPLItop0_01 = 0.5; % Replace with actual value
  1250. % GroupNameWPLI = {'Group 1', 'Group 2', 'Group 3'}; % Replace with actual group names
  1251. % Initialize fraction of channel pairs exceeding the threshold
  1252. numGroups = 3; % Number of stimulation groups
  1253. numBands = size(BOI2Hz, 2)-1; % Number of frequency bands
  1254. fractionExceed = zeros(numGroups, numBands);
  1255. % averagePermutatedWPLI_2to55_top0_0001 = 0.0894
  1256. % averagePermutatedWPLI_2to55_top0_001 = 0.04
  1257. averagePermutatedWPLIvalue = 0.1202; % top quartile threshold %
  1258. % Count total non-empty cells
  1259. numChannels = size(AvgRealWPLI, 1);
  1260. totalPairs = 0;
  1261. for i = 1:numChannels
  1262. for j = 1:numChannels
  1263. if ~isempty(AvgRealWPLI{i, j})
  1264. totalPairs = totalPairs + 1;
  1265. end
  1266. end
  1267. end
  1268. % Loop through stimulation groups and frequency bands
  1269. for group = 1:numGroups
  1270. for iBand = 1:numBands
  1271. exceedCount = 0;
  1272. % Loop over all channel pairs
  1273. for i = 1:numChannels
  1274. for j = 1:numChannels
  1275. if isempty(AvgRealWPLI{i, j})
  1276. continue; % Skip empty cells
  1277. end
  1278. % Check if the value for the group and band exceeds the threshold
  1279. if peakWPLI{i, j}(group, iBand) > averagePermutatedWPLIvalue %averagePermutatedWPLItop0_01
  1280. exceedCount = exceedCount + 1;
  1281. end
  1282. end
  1283. end
  1284. % Calculate the fraction for this group and band
  1285. fractionExceed(group, iBand) = exceedCount / totalPairs;
  1286. end
  1287. end
  1288. figure;
  1289. for iGroup = 1:numGroups
  1290. subplot(1, numGroups, iGroup);
  1291. barHandle = bar(fractionExceed(iGroup, :));
  1292. title(GroupNameWPLI{iGroup}); % Use iGroup name as title
  1293. xlabel('Frequency Band');
  1294. ylabel('Fraction of Channel Pairs');
  1295. % Customize x-axis labels with both BandName and BandHzName
  1296. xticks(1:numBands);
  1297. xticklabels(arrayfun(@(iBand) [BandName{iBand}, ' (', BandHzName2Hz{iBand}, ')'], ...
  1298. 1:numBands, 'UniformOutput', false));
  1299. xtickangle(45);
  1300. ylim([0 1]); % Fractions are between 0 and 1
  1301. % Add exceed count conditionally above or below the bar
  1302. barHeights = fractionExceed(iGroup, :);
  1303. for iBand = 1:numBands
  1304. exceedCount = round(barHeights(iBand) * totalPairs); % Calculate exceed count
  1305. if barHeights(iBand) < 0.2
  1306. % Place count above the bar
  1307. text(iBand, barHeights(iBand) + 0.04, num2str(exceedCount), ...
  1308. 'HorizontalAlignment', 'center', 'FontSize', 10);
  1309. else
  1310. % Place count below the top of the bar
  1311. text(iBand, barHeights(iBand) - 0.01, num2str(exceedCount), ...
  1312. 'HorizontalAlignment', 'center', 'FontSize', 10, 'VerticalAlignment', 'top');
  1313. end
  1314. end
  1315. % Store percentage of Alpha band exceedance for this iGroup
  1316. alphaExceedPercent(iGroup) = barHeights(3) * 100; % Alpha band is iBand = 4
  1317. end
  1318. % Add super title
  1319. sgtitle(['Fraction of Channel Pairs Exceeding top quartile WPLI value: ' num2str(averagePermutatedWPLIvalue)]);
  1320. % Calculate and display the average percentage of Alpha band exceedance
  1321. averageAlphaPercent = mean(alphaExceedPercent);
  1322. disp(['Average Percentage of Channel Pairs Exceeding Threshold for Alpha Band: ', num2str(averageAlphaPercent), '%']);
  1323. % Display percentages for each group
  1324. for group = 1:numGroups
  1325. disp(['Group ', GroupNameWPLI{group}, ': ', num2str(alphaExceedPercent(group)), '%']);
  1326. end
  1327. %% Plot bar graph of # channels exceeding p=0.0001 WPLI threshold (New Fig4C MS) MKA 2025-02-06
  1328. BOI2HzFiveBands=[2 4 8 13 30; 4 8 13 30 37];
  1329. BandNameFiveBands={'Delta','Theta','Alpha','Beta','Slow Gamma'};
  1330. BandHzName2HzFiveBands={'2-4 Hz','4-8 Hz','8-13Hz','13-30 Hz','30-37 Hz'};
  1331. numBands = size(BOI2HzFiveBands, 2);
  1332. %%% Initialize storage for peak frequencies
  1333. peakWPLI_Freq = cell(numChannels, numChannels);
  1334. %%% Loop over channel pairs
  1335. for i = 1:numChannels
  1336. for j = 1:numChannels
  1337. if isempty(AvgRealWPLI{i, j})
  1338. continue; % Skip empty cells
  1339. end
  1340. AvgRealWPLI_ijChPair = AvgRealWPLI{i, j}; % AvgRealWPLI_ijChPair = 3x197 double (stimulation groups x frequencies)
  1341. % Preallocate storage for this channel pair
  1342. peakWPLI{i, j} = zeros(size(AvgRealWPLI_ijChPair, 1), numBands); % stim groups x frequency bands
  1343. peakWPLI_Freq{i, j} = zeros(size(AvgRealWPLI_ijChPair, 1), numBands); % Store peak frequency
  1344. % Loop over bands of interest
  1345. for iBand = 1:numBands
  1346. % Find indices corresponding to the frequency band
  1347. iBandFreqIndices = Fplot >= BOI2HzFiveBands(1, iBand) & Fplot <= BOI2HzFiveBands(2, iBand);
  1348. Fplot_iBand = Fplot(iBandFreqIndices); % Extract frequencies in this band
  1349. % Get peak WPLI for each stimulation group
  1350. for iGroup = 1:size(AvgRealWPLI_ijChPair, 1)
  1351. AvgRealWPLI_ijChPair_iBand_iGroup = AvgRealWPLI_ijChPair(iGroup, iBandFreqIndices); %
  1352. % Find peak WPLI value
  1353. [peakWPLI{i, j}(iGroup, iBand), maxIdx] = max(AvgRealWPLI_ijChPair_iBand_iGroup);
  1354. % Store corresponding frequency
  1355. peakWPLI_Freq{i, j}(iGroup, iBand) = Fplot_iBand(maxIdx);
  1356. end
  1357. end
  1358. end
  1359. end
  1360. % Result is stored in peakWPLI{i, j}(group, iBand), where:
  1361. % i, j = channel indices
  1362. % group = stimulation group
  1363. % iBand = frequency band index
  1364. %%% Plot fraction of channel pairs with WPLI greater than p=0.01 permutated WPLI
  1365. % Example variables (replace these with actual data)
  1366. % averagePermutatedWPLItop0_01 = 0.5; % Replace with actual value
  1367. % GroupNameWPLI = {'Group 1', 'Group 2', 'Group 3'}; % Replace with actual group names
  1368. % Initialize fraction of channel pairs exceeding the threshold
  1369. numGroups = 3; % Number of stimulation groups
  1370. numBands = size(BOI2HzFiveBands, 2); % Number of frequency bands
  1371. fractionExceed = zeros(numGroups, numBands);
  1372. % averagePermutatedWPLI_2to55_top0_0001 = 0.0894
  1373. % averagePermutatedWPLI_2to55_top0_001 = 0.04
  1374. averagePermutatedWPLIvalue = averagePermutatedWPLI_2to55_top0_0001; % top quartile threshold % averagePermutatedWPLI_2to55_top0_001;
  1375. % Count total non-empty cells
  1376. numChannels = size(AvgRealWPLI, 1);
  1377. totalPairs = 0;
  1378. %%% count the number of total channel pairs based on number of AvgRealWPLI values (use mask)
  1379. for i = 1:numChannels
  1380. for j = 1:numChannels
  1381. if ~isempty(AvgRealWPLI{i, j})
  1382. totalPairs = totalPairs + 1;
  1383. end
  1384. end
  1385. end
  1386. totalPairsMaskSum = sum(~cellfun(@isempty, AvgRealWPLI), 'all');
  1387. % Preallocate arrays
  1388. peakWPLIarray = nan(numChannels, numChannels, numGroups, numBands);
  1389. exceedMask = false(numChannels, numChannels, numGroups, numBands);
  1390. fractionExceed = zeros(numGroups, numBands);
  1391. % Create a logical mask for non-empty channel pairs
  1392. validPairs = ~cellfun(@isempty, AvgRealWPLI);
  1393. % Loop through stimulation groups and frequency bands
  1394. for iGroup = 1:numGroups
  1395. for iBand = 1:numBands
  1396. % % Extract peakWPLI values into an array (set NaN for empty pairs)
  1397. % peakWPLIarray(validPairs, iGroup, iBand) = cellfun(@(x) x(iGroup, iBand), peakWPLI(validPairs));
  1398. % peakWPLIarray(~validPairs) = NaN; % Ensure empty cells remain NaN
  1399. % Extract peakWPLI values into an array (set NaN for empty pairs)
  1400. tempValues = nan(numChannels, numChannels); % Temporary storage
  1401. tempValues(validPairs) = cellfun(@(x) x(iGroup, iBand), peakWPLI(validPairs));
  1402. % Store in preallocated array
  1403. peakWPLIarray(:, :, iGroup, iBand) = tempValues;
  1404. % Create a logical mask for values exceeding the threshold
  1405. exceedMask(:, :, iGroup, iBand) = peakWPLIarray(:, :, iGroup, iBand) > averagePermutatedWPLIvalue;
  1406. % Count the number of exceeding pairs
  1407. exceedCount = sum(exceedMask(:, :, iGroup, iBand), 'all');
  1408. % Calculate the fraction
  1409. fractionExceed(iGroup, iBand) = exceedCount / totalPairs;
  1410. end
  1411. end
  1412. % Loop through stimulation groups and frequency bands
  1413. for group = 1:numGroups
  1414. for iBand = 1:numBands
  1415. exceedCount = 0;
  1416. % Loop over all channel pairs
  1417. for i = 1:numChannels
  1418. for j = 1:numChannels
  1419. if isempty(AvgRealWPLI{i, j})
  1420. continue; % Skip empty cells
  1421. end
  1422. % Check if the value for the group and band exceeds the threshold
  1423. if peakWPLI{i, j}(group, iBand) > averagePermutatedWPLIvalue %averagePermutatedWPLItop0_01
  1424. exceedCount = exceedCount + 1;
  1425. end
  1426. end
  1427. end
  1428. % Calculate the fraction for this group and band
  1429. fractionExceed(group, iBand) = exceedCount / totalPairs;
  1430. end
  1431. end
  1432. %%% Plot High Functional Connectivity Bar Plot *** fig4c
  1433. figure;
  1434. for iGroup = 1:numGroups
  1435. subplot(1, numGroups, iGroup);
  1436. barHandle = bar(fractionExceed(iGroup, :));
  1437. if iGroup ~= 3
  1438. title(GroupNameWPLI{iGroup}); % Use iGroup name as title
  1439. else
  1440. title("Light"); % Use "Light" instead of "LightRT"
  1441. end
  1442. if iGroup ==2
  1443. xlabel('Frequency Band');
  1444. end
  1445. if iGroup == 1
  1446. ylabel('Fraction of Channel Pairs');
  1447. end
  1448. % Customize x-axis labels with both BandName and BandHzName
  1449. iBand = 1;
  1450. xticks(1:numBands);
  1451. xticklabels(arrayfun(@(iBand) [BandNameFiveBands{iBand}], ...
  1452. 1:numBands, 'UniformOutput', false));
  1453. xtickangle(45);
  1454. ylim([0 0.55]); % Fractions are between 0 and 1
  1455. % Add exceed count conditionally above or below the bar
  1456. barHeights = fractionExceed(iGroup, :);
  1457. for iBand = 1:numBands
  1458. exceedCount = round(barHeights(iBand) * totalPairs); % Calculate exceed count
  1459. if barHeights(iBand) < 0.2
  1460. % Place count above the bar
  1461. text(iBand, barHeights(iBand) + 0.06500000000000000000000001, num2str(exceedCount), ...
  1462. 'HorizontalAlignment', 'center', 'FontSize', 10, 'VerticalAlignment', 'top'); %
  1463. else
  1464. % Place count below the top of the bar
  1465. text(iBand, barHeights(iBand) + 0.065, num2str(exceedCount), ...
  1466. 'HorizontalAlignment', 'center', 'FontSize', 10, 'VerticalAlignment', 'top');
  1467. end
  1468. end
  1469. % Store percentage of Alpha band exceedance for this iGroup
  1470. alphaExceedPercent(iGroup) = barHeights(3) * 100; % Alpha band is iBand = 4
  1471. end
  1472. % Add super title
  1473. sgtitle(['Fraction of Channel Pairs Out of 492 Exceeding permutated p=0.0001 WPLI value: ' num2str(averagePermutatedWPLIvalue)], 'FontSize', 6);
  1474. % Calculate and display the average percentage of Alpha band exceedance
  1475. averageAlphaPercent = mean(alphaExceedPercent);
  1476. disp(['Average Percentage of Channel Pairs Exceeding Threshold for Alpha Band: ', num2str(averageAlphaPercent), '%']);
  1477. % Display percentages for each group
  1478. for group = 1:numGroups
  1479. disp(['Group ', GroupNameWPLI{group}, ': ', num2str(alphaExceedPercent(group)), '%']);
  1480. end
  1481. % Set figure size
  1482. set(gcf, 'PaperUnits', 'inches', 'PaperPosition', [0 0 4 2.5]); % 6x4 inches
  1483. % Save as SVG
  1484. saveDateFig4C = datestr(datetime, 'yy-mm-dd_HHMMSSFFF');
  1485. saveas(gcf, ['Fig4C_WPLI_exceed' saveDateFig4C '.svg']);
  1486. saveas(gcf, ['Fig4C_WPLI_exceed' saveDateFig4C '.png']);
  1487. %% Get Alpha Peak Amplitude and Alpha Peak Frequency *** fig4b
  1488. iWPLICh = 1; % 1 = FP1
  1489. jWPLICh = 5; % 5 = FC1
  1490. for iStimGroup = 1:length(SubjGWPLI)
  1491. groupMean = mean(DataPlot{iStimGroup});
  1492. alphaMean = groupMean(13:23);
  1493. [alphaPeakAmplitude, alphaPeakIndex] = max(groupMean(13:23));
  1494. alphaFreqs = Fplot(13:23); % 13 to 23 should be 8 to 13 Hz
  1495. alphaPeakFrequency = alphaFreqs(alphaPeakIndex);
  1496. alphaPeakAmplitudeList(iWPLICh,jWPLICh,iStimGroup) = alphaPeakAmplitude;
  1497. alphaPeakFrequencyList(iWPLICh,jWPLICh,iStimGroup) = alphaPeakFrequency;
  1498. end
  1499. % Create the actual plot fig4b (commented for getting alpha measurements) UNCOMMENT TO PLOT!
  1500. figure;
  1501. xrange = [0 37]; % x axis range (frequency)
  1502. subplot('Position',[0.1 0.1 0.88 0.88])
  1503. % Cut data down to 55Hz
  1504. % Fplot55=Fplot(1:107);
  1505. % DataPlot55 = cellfun(@(x) x(1:107), DataPlot, 'UniformOutput', false); % Works like this-> DataPlot55=DataPlot(1:107);
  1506. RateHist_GroupPlot(Fplot,DataPlot,FlickerColor,ParamWPLI); %% Lu's function
  1507. % Add a dotted line at averagePermutatedWPLIvalue
  1508. yline(averagePermutatedWPLIvalue, '--', 'Color', 'k', 'LineWidth', 1.2);
  1509. text(range(xrange)/2-5,ParamWPLI.Ytick(end),[EEGch{iWPLICh} '-' EEGch{jWPLICh}],'FontSize', 10);
  1510. % set(gca,'xlim',[0 100],'xtick',[0:20:120],'ylim',[0 0.3],'ytick',ParamWPLI.Ytick);
  1511. set(gca,'xlim',xrange,'xtick',[1 4 8 13 30 37 40 50 60 80 100],'ylim',[0 0.15],'ytick',[0 0.05 0.1 0.15]);
  1512. xlabel('Frequency Hz')
  1513. % xlim(xrange);
  1514. ylabel('WPLI')
  1515. title(['EEG Channel Pair: ' EEGch{iWPLICh} '-' EEGch{jWPLICh}]);
  1516. ax=gca;
  1517. ax.XGrid = 'on';
  1518. ax.YGrid = 'off';
  1519. ax.FontSize = 7; % Set the desired font size for tick marks
  1520. LuFontStandard
  1521. papersizePX=[0 0 8 8];
  1522. papersizePX=1.3*[0 0 5 3.2]; % 09/09/24
  1523. set(gcf, 'PaperUnits', 'centimeters');
  1524. set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  1525. saveas(gcf,[SubSaveFig '40LightRandBands_' EEGch{iWPLICh} '-' EEGch{jWPLICh}],'svg');
  1526. % saveas(gcf,[SubSaveFig '40LightRandBands_' EEGch{iWPLICh} '-' EEGch{jWPLICh}],'tiff');
  1527. close all
  1528. nonzerosGroup1= nonzeros(alphaPeakAmplitudeList(:,:,1));
  1529. nonzerosGroup2= nonzeros(alphaPeakAmplitudeList(:,:,2));
  1530. nonzerosGroup3= nonzeros(alphaPeakAmplitudeList(:,:,3));
  1531. AlphaPeakSaveFilename = fullfile(SaveFolder, 'AlphaPeakAmpAndFreq') %#ok<NOPTS>
  1532. save(AlphaPeakSaveFilename,"alphaPeakAmplitudeList","alphaPeakFrequencyList");
  1533. % papersizePX=[0 0 6*length(EEGchInd) 6*length(EEGchInd)];
  1534. % set(gcf, 'PaperUnits', 'centimeters');
  1535. % set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  1536. % saveas(gcf,[SaveTemp num2str(FBand(1)) '-' num2str(FBand(2)) 'HzAllCh40Random' TrialTypeName{iCom}],'pdf');
  1537. % saveas(gcf,[SaveTemp num2str(FBand(1)) '-' num2str(FBand(2)) 'HzAllCh40Random' TrialTypeName{iCom}],'png');
  1538. % saveas(gcf,[SaveTemp num2str(FBand(1)) '-' num2str(FBand(2)) 'HzAllCh40Random' TrialTypeName{iCom} '.eps'],'epsc');
  1539. close all
  1540. close all
  1541. %% Frequency Band Definition, for maps - May need to start running from here for Fig4
  1542. BOI=[1 4 8 13 30 39.5 43;4 8 13 30 37 41.5 100];
  1543. BandName={'Delta','Theta','Alpha','Beta','Gamma-1','Gamma-E','Gamma-2'};
  1544. % Gamma-E is 40
  1545. BandHzName={'1-4 Hz','4-8 Hz','8-13Hz','13-30 Hz','30-37 Hz','39-41 Hz','43-100 Hz'};
  1546. FreqFunc{1}=@nanmean;
  1547. FreqFunc{2}=@nanmedian;
  1548. FreqFunc{3}=@nanmax;
  1549. FreqFuncNames={'mean','median','peak'}; % Check three different spots in band - previously called 'FunGroupName'
  1550. %comparison between groups w/ freq data
  1551. diffTPmap=zeros(length(ChanEEGLab),length(ChanEEGLab),size(BOI,2),length(TrialType)); % for T test
  1552. diffRPmap=diffTPmap;
  1553. diffTmap=diffTPmap;
  1554. rSpear=diffTPmap;
  1555. pSpear=diffTPmap;
  1556. % % TNodeTh=10;
  1557. % FC_BrainEEGLu(ChanPos,AdjWeight,NodeWeight,Param)
  1558. %% Preallocate ChannelPairName Cell Array
  1559. signifChPairNameTopo = cell(3,3,3,3,2);
  1560. signifChPairNameTopo{3,3,3,3,2} = [];
  1561. % WPLIBResults = struct();
  1562. %% FC parameters for plotting; (may contain p-value variable)
  1563. clear FCpara
  1564. FCpara.ColorMap=colorMapPN; %%%Color map for correlation link
  1565. %% Create orange indigo colormap
  1566. % Number of colors in the colormap
  1567. n = 64;
  1568. % Define orange and indigo RGB values
  1569. orange = [1, 0.75, 0];
  1570. indigo = [0.2, 0.1, 1];
  1571. % Create a colormap by interpolating between orange and indigo
  1572. custom_cmap = [linspace(indigo(1), orange(1), n)', ...
  1573. linspace(indigo(2), orange(2), n)', ...
  1574. linspace(indigo(3), orange(3), n)'];
  1575. % Apply the custom colormap
  1576. % colormap(custom_cmap);
  1577. FCpara.ColorMap=custom_cmap; %%%Color map for correlation link
  1578. %% Other parameters
  1579. FCpara.NodeColor=[0.8 0.8 0.8]; %%%Node Color of Nodes for FC, not important, it is actually defined in ChanPos
  1580. FCpara.Clim=[-1 1]; %%%Color Limit FC
  1581. % FCpara.MarkerSize=8; %%%MarkerSize of scatter
  1582. % FCpara.EdgeColor=[1 0 0]; %%% This is not needed as ColorMap field and Clim field would determine the edge color
  1583. FCpara.EdgeTh=0.1; %%%
  1584. FCpara.NodeTh=0.05; %%%
  1585. FCparaT=FCpara; %%%%T test parameters
  1586. FCparaT.EdgeTh=3; %%%
  1587. FCparaT.NodeTh=0.001; %%%
  1588. FCspear=FCpara; %%%%Spearman r parameters
  1589. FCspear.EdgeTh=0.1; %%% This is the plotting threshold
  1590. FCspear.NodeTh=0.001; %%%
  1591. WPLIEdgeTh=0.6;
  1592. WPLINodeTh=0.05;
  1593. TEdgeTh=1;
  1594. TNodeTh=0.001;
  1595. pTEdgeTh=0.05;
  1596. pSpearEdgeTh=0.05; %8/8/24 0.05 to 0.1
  1597. rSpearEdgeTh=0.1;
  1598. rSpearNodeTh=0.001;
  1599. ScaleWPLI=0.5;
  1600. ScaleT=1;
  1601. ScaleSpear=0.2;
  1602. % BOI=[5;15];
  1603. BOI=[1 4 8 13 30 39.5 43;4 8 13 30 37 41.5 100];
  1604. BandName={'Delta','Theta','Alpha','Beta','Gamma-1','Gamma-E','Gamma-2'};
  1605. BandHzName={'1-4 Hz','4-8 Hz','8-13Hz','13-30 Hz','30-37 Hz','39-41 Hz','43-100 Hz'};
  1606. %% FC calculation - WPLI-Behavior - Spearman Topo Plots *** (old fig4c) fig4d fig4e fig4f
  1607. FCpara.EdgeTh=0.1; %%%
  1608. BOI=[2 4 8 8 10 13 30 39 43;4 8 13 10 13 30 37 41 100]; %[1 4 8 8 10 13 30 39.5 43;4 8 13 10 13 30 37 41.5 100];
  1609. BandName={'Delta','Theta','Alpha','LowAlpha','HighAlpha','Beta','Gamma-1','Gamma-E','Gamma-2'};
  1610. BandHzName={'2-4 Hz','4-8 Hz','8-13Hz','8-10Hz','10-13Hz','13-30 Hz','30-37 Hz','39-41 Hz','43-100 Hz'};
  1611. for iFreqFunc=[3] %1:length(FreqFuncNames) % mean = 1, median =2, peak = 3
  1612. todayDate = datestr(now, 'yymmdd');
  1613. for iTrialType=3 %1:length(TrialType) % 1=Hit, 2=Miss, 3=HitANDMiss
  1614. %% Create save folder
  1615. SaveTemp=[SubSaveWPLI TrialTypeName{iTrialType} '\'];
  1616. SaveTemp=[SaveTemp todayDate '_' FreqFuncNames{iFreqFunc} '_p' num2str(pSpearEdgeTh) '_colored\' ];
  1617. mkdir(SaveTemp)
  1618. %% WPLI figure - ALL groups (6 groups) & ALL BOI - *** (old fig4c) fig4d Complete version - see MS version below
  1619. % Fig 4 C and D are composed of panels created by this section, cut and pasted together in
  1620. % illustrator
  1621. % FCpara.EdgeTh=0.06;
  1622. % FCpara.EdgeTh=0.15; % previous arbitray threshold
  1623. FCpara.EdgeTh=0.1202; % top quartile threshold
  1624. % FCpara.EdgeTh=averagePermutatedWPLI_2to55_top0_0001; % = 0.0894 - top 0.0001 permutation threshold
  1625. % FCpara.EdgeTh=averagePermutatedWPLI_2to55_top0_001; % = 0.0484 - top 0.001 permutation threshold
  1626. FCgroups = [1,3,6]; % [1,3,6] = [40hz , Random, LightRT]
  1627. figure;
  1628. nGroups = length(FCgroups); % (1:3) Just first 3 groups
  1629. for iBOI=1:size(BOI,2)-1
  1630. NeedI=find(Fplot>=BOI(1,iBOI)&Fplot<BOI(2,iBOI));
  1631. for iStimGroup=1:length(FCgroups)
  1632. if iFreqFunc==3
  1633. Tdata{iStimGroup,iBOI,iTrialType}=squeeze(FreqFunc{iFreqFunc}(WPLIall(SubjG{FCgroups(iStimGroup)},NeedI,EEGchInd,EEGchInd,iTrialType),[],2));
  1634. MapGroup{iStimGroup,iBOI,iTrialType}=squeeze(nanmean(FreqFunc{iFreqFunc}(WPLIall(SubjG{FCgroups(iStimGroup)},NeedI,EEGchInd,EEGchInd,iTrialType),[],2),1));
  1635. else
  1636. Tdata{iStimGroup,iBOI,iTrialType}=squeeze(FreqFunc{iFreqFunc}(WPLIall(SubjG{FCgroups(iStimGroup)},NeedI,EEGchInd,EEGchInd,iTrialType),2));
  1637. MapGroup{iStimGroup,iBOI,iTrialType}=squeeze(nanmean(FreqFunc{iFreqFunc}(WPLIall(SubjG{FCgroups(iStimGroup)},NeedI,EEGchInd,EEGchInd,iTrialType),2),1));
  1638. end
  1639. EdgeColor=[0.8 0.8 0.8];
  1640. % FC_BrainEEGLu(ChanPosColin27,MapGroup{iG,iFF,iCom},[],EdgeTh,NodeTh,EdgeColor,[])
  1641. % axis off
  1642. subplotLU(nGroups,size(BOI,2),iStimGroup,iBOI);
  1643. % WPLIsForPlot(iStimGroup,iBOI,iTrialType) = MapGroup{iStimGroup,iBOI,iTrialType};
  1644. % nChPairAboveThreshold(iStimGroup,iBOI,iTrialType) = sum(WPLIsForPlot>FCpara.EdgeTh);
  1645. % Calculate the number of channel pairs with WPLI > EdgeTh
  1646. currentMap = MapGroup{iStimGroup, iBOI, iTrialType};
  1647. nChPairAboveThreshold = sum(currentMap(:) > FCpara.EdgeTh);
  1648. AboveThreholdChannelPairAllGroups(iStimGroup,iBOI,iTrialType) = nChPairAboveThreshold;
  1649. %%% The plotting function
  1650. FC_BrainEEGLu(ChanPosColin27,MapGroup{iStimGroup,iBOI,iTrialType},[],FCpara)
  1651. axis off
  1652. % Add text below the plot
  1653. text(0.5, -0.15, sprintf('Pairs > Th: %d', nChPairAboveThreshold), ...
  1654. 'Units', 'normalized', 'HorizontalAlignment', 'center', 'FontSize', 10);
  1655. if iStimGroup==nGroups
  1656. xlabel(BandName{iBOI});
  1657. text(-0.1,0,0.1,[BandName{iBOI} ' (' BandHzName{iBOI} ')'],'horizontalalignment','center','verticalalignment','top','fontsize',10)
  1658. end
  1659. if iBOI==1
  1660. ylabel(GroupName{FCgroups(iStimGroup)})
  1661. yt=text(0,0.1,0.1,GroupName{FCgroups(iStimGroup)},'horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  1662. end
  1663. end
  1664. end
  1665. % Define the edge threshold title
  1666. edgeThresholdTitle = sprintf('Edge Threshold: %.4f', FCpara.EdgeTh);
  1667. % Add the title displaying the edge threshold
  1668. sgtitle(edgeThresholdTitle, 'FontSize', 6, 'FontWeight', 'bold', 'Interpreter', 'none');
  1669. % % Add the title displaying the edge threshold and position it higher
  1670. % titleHandle = sgtitle(edgeThresholdTitle, 'FontSize', 12);
  1671. % titleHandle.Position = [0.5, 0.98, 0]; % [x, y, z] position in normalized figure units
  1672. % subplot('position',[0.5 0.51 0.3 0.01]);
  1673. % b=colorbar('southoutside');
  1674. % set(gca,'xtick',[],'ytick',[])
  1675. % set(b,'position',[0.5 0.5 0.3 0.03],'Limits',[0 1],'Ticks',[0 1],'Ticklabels',PowerLab);
  1676. % xlabel(b,'Log Normalized Power')
  1677. LuFontStandard;
  1678. papersizePX=[0 0 6*size(BOI,2) 6*nGroups+3];
  1679. set(gcf, 'PaperUnits', 'centimeters');
  1680. set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  1681. % Adjust title spacing relative to the subplots
  1682. t = sgtitle(edgeThresholdTitle, 'FontSize', 6);
  1683. % t.Position(2) = t.Position(2) + 0.03; % Raise the title slightly
  1684. saveWPLIFigName = ['WPLI40_' num2str(nGroups) 'Groups_' num2str(100*FCpara.EdgeTh) 'E-2EdgeThrs620.png'];
  1685. % saveas(gcf,[SaveTemp saveWPLIFigName],'pdf');
  1686. saveas(gcf,[SaveTemp saveWPLIFigName],'png');
  1687. saveas(gcf,[SaveTemp saveWPLIFigName],'svg');
  1688. % saveas(gcf,[SaveTemp saveWPLIFigName],'epsc');
  1689. FCpara.EdgeTh=0.1; % reset
  1690. %% WPLI figure *** fig4d MS version
  1691. % Fig 4 C and D are composed of panels created by this section, cut and pasted together in
  1692. % illustrator
  1693. % Lower Alpha and Upper Alpha only
  1694. BOI=[8 10; 10 13]; %[1 4 8 8 10 13 30 39.5 43;4 8 13 10 13 30 37 41.5 100];
  1695. BandName={'LowerAlpha','UpperAlpha'};
  1696. BandHzName={'8-10Hz','10-13Hz'};
  1697. % FCpara.EdgeTh=0.06;
  1698. % FCpara.EdgeTh=0.15; % previous arbitray threshold
  1699. FCpara.EdgeTh=0.1202; % top quartile threshold
  1700. % FCpara.EdgeTh=averagePermutatedWPLI_2to55_top0_0001; % = 0.0894 - top 0.0001 permutation threshold
  1701. % FCpara.EdgeTh=averagePermutatedWPLI_2to55_top0_001; % = 0.0484 - top 0.001 permutation threshold
  1702. FCgroups = [1,3,6]; % [1,3,6] = [40hz , Random, LightRT]
  1703. figure;
  1704. nGroups = length(FCgroups); % (1:3) Just first 3 groups
  1705. for iBOI=1:size(BOI,2) % just lower alpha and upper alpha
  1706. NeedI=find(Fplot>=BOI(1,iBOI)&Fplot<BOI(2,iBOI));
  1707. for iStimGroup=1:length(FCgroups)
  1708. if iFreqFunc==3
  1709. Tdata{iStimGroup,iBOI,iTrialType}=squeeze(FreqFunc{iFreqFunc}(WPLIall(SubjG{FCgroups(iStimGroup)},NeedI,EEGchInd,EEGchInd,iTrialType),[],2));
  1710. MapGroup{iStimGroup,iBOI,iTrialType}=squeeze(nanmean(FreqFunc{iFreqFunc}(WPLIall(SubjG{FCgroups(iStimGroup)},NeedI,EEGchInd,EEGchInd,iTrialType),[],2),1));
  1711. else
  1712. Tdata{iStimGroup,iBOI,iTrialType}=squeeze(FreqFunc{iFreqFunc}(WPLIall(SubjG{FCgroups(iStimGroup)},NeedI,EEGchInd,EEGchInd,iTrialType),2));
  1713. MapGroup{iStimGroup,iBOI,iTrialType}=squeeze(nanmean(FreqFunc{iFreqFunc}(WPLIall(SubjG{FCgroups(iStimGroup)},NeedI,EEGchInd,EEGchInd,iTrialType),2),1));
  1714. end
  1715. EdgeColor=[0.8 0.8 0.8];
  1716. % FC_BrainEEGLu(ChanPosColin27,MapGroup{iG,iFF,iCom},[],EdgeTh,NodeTh,EdgeColor,[])
  1717. % axis off
  1718. % subplotLU(nGroups,size(BOI,2),iStimGroup,iBOI); % old 3x2
  1719. % subplotLU(1, nGroups*size(BOI,2), 1, iStimGroup+3*(iBOI-1)); % 1x6 grid
  1720. subplotLU(size(BOI,2), nGroups, iBOI, iStimGroup); % 2x3 grid
  1721. % WPLIsForPlot(iStimGroup,iBOI,iTrialType) = MapGroup{iStimGroup,iBOI,iTrialType};
  1722. % nChPairAboveThreshold(iStimGroup,iBOI,iTrialType) = sum(WPLIsForPlot>FCpara.EdgeTh);
  1723. % Calculate the number of channel pairs with WPLI > EdgeTh
  1724. currentMap = MapGroup{iStimGroup, iBOI, iTrialType};
  1725. nChPairAboveThreshold = sum(currentMap(:) > FCpara.EdgeTh);
  1726. AboveThreholdChannelPairAllGroups(iStimGroup,iBOI,iTrialType) = nChPairAboveThreshold;
  1727. %%% The plotting function
  1728. FC_BrainEEGLu(ChanPosColin27,MapGroup{iStimGroup,iBOI,iTrialType},[],FCpara)
  1729. axis off
  1730. % Add text below the plot
  1731. % text(0.5, -0.15, sprintf('Pairs > Th: %d', nChPairAboveThreshold), ...
  1732. % 'Units', 'normalized', 'HorizontalAlignment', 'center', 'FontSize', 10);
  1733. if iStimGroup==2
  1734. xlabel(BandName{iBOI});
  1735. text(0.125,0,0.1,[BandName{iBOI} ' (' BandHzName{iBOI} ')'],'horizontalalignment','center','verticalalignment','top','fontsize',18)
  1736. end
  1737. xlabel(GroupName{FCgroups(iStimGroup)})
  1738. yt=text(-0.125,0,0.1,GroupName{FCgroups(iStimGroup)},'horizontalalignment','center','verticalalignment','bottom','fontsize',14);
  1739. end
  1740. end
  1741. % Define the edge threshold title
  1742. % edgeThresholdTitle = sprintf('Edge Threshold: %.4f', FCpara.EdgeTh);
  1743. % % Add the title displaying the edge threshold
  1744. % sgtitle(edgeThresholdTitle, 'FontSize', 6, 'FontWeight', 'bold', 'Interpreter', 'none');
  1745. % t.Position(2) = t.Position(2) + 0.05; % Move title slightly up
  1746. % % Add the title displaying the edge threshold and position it higher
  1747. % titleHandle = sgtitle(edgeThresholdTitle, 'FontSize', 12);
  1748. % titleHandle.Position = [0.5, 0.98, 0]; % [x, y, z] position in normalized figure units
  1749. % subplot('position',[0.5 0.51 0.3 0.01]);
  1750. % b=colorbar('southoutside');
  1751. % set(gca,'xtick',[],'ytick',[])
  1752. % set(b,'position',[0.5 0.5 0.3 0.03],'Limits',[0 1],'Ticks',[0 1],'Ticklabels',PowerLab);
  1753. % xlabel(b,'Log Normalized Power')
  1754. LuFontStandard;
  1755. % papersizePX=[0 0 6*size(BOI,2)*4 6/3*nGroups+3]; % 1x6
  1756. % 3x2 grid: papersizePX=[0 0 6*size(BOI,2) 6*nGroups+3];
  1757. papersizePX=[0 0 nGroups*6 size(BOI,2)*6+2.5]; % 2x3 : [0 0 width{x} heigth{y}]
  1758. set(gcf, 'PaperUnits', 'centimeters');
  1759. set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  1760. % Adjust title spacing relative to the subplots
  1761. % t = sgtitle(edgeThresholdTitle, 'FontSize', 6);
  1762. % t.Position(2) = t.Position(2) + 0.03; % Raise the title slightly
  1763. % saveas(gcf,[SaveTemp saveWPLIFigName],'pdf');
  1764. saveDateFig4D = datestr(datetime, 'yy-mm-dd_HHMMSSFFF');
  1765. saveWPLIFigName = ['WPLI40_' num2str(nGroups) 'Groups_' saveDateFig4D '_' num2str(100*FCpara.EdgeTh) 'E-2EdgeThrs620.png'];
  1766. saveas(gcf,['Fig4Panels/' saveWPLIFigName ],'png');
  1767. % saveas(gcf,['Fig4Panels/' saveWPLIFigName],'svg');
  1768. % saveas(gcf,[SaveTemp saveWPLIFigName],'epsc');
  1769. FCpara.EdgeTh=0.1; % reset
  1770. %% Perform Chi-Squared on proportions on elevated channels: 40vL & 40VRand
  1771. % Data
  1772. total_pairs = 496;
  1773. elevated_40Hz = 93;
  1774. elevated_Light = 32;
  1775. % Define contingency table
  1776. table_40Hz_Light = [elevated_40Hz, total_pairs - elevated_40Hz;
  1777. elevated_Light, total_pairs - elevated_Light];
  1778. % % Perform chi-squared test and compute effect size
  1779. % [chi2_40Hz_Light, p_40Hz_Light, V_40Hz_Light] = analyzeChiSquared(table_40Hz_Light, total_pairs);
  1780. %
  1781. % % Display results
  1782. % fprintf('Results for 40Hz vs Light:\n');
  1783. % fprintf('Chi-squared (X^2): %.2f\n', chi2_40Hz_Light);
  1784. % fprintf('p-value: %.4f\n', p_40Hz_Light);
  1785. % fprintf('Cramér''s V (Effect size): %.4f\n', V_40Hz_Light);
  1786. %% Chi Squared Tests: 40vLight, 40vRandom for Lower and Upper Alpha Fig4D stats - OLD
  1787. % Contingency tables for the tests
  1788. % Lower Alpha:
  1789. % 40Hz vs Light
  1790. observed_LowerAlpha_40vLight = [93, 403;
  1791. 32, 464];
  1792. % Lower Alpha: 40Hz vs Random
  1793. observed_LowerAlpha_40vRandom = [93, 403;
  1794. 7, 489];
  1795. % Upper Alpha:
  1796. % 40Hz vs Light
  1797. observed_UpperAlpha_40vLight = [63, 433;
  1798. 7, 489];
  1799. % Upper Alpha: 40Hz vs Random
  1800. observed_UpperAlpha_40vRandom = [63, 433;
  1801. 267, 229];
  1802. % Perform chi-squared tests
  1803. fprintf('Lower Alpha (40Hz vs Light):\n');
  1804. [chi2_LA_40vLight, p_LA_40vLight, V_LA_40vLight, dof_LA_40vLight, total_LA_40vLight] = chi_squared_test(observed_LowerAlpha_40vLight);
  1805. fprintf('\nLower Alpha (40Hz vs Random):\n');
  1806. [chi2_LA_40vRandom, p_LA_40vRandom, V_LA_40vRandom, dof_LA_40vRandom, total_LA_40vRandom] = chi_squared_test(observed_LowerAlpha_40vRandom);
  1807. fprintf('\nUpper Alpha (40Hz vs Light):\n');
  1808. [chi2_UA_40vLight, p_UA_40vLight, V_UA_40vLight, dof_UA_40vLight, total_UA_40vLight] = chi_squared_test(observed_UpperAlpha_40vLight);
  1809. fprintf('\nUpper Alpha (40Hz vs Random):\n');
  1810. [chi2_UA_40vRandom, p_UA_40vRandom, V_UA_40vRandom, dof_UA_40vRandom, total_UA_40vRandom] = chi_squared_test(observed_UpperAlpha_40vRandom);
  1811. % FDR correction:
  1812. % Lower Alpha p-values
  1813. p_values_LowerAlpha = [p_LA_40vLight, p_LA_40vRandom];
  1814. % Upper Alpha p-values
  1815. p_values_UpperAlpha = [p_UA_40vLight, p_UA_40vRandom];
  1816. % Apply FDR correction to Lower Alpha and Upper Alpha
  1817. adjusted_p_LowerAlpha = fdr_correction(p_values_LowerAlpha);
  1818. adjusted_p_UpperAlpha = fdr_correction(p_values_UpperAlpha);
  1819. % Display results
  1820. disp('FDR-corrected p-values for Lower Alpha:');
  1821. disp(adjusted_p_LowerAlpha);
  1822. disp('FDR-corrected p-values for Upper Alpha:');
  1823. disp(adjusted_p_UpperAlpha);
  1824. % Display FDR-corrected p-values for Lower Alpha in scientific notation
  1825. disp('FDR-corrected p-values for Lower Alpha (scientific notation):');
  1826. fprintf('%.15e\n', adjusted_p_LowerAlpha);
  1827. % Display FDR-corrected p-values for Upper Alpha in scientific notation
  1828. disp('FDR-corrected p-values for Upper Alpha (scientific notation):');
  1829. fprintf('%.15e\n', adjusted_p_UpperAlpha);
  1830. %% Chi-squared test of top quartile ch pairs & permutated channel pairs - no longer needed 1/16/24
  1831. % % Assuming AboveThresholdChannelPairAllGroupsTopQuart and
  1832. % % AboveThresholdChannelPairAllGroupsPermutated are already loaded.
  1833. %
  1834. % % Set denominator for proportions
  1835. % denominator = 496;
  1836. %
  1837. % % Extract dimensions
  1838. % dims = size(AboveThreholdChannelPairAllGroupsTopQuart);
  1839. %
  1840. % % Initialize matrices for storing results
  1841. % chi2_stat = zeros(dims); % Chi-squared statistic
  1842. % p_value = zeros(dims); % P-value
  1843. % h_test = zeros(dims); % Hypothesis test result (1: reject null, 0: fail to reject)
  1844. %
  1845. % % Loop through each element
  1846. % for i = 1:dims(1)
  1847. % for j = 1:dims(2)
  1848. % for k = 1:dims(3)
  1849. % % Observed data
  1850. % obs1 = AboveThreholdChannelPairAllGroupsTopQuart(i,j,k);
  1851. % obs2 = AboveThreholdChannelPairAllGroupsPermutated(i,j,k);
  1852. %
  1853. % % Proportions
  1854. % prop1 = obs1 / denominator;
  1855. % prop2 = obs2 / denominator;
  1856. %
  1857. % % Pooled proportion under null hypothesis
  1858. % pooled_p = (obs1 + obs2) / (2 * denominator);
  1859. %
  1860. % % Expected counts under null hypothesis
  1861. % exp1 = pooled_p * denominator;
  1862. % exp2 = pooled_p * denominator;
  1863. %
  1864. % % Chi-squared statistic for this pair
  1865. % chi2_stat(i,j,k) = ((obs1 - exp1)^2 / exp1) + ((obs2 - exp2)^2 / exp2);
  1866. %
  1867. % % Degrees of freedom
  1868. % df = 1;
  1869. %
  1870. % % Compute p-value
  1871. % p_value(i,j,k) = 1 - chi2cdf(chi2_stat(i,j,k), df);
  1872. %
  1873. % % Hypothesis test: reject null if p < 0.05
  1874. % h_test(i,j,k) = p_value(i,j,k) < 0.05;
  1875. % end
  1876. % end
  1877. % end
  1878. %
  1879. % % Display results for inspection
  1880. % disp('Chi-squared statistics:');
  1881. % disp(chi2_stat);
  1882. %
  1883. % disp('P-values (3D matrix):');
  1884. % disp(p_value);
  1885. %
  1886. % disp('Hypothesis test results (3D matrix, 1: reject null, 0: fail to reject):');
  1887. % disp(h_test);
  1888. %% WPLI figure (Both Controls)
  1889. % figure;
  1890. % for iBOI=1:size(BOI,2)
  1891. % NeedI=find(Fplot>=BOI(1,iBOI)&Fplot<BOI(2,iBOI));
  1892. % for iStimGroup=1:length(SubjG)
  1893. % if iFreqFunc==3
  1894. % Tdata{iStimGroup,iBOI,iTrialType}=squeeze(FreqFunc{iFreqFunc}(WPLIall(SubjG{iStimGroup},NeedI,EEGchInd,EEGchInd,iTrialType),[],2));
  1895. % MapGroup{iStimGroup,iBOI,iTrialType}=squeeze(nanmean(FreqFunc{iFreqFunc}(WPLIall(SubjG{iStimGroup},NeedI,EEGchInd,EEGchInd,iTrialType),[],2),1));
  1896. % else
  1897. % Tdata{iStimGroup,iBOI,iTrialType}=squeeze(FreqFunc{iFreqFunc}(WPLIall(SubjG{iStimGroup},NeedI,EEGchInd,EEGchInd,iTrialType),2));
  1898. % MapGroup{iStimGroup,iBOI,iTrialType}=squeeze(nanmean(FreqFunc{iFreqFunc}(WPLIall(SubjG{iStimGroup},NeedI,EEGchInd,EEGchInd,iTrialType),2),1));
  1899. % end
  1900. %
  1901. % EdgeColor=[0.8 0.8 0.8];
  1902. % % FC_BrainEEGLu(ChanPosColin27,MapGroup{iG,iFF,iCom},[],EdgeTh,NodeTh,EdgeColor,[])
  1903. % % axis off
  1904. %
  1905. % subplotLU(2,size(BOI,2),iStimGroup,iBOI);
  1906. % FC_BrainEEGLu(ChanPosColin27,MapGroup{iStimGroup,iBOI,iTrialType},[],FCpara)
  1907. % axis off
  1908. % if iStimGroup==2
  1909. % xlabel(BandName{iBOI});
  1910. % text(-0.1,0,0.1,[BandName{iBOI} ' (' BandHzName{iBOI} ')'],'horizontalalignment','center','verticalalignment','top','fontsize',10)
  1911. % end
  1912. %
  1913. % if iBOI==1
  1914. % ylabel(GroupName{iStimGroup})
  1915. % yt=text(0,0.1,0.1,GroupName{iStimGroup},'horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  1916. % end
  1917. % end
  1918. % end
  1919. % % subplot('position',[0.5 0.51 0.3 0.01]);
  1920. % % b=colorbar('southoutside');
  1921. % % set(gca,'xtick',[],'ytick',[])
  1922. % % set(b,'position',[0.5 0.5 0.3 0.03],'Limits',[0 1],'Ticks',[0 1],'Ticklabels',PowerLab);
  1923. % % xlabel(b,'Log Normalized Power')
  1924. % LuFontStandard;
  1925. % papersizePX=[0 0 6*size(BOI,2) 6*2+3];
  1926. % set(gcf, 'PaperUnits', 'centimeters');
  1927. % set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  1928. %
  1929. % WPLI40BothControlsFigName = 'WPLI40BothControls620';
  1930. % % saveas(gcf,[SaveTemp WPLI40BothControlsFigName],'pdf');
  1931. % saveas(gcf,[SaveTemp WPLI40BothControlsFigName],'png');
  1932. % % saveas(gcf,[SaveTemp WPLI40BothControlsFigName],'epsc');
  1933. %% WPLI Difference Figure (VARIABLES NEED TO BE RENAMED IN THIS SECTION)
  1934. % figure;
  1935. % pTEdgeTh=0.1;
  1936. % %comparison between groups w/ freq data
  1937. % diffTPmap=zeros(length(ChanEEGLab),length(ChanEEGLab),size(BOI,2),length(TrialType)); % for T test
  1938. % diffRPmap=diffTPmap;
  1939. % diffTmap=diffTPmap;
  1940. %
  1941. % % Group selection: 1=40Hz, 2=Light, 3=Random, 4=LightRT
  1942. % Group1 = 1;
  1943. % Group2 = 2;
  1944. %
  1945. % for iBOI=1:size(BOI,2)-1 % minus 1 to remove BOI with 60 Hz
  1946. % for iCh=1:length(ChanEEGLab)
  1947. % for jCh=iCh+1:length(ChanEEGLab)
  1948. % [~,diffTPmap(iCh,jCh,iBOI,iTrialType),~,stats]=ttest2(Tdata{Group1,iBOI,iTrialType}(:,iCh,jCh),Tdata{Group2,iBOI,iTrialType}(:,iCh,jCh));
  1949. % [diffRPmap(iCh,jCh,iBOI,iTrialType),~,~]=ranksum(Tdata{Group1,iBOI,iTrialType}(:,iCh,jCh),Tdata{Group2,iBOI,iTrialType}(:,iCh,jCh));
  1950. % diffTmap(iCh,jCh,iBOI,iTrialType)=stats.tstat;
  1951. % end
  1952. % end
  1953. % subplotLU(1,size(BOI,2),1,iBOI);
  1954. % Adj=diffTmap(:,:,iBOI,iTrialType);
  1955. % AdjP=diffTPmap(:,:,iBOI,iTrialType);
  1956. % Adj(AdjP>pTEdgeTh)=0;
  1957. % FC_BrainEEGLu(ChanPosColin27,Adj,[],FCparaT)
  1958. % axis off
  1959. %
  1960. % if iBOI==1
  1961. % yt=text(-0.0,0.1,0.1,['T, Sig-Diff FC, ' GroupName{Group1} '-' GroupName{Group2} ],'horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  1962. % end
  1963. %
  1964. % % % subplotLU(2,size(BOI,2),2,iFF);
  1965. % % % xlabel(BName{iFF});
  1966. % % % text(0,-0.55,[BName{iFF} ' (' BName2{iFF} ')'],'horizontalalignment','center','verticalalignment','top','fontsize',10)
  1967. % % % if iFF==1
  1968. % % % a=ylabel('40Hz-Random');
  1969. % % % % a.Position=[0.01 0.5 0.03 0.4];
  1970. % % % % a.verticalalignment='middle';
  1971. % % % % set(a,'Position',[0.01 0.5 0.03 0.4],'Verticalalignment','middle')
  1972. % % % set(a,'Verticalalignment','middle')
  1973. % % % yt=text(-120,0,'P, 40-Rand.','horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  1974. % % % end
  1975. %
  1976. % % xlabel(BName{iFF});
  1977. % text(-0.1,0,0.1,[BandName{iBOI} ' (' BandHzName{iBOI} ')'],'horizontalalignment','center','verticalalignment','top','fontsize',10)
  1978. % end
  1979. %
  1980. % LuFontStandard;
  1981. % papersizePX=[0 0 6*size(BOI,2) 6+2];
  1982. % set(gcf, 'PaperUnits', 'centimeters');
  1983. % set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  1984. % WPLIDiffFigName = ['WPLIDiff ' GroupName{Group1} '-' GroupName{Group2} ' p' num2str(pTEdgeTh)];
  1985. % sgtitle(WPLIDiffFigName)
  1986. %
  1987. %
  1988. % % saveas(gcf,[SaveTemp WPLIDiffBothControlsfigName],'pdf');
  1989. % saveas(gcf,[SaveTemp WPLIDiffFigName '.png'],'png');
  1990. % % saveas(gcf,[SaveTemp WPLIDiffBothControlsfigName],'epsc');
  1991. %% Define Group Set Names - beginning of fig4e fig4f
  1992. GroupSetsName{1} ='40andL'; % 40 and Light
  1993. GroupSetsName{2}='40andR'; % 40 and Random
  1994. GroupSetsName{3}='All3Groups'; % All three groups
  1995. %% Loop thru all groups *** fig4e fig4f
  1996. pSpearEdgeTh = 0.1;
  1997. for iGroupSet = [1 2]% 1:length(GroupSetsName)
  1998. if iGroupSet == 1
  1999. %% 40&LightRT
  2000. % included subs
  2001. Group1 = 1; % 1= 40Hz flicker group
  2002. Group2 = 6; % 2 = LightRT group
  2003. dataName = [TrialTypeName{iTrialType} ' ' FreqFuncNames{iFreqFunc} ' p' num2str(pSpearEdgeTh) ' ' GroupName{Group1} GroupName{Group2}];
  2004. IncludedSubj=union(SubjG{Group1},SubjG{Group2}); %40 + LightRT
  2005. elseif iGroupSet == 2
  2006. %% 40andR
  2007. Group1 = 1; % 1= 40Hz flicker group
  2008. Group2 = 3; % 3 = Random group
  2009. dataName = [TrialTypeName{iTrialType} ' ' FreqFuncNames{iFreqFunc} ' p' num2str(pSpearEdgeTh) ' ' GroupName{Group1} GroupName{Group2}];
  2010. IncludedSubj=union(SubjG{Group1},SubjG{Group2}); %40 + Random
  2011. elseif iGroupSet == 3
  2012. %% All3Groups
  2013. Group1 = 1; % 1= 40Hz flicker group
  2014. Group2 = 6; % 6 = LightRT
  2015. Group3 = 3; % 3 = Random group
  2016. dataName = [TrialTypeName{iTrialType} FreqFuncNames{iFreqFunc} GroupName{Group1} GroupName{Group2} GroupName{Group3}];
  2017. IncludedSubj=union(SubjG{Group1},SubjG{Group2},SubjG{Group3}); %40 + LightRT + Random
  2018. end
  2019. %% Topo plot ALPHA GENERAL
  2020. %% Define Bands of Interest
  2021. % BOI=[8 8 10; 13 10 13];
  2022. % BandName={'Alpha','LowAlpha','HighAlpha'};
  2023. % BandHzName={'8-13Hz','8-10Hz','10-13Hz'};
  2024. BOI=[ 8 10; 10 13];
  2025. BandName={'LowAlpha','HighAlpha'};
  2026. BandHzName={'8-10Hz','10-13Hz'};
  2027. %% Spearman Acc-FC
  2028. SpearmanAccFctopo = figure;
  2029. for iBOI=1:size(BOI,2)
  2030. NeedI=find(Fplot>=BOI(1,iBOI)&Fplot<=BOI(2,iBOI));
  2031. if iFreqFunc==3
  2032. WPLItemp=squeeze(FreqFunc{iFreqFunc}(WPLIall(IncludedSubj,NeedI,EEGchInd,EEGchInd,iTrialType),[],2));
  2033. else
  2034. WPLItemp=squeeze(FreqFunc{iFreqFunc}(WPLIall(IncludedSubj,NeedI,EEGchInd,EEGchInd,iTrialType),2));
  2035. end
  2036. Acctemp=Acc(IncludedSubj);
  2037. %% Include only subjects with high accuracy (>80%) - Accuracy cuts
  2038. highAccThreshold = 0.8;
  2039. highAccSubjectsIndex = Acctemp>highAccThreshold;
  2040. WPLItemp = WPLItemp(highAccSubjectsIndex,:,:);
  2041. Acctemp80 = Acctemp(highAccSubjectsIndex);
  2042. %% loop thru each channel - calculate Spearman correlation R and p-value
  2043. for iCh=1:length(ChanEEGLab)-1
  2044. % for jCh=2:length(ChanEEGLab)
  2045. [rSpear(iCh,iCh+1:end,iBOI,iTrialType),pSpear(iCh,iCh+1:end,iBOI,iTrialType)]=corr(squeeze(WPLItemp(:,iCh,iCh+1:end)),Acctemp80,'type','spearman','rows','pairwise');
  2046. % end
  2047. end
  2048. subplotLU(1,size(BOI,2),1,iBOI); %can change to 2 for second row
  2049. Adj=rSpear(:,:,iBOI,iTrialType); %32x32 of r values?
  2050. AdjP=pSpear(:,:,iBOI,iTrialType); %32x32 of p values?
  2051. Adj(AdjP>pSpearEdgeTh)=0; % deletes all r values of ch-pairs with p-value > threshold
  2052. signifchPairRAccTopo = Adj~=0; % creates 32x32 logical of significant channel pairs
  2053. signifChPairNameTopo{iFreqFunc,iTrialType,iGroupSet,iBOI,1} = getListofSignifChPairs(signifchPairRAccTopo,EEGch); % get list of channel names not equal to zero
  2054. %% Now, plot only for significant positive WPLI-accuracy correlations
  2055. % posWPLIAccChPairFolder = ['PositiveAccCorr\' GroupSetsName{iGroupSet} '\' BName{iBOI} '\'];
  2056. % mkdir([SaveTemp posWPLIAccChPairFolder])
  2057. % for iCh = 1:length(EEGchInd)
  2058. % for jCh = iCh+1:length(EEGchInd)
  2059. % % Check if the channel pair has a positive significant correlation
  2060. % if Adj(iCh, jCh) > 0
  2061. % clear DataPlot;
  2062. %
  2063. % for iStimGroup = 1:length(SubjGWPLI)
  2064. % DataPlot{iStimGroup} = squeeze(WPLIall(SubjGWPLI{iStimGroup}, :, EEGchInd(iCh), EEGchInd(jCh), iTrialType));
  2065. % Invalid = isnan(DataPlot{iStimGroup}(:, 1));
  2066. % DataPlot{iStimGroup}(Invalid, :) = [];
  2067. % end
  2068. %
  2069. % if isempty(DataPlot{1}) || isempty(DataPlot{2})
  2070. % continue;
  2071. % end
  2072. %
  2073. % iPlot = iPlot + 1;
  2074. %
  2075. % % Plot WPLI data for the channel pair
  2076. % figure;
  2077. % xrange = [0 50]; % Frequency range for x-axis
  2078. % subplot('Position', [0.1 0.1 0.88 0.88]);
  2079. % RateHist_GroupPlot(Fplot, DataPlot, FlickerColor, ParamWPLI); % Custom function for plotting
  2080. %
  2081. % text(range(xrange)/2 - 5, ParamWPLI.Ytick(end), [EEGch{iCh} '-' EEGch{jCh}], 'FontSize', 10); % Add channel pair label
  2082. % set(gca, 'xlim', xrange, 'xtick', [1 4 8 13 30 40 50 60 80 100], 'ylim', [0 0.2], 'ytick', ParamWPLI.Ytick);
  2083. % xlabel('Frequency (Hz)');
  2084. % ylabel('WPLI');
  2085. %
  2086. % % Formatting for grid and font size
  2087. % ax = gca;
  2088. % ax.XGrid = 'on';
  2089. % ax.YGrid = 'off';
  2090. % ax.FontSize = 7;
  2091. %
  2092. % % Set paper size and save the figure
  2093. % LuFontStandard; % Custom function for standard fonts
  2094. % papersizePX = 1.3 * [0 0 5 3.2]; % Custom figure size
  2095. % set(gcf, 'PaperUnits', 'centimeters');
  2096. % set(gcf, 'PaperPosition', papersizePX, 'PaperSize', papersizePX(3:4));
  2097. %
  2098. % % Save the figure in SVG format
  2099. % % saveas(gcf, [SaveTemp 'PositiveAccCorr\' '40LightRandBands_' EEGch{iCh} '-' EEGch{jCh}], 'svg');
  2100. % saveas(gcf, [SaveTemp posWPLIAccChPairFolder 'WPLI_40LR_' EEGch{iCh} '-' EEGch{jCh}], 'png');
  2101. % end
  2102. % end
  2103. % end
  2104. % figure(SpearmanAccFctopo);
  2105. %% PLotting and axises ****fig4e
  2106. FC_BrainEEGLu(ChanPosColin27,Adj,[sum(Adj,1)/2],FCspear) % plotting function! (3rd parameter is node weight)
  2107. nSignifChPairs = sum(AdjP<pSpearEdgeTh & AdjP>0, 'all');
  2108. axis off
  2109. if iBOI==1
  2110. % yt=text(0,0.1,0.1,'Sig-Corr. FC-Acc','horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  2111. end
  2112. % text(-0.1,0,0.1,[BandName{iBOI} ' (' BandHzName{iBOI} ')'],'horizontalalignment','center','verticalalignment','top','fontsize',10)
  2113. text(-0.11,0,-0.1,['nChpairs=' num2str(nSignifChPairs)],'horizontalalignment','center','verticalalignment','top','fontsize',10)
  2114. if iBOI==2
  2115. % title(dataName, 'Units', 'normalized', 'Position', [0.5, 0.9, 0]) % MKA 9/23
  2116. % title(dataName, 'Units', 'normalized', 'Position', [0, 0.9, 0])
  2117. end
  2118. %% Display the list of significant channel pairs below the subplot
  2119. % signifChPairNames = getListofSignifChPairs(signifchPairRAccTopo, EEGch);
  2120. % % text(0.5, -0.2, strjoin(signifChPairNames, ', '), 'Units', 'normalized', 'HorizontalAlignment', 'center', 'FontSize', 8);
  2121. % % verticalSignifChPairs = strjoin(signifChPairNames, '\n');
  2122. % % text(-0.11, -0.2, verticalSignifChPairs, 'Units', 'normalized', 'HorizontalAlignment', 'center', 'FontSize', 8);
  2123. %
  2124. % % New subplot for channel pair names (as three columns)
  2125. % subplotLU(2, size(BOI, 2), 2, iBOI); % New row (2nd row) for the names
  2126. %
  2127. % % Divide the list into three columns
  2128. % numNames = length(signifChPairNames);
  2129. % numPerCol = ceil(numNames / 3);
  2130. %
  2131. % % Split the list of significant channel pairs into three columns
  2132. % col1 = signifChPairNames(1:numPerCol);
  2133. % col2 = signifChPairNames(numPerCol+1:min(2*numPerCol, numNames));
  2134. % col3 = signifChPairNames(2*numPerCol+1:end);
  2135. %
  2136. % % Prepare the text to display in columns
  2137. % colText = sprintf('%s\n', col1{:});
  2138. % colText2 = sprintf('%s\n', col2{:});
  2139. % colText3 = sprintf('%s\n', col3{:});
  2140. %
  2141. % % Display the three columns of significant channel pairs
  2142. % text(0.2, 0.5, colText, 'Units', 'normalized', 'HorizontalAlignment', 'left', 'FontSize', 6);
  2143. % text(0.5, 0.5, colText2, 'Units', 'normalized', 'HorizontalAlignment', 'left', 'FontSize', 6);
  2144. % text(0.8, 0.5, colText3, 'Units', 'normalized', 'HorizontalAlignment', 'left', 'FontSize', 6);
  2145. % axis off
  2146. end
  2147. %% Set up and save figure
  2148. LuFontStandard;
  2149. papersizePX=[0 0 6*size(BOI,2) 6+2]; % size of paper - height and width of brain plot 9/10. gain increase the papersize to try to increase the resolution
  2150. set(gcf, 'PaperUnits', 'centimeters');
  2151. set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  2152. SpearmanFCAccName1 = ['SpearmanFCAcc_ALPHA' '_' dataName];
  2153. % saveas(gcf,[SaveTemp SpearmanFCAccName1],'pdf');
  2154. saveas(gcf,[SaveTemp SpearmanFCAccName1 '.png'],'png');
  2155. print(gcf, '-dsvg', [SaveTemp SpearmanFCAccName1 '.svg'] , '-r0'); % -r0 ensures full vector output.
  2156. % saveas(gcf,[SaveTemp SpearmanFCAccName1],'svg'); % -r0 ensures full vector output.);
  2157. % saveas(gcf,[SaveTemp SpearmanFCAccName1 '.eps'],'epsc');
  2158. close all
  2159. %% Spearman RT-FC - Topo plot GENERAL ALPHA ***fig4f
  2160. SpearmanRTFCtopo = figure;
  2161. for iBOI=1:size(BOI,2)
  2162. NeedI=find(Fplot>=BOI(1,iBOI)&Fplot<=BOI(2,iBOI));
  2163. if iFreqFunc==3
  2164. WPLItemp=squeeze(FreqFunc{iFreqFunc}(WPLIall(IncludedSubj,NeedI,EEGchInd,EEGchInd,iTrialType),[],2));
  2165. else
  2166. WPLItemp=squeeze(FreqFunc{iFreqFunc}(WPLIall(IncludedSubj,NeedI,EEGchInd,EEGchInd,iTrialType),2));
  2167. end
  2168. RTtemp=SubjsAvgRT(IncludedSubj);
  2169. %% Include only subjects with high accuracy (>80%) - Accuracy cuts
  2170. highAccThreshold = 0.8;
  2171. highAccSubjectsIndex = Acctemp>highAccThreshold;
  2172. WPLItemp = WPLItemp(highAccSubjectsIndex,:,:);
  2173. RTtemp = RTtemp(highAccSubjectsIndex);
  2174. %%
  2175. for iCh=1:length(ChanEEGLab)-1
  2176. % for jCh=2:length(ChanEEGLab)
  2177. [rSpear(iCh,iCh+1:end,iBOI,iTrialType),pSpear(iCh,iCh+1:end,iBOI,iTrialType)]=corr(squeeze(WPLItemp(:,iCh,iCh+1:end)),RTtemp,'type','spearman','rows','pairwise');
  2178. % end
  2179. end
  2180. subplotLU(1,size(BOI,2),1,iBOI); %can change to 2 for second row
  2181. Adj=rSpear(:,:,iBOI,iTrialType); %32x32 of r values?
  2182. AdjP=pSpear(:,:,iBOI,iTrialType); %32x32 of p values?
  2183. Adj(AdjP>pSpearEdgeTh)=0;
  2184. signifchPairRRTTopo = Adj~=0; % creates 32x32 logical of significant channel pairs
  2185. signifChPairNameTopo{iFreqFunc,iTrialType,iGroupSet,iBOI,2} = getListofSignifChPairs(signifchPairRRTTopo,EEGch); % get list of channel names not equal to zero
  2186. %% Now, plot only for significant negative WPLI-accuracy correlations
  2187. % negWPLIRTChPairFolder = ['NegativeRTCorr\' GroupSetsName{iGroupSet} '\' BName{iBOI} '\'];
  2188. % mkdir([SaveTemp negWPLIRTChPairFolder])
  2189. % for iCh = 1:length(EEGchInd)
  2190. % for jCh = iCh+1:length(EEGchInd)
  2191. % % Check if the channel pair has a positive significant correlation
  2192. % if Adj(iCh, jCh) < 0
  2193. % clear DataPlot;
  2194. %
  2195. % for iStimGroup = 1:length(SubjGWPLI)
  2196. % DataPlot{iStimGroup} = squeeze(WPLIall(SubjGWPLI{iStimGroup}, :, EEGchInd(iCh), EEGchInd(jCh), iTrialType));
  2197. % Invalid = isnan(DataPlot{iStimGroup}(:, 1));
  2198. % DataPlot{iStimGroup}(Invalid, :) = [];
  2199. % end
  2200. %
  2201. % if isempty(DataPlot{1}) || isempty(DataPlot{2})
  2202. % continue;
  2203. % end
  2204. %
  2205. % iPlot = iPlot + 1;
  2206. %
  2207. % % Plot WPLI data for the channel pair
  2208. % figure;
  2209. % xrange = [0 50]; % Frequency range for x-axis
  2210. % subplot('Position', [0.1 0.1 0.88 0.88]);
  2211. % RateHist_GroupPlot(Fplot, DataPlot, FlickerColor, ParamWPLI); % Custom function for plotting
  2212. %
  2213. % text(range(xrange)/2 - 5, ParamWPLI.Ytick(end), [EEGch{iCh} '-' EEGch{jCh}], 'FontSize', 10); % Add channel pair label
  2214. % set(gca, 'xlim', xrange, 'xtick', [1 4 8 13 30 40 50 60 80 100], 'ylim', [0 0.2], 'ytick', ParamWPLI.Ytick);
  2215. % xlabel('Frequency (Hz)');
  2216. % ylabel('WPLI');
  2217. %
  2218. % % Formatting for grid and font size
  2219. % ax = gca;
  2220. % ax.XGrid = 'on';
  2221. % ax.YGrid = 'off';
  2222. % ax.FontSize = 7;
  2223. %
  2224. % % Set paper size and save the figure
  2225. % LuFontStandard; % Custom function for standard fonts
  2226. % papersizePX = 1.3 * [0 0 5 3.2]; % Custom figure size
  2227. % set(gcf, 'PaperUnits', 'centimeters');
  2228. % set(gcf, 'PaperPosition', papersizePX, 'PaperSize', papersizePX(3:4));
  2229. %
  2230. % % Save the figure in SVG format
  2231. % % saveas(gcf, [SaveTemp 'PositiveAccCorr\' '40LightRandBands_' EEGch{iCh} '-' EEGch{jCh}], 'svg');
  2232. % saveas(gcf, [SaveTemp negWPLIRTChPairFolder 'WPLI_40LR_' EEGch{iCh} '-' EEGch{jCh}], 'png');
  2233. % end
  2234. % end
  2235. % end
  2236. % figure(SpearmanRTFCtopo);
  2237. %% Plot & axises ***fig4f
  2238. FC_BrainEEGLu(ChanPosColin27,Adj,[sum(Adj,1)/2],FCspear)
  2239. nSignifChPairs = sum(AdjP<pSpearEdgeTh & AdjP>0, 'all');
  2240. axis off
  2241. if iBOI==1
  2242. % yt=text(0,0.1,0.1,'Sig-Corr. FC-RT','horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  2243. end
  2244. % text(-0.1,0,0.1,[BandName{iBOI} ' (' BandHzName{iBOI} ')'],'horizontalalignment','center','verticalalignment','top','fontsize',10)
  2245. text(-0.11,0,-0.1,['nChpairs=' num2str(nSignifChPairs)],'horizontalalignment','center','verticalalignment','top','fontsize',10)
  2246. if iBOI==2
  2247. % title(dataName, 'Units', 'normalized', 'Position', [0.5, 0.9, 0]) % MKA 9/23
  2248. % title(dataName, 'Units', 'normalized', 'Position', [0, 0.9, 0])
  2249. end
  2250. %% Display the list of significant channel pairs below the subplot
  2251. % signifChPairNames = getListofSignifChPairs(signifchPairRRTTopo, EEGch);
  2252. %
  2253. % % New subplot for channel pair names (as three columns)
  2254. % subplotLU(2, size(BOI, 2), 2, iBOI); % New row (2nd row) for the names
  2255. %
  2256. % % Divide the list into three columns
  2257. % numNames = length(signifChPairNames);
  2258. % numPerCol = ceil(numNames / 3);
  2259. %
  2260. % % Split the list of significant channel pairs into three columns
  2261. % col1 = signifChPairNames(1:numPerCol);
  2262. % col2 = signifChPairNames(numPerCol+1:min(2*numPerCol, numNames));
  2263. % col3 = signifChPairNames(2*numPerCol+1:end);
  2264. %
  2265. % % Prepare the text to display in columns
  2266. % colText = sprintf('%s\n', col1{:});
  2267. % colText2 = sprintf('%s\n', col2{:});
  2268. % colText3 = sprintf('%s\n', col3{:});
  2269. %
  2270. % % Display the three columns of significant channel pairs
  2271. % text(0.2, 0.5, colText, 'Units', 'normalized', 'HorizontalAlignment', 'left', 'FontSize', 6);
  2272. % text(0.5, 0.5, colText2, 'Units', 'normalized', 'HorizontalAlignment', 'left', 'FontSize', 6);
  2273. % text(0.8, 0.5, colText3, 'Units', 'normalized', 'HorizontalAlignment', 'left', 'FontSize', 6);
  2274. % axis off
  2275. end
  2276. LuFontStandard;
  2277. papersizePX=[0 0 6*size(BOI,2) 6+2];
  2278. set(gcf, 'PaperUnits', 'centimeters');
  2279. set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  2280. SpearmanFCRTName2 = ['SpearmanFCRT_ALPHA' '_' dataName];
  2281. % saveas(gcf,[SaveTemp SpearmanFCRTName2],'pdf');
  2282. saveas(gcf,[SaveTemp SpearmanFCRTName2 '.png'],'png');
  2283. saveas(gcf,[SaveTemp SpearmanFCRTName2 '.svg'],'svg');
  2284. % saveas(gcf,[SaveTemp SpearmanFCRTName2 '.eps'],'epsc');
  2285. close all
  2286. end
  2287. %% Spearman Acc-FC - Topo plot (pre 6/18/24) 40 vs Light (All BANDS)
  2288. % figure;
  2289. %
  2290. % Group1 = 1; % 1= 40Hz flicker group
  2291. % Group2 = 2; % 2 = Light group
  2292. % dataName = [TrialTypeName{iTrialType} FreqFuncNames{iFreqFunc} GroupName{Group1} GroupName{Group2}];
  2293. %
  2294. % IncludedSubj=union(SubjG{Group1},SubjG{Group2}); %40 vs Light
  2295. % for iBOI=1:size(BOI,2)
  2296. % NeedI=find(Fplot>=BOI(1,iBOI)&Fplot<=BOI(2,iBOI));
  2297. % if iFreqFunc==3
  2298. % WPLItemp=squeeze(FreqFunc{iFreqFunc}(WPLIall(IncludedSubj,NeedI,EEGchInd,EEGchInd,iTrialType),[],2));
  2299. % else
  2300. % WPLItemp=squeeze(FreqFunc{iFreqFunc}(WPLIall(IncludedSubj,NeedI,EEGchInd,EEGchInd,iTrialType),2));
  2301. % end
  2302. % Acctemp=Acc(IncludedSubj);
  2303. %
  2304. % for iCh=1:length(ChanEEGLab)-1
  2305. % % for jCh=2:length(ChanEEGLab)
  2306. % WPLIijChPair = squeeze(WPLItemp(:,iCh,iCh+1:end));
  2307. % [rSpear(iCh,iCh+1:end,iBOI,iTrialType),pSpear(iCh,iCh+1:end,iBOI,iTrialType)]=corr(WPLIijChPair,Acctemp,'type','spearman','rows','pairwise');
  2308. % % end
  2309. % end
  2310. % subplotLU(1,size(BOI,2),1,iBOI);
  2311. %
  2312. % Adj=rSpear(:,:,iBOI,iTrialType);
  2313. % AdjP=pSpear(:,:,iBOI,iTrialType);
  2314. % Adj(AdjP>pSpearEdgeTh)=0;
  2315. % FC_BrainEEGLu(ChanPosColin27,Adj,[],FCspear)
  2316. % axis off
  2317. % if iBOI==1
  2318. % yt=text(0,0.1,0.1,'Sig-Corr. FC-Acc','horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  2319. % end
  2320. % text(-0.1,0,0.1,[BandName{iBOI} ' (' BandHzName{iBOI} ')'],'horizontalalignment','center','verticalalignment','top','fontsize',10)
  2321. %
  2322. % if iBOI==4
  2323. % title(dataName, 'Units', 'normalized', 'Position', [0.5, 0.9, 0])
  2324. % end
  2325. % end
  2326. %
  2327. % LuFontStandard;
  2328. % papersizePX=[0 0 6*size(BOI,2) 6+2];
  2329. % set(gcf, 'PaperUnits', 'centimeters');
  2330. % set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  2331. %
  2332. %
  2333. % SpearmanFCAccName1 = ['SpearmanFCAcc' dataName];
  2334. % saveas(gcf,[SaveTemp SpearmanFCAccName1],'pdf');
  2335. % saveas(gcf,[SaveTemp SpearmanFCAccName1],'png');
  2336. % saveas(gcf,[SaveTemp SpearmanFCAccName1 '.eps'],'epsc');
  2337. % close all
  2338. %% Spearman Acc-FC - Topo plot (pre 6/18/24) (40 vs Light) ALPHA
  2339. % BOI=[8 8 10; 13 10 13];
  2340. % BandName={'Alpha','LowAlpha','HighAlpha'};
  2341. % % Gamma-E is 40
  2342. % BandHzName={'8-13Hz','8-10Hz','10-13Hz',};
  2343. % figure;
  2344. %
  2345. % Group1 = 1; % 1= 40Hz flicker group
  2346. % Group2 = 2; % 2 = Light group
  2347. % dataName = [TrialTypeName{iTrialType} FreqFuncNames{iFreqFunc} GroupName{Group1} GroupName{Group2}];
  2348. %
  2349. % IncludedSubj=union(SubjG{Group1},SubjG{Group2}); %40 + Light
  2350. % for iBOI=1:size(BOI,2)
  2351. % NeedI=find(Fplot>=BOI(1,iBOI)&Fplot<=BOI(2,iBOI));
  2352. % if iFreqFunc==3
  2353. % WPLItemp=squeeze(FreqFunc{iFreqFunc}(WPLIall(IncludedSubj,NeedI,EEGchInd,EEGchInd,iTrialType),[],2));
  2354. % else
  2355. % WPLItemp=squeeze(FreqFunc{iFreqFunc}(WPLIall(IncludedSubj,NeedI,EEGchInd,EEGchInd,iTrialType),2));
  2356. % end
  2357. % Acctemp=Acc(IncludedSubj);
  2358. %
  2359. % for iCh=1:length(ChanEEGLab)-1
  2360. % % for jCh=2:length(ChanEEGLab)
  2361. % [rSpear(iCh,iCh+1:end,iBOI,iTrialType),pSpear(iCh,iCh+1:end,iBOI,iTrialType)]=corr(squeeze(WPLItemp(:,iCh,iCh+1:end)),Acctemp,'type','spearman','rows','pairwise');
  2362. % % end
  2363. % end
  2364. % subplotLU(1,size(BOI,2),1,iBOI);
  2365. %
  2366. % Adj=rSpear(:,:,iBOI,iTrialType);
  2367. % AdjP=pSpear(:,:,iBOI,iTrialType);
  2368. % Adj(AdjP>pSpearEdgeTh)=0;
  2369. % FC_BrainEEGLu(ChanPosColin27,Adj,[],FCspear)
  2370. %
  2371. % axis off
  2372. % if iBOI==1
  2373. % yt=text(0,0.1,0.1,'Sig-Corr. FC-Acc','horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  2374. % end
  2375. % text(-0.1,0,0.1,[BandName{iBOI} ' (' BandHzName{iBOI} ')'],'horizontalalignment','center','verticalalignment','top','fontsize',10)
  2376. %
  2377. % if iBOI==2
  2378. % title(dataName, 'Units', 'normalized', 'Position', [0.5, 0.9, 0])
  2379. % end
  2380. % end
  2381. %
  2382. % LuFontStandard;
  2383. % papersizePX=[0 0 6*size(BOI,2) 6+2];
  2384. % set(gcf, 'PaperUnits', 'centimeters');
  2385. % set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  2386. %
  2387. %
  2388. % SpearmanFCAccName1 = ['SpearmanFCAccALPHA_p' num2str(pSpearEdgeTh*100) '_' dataName];
  2389. % % saveas(gcf,[SaveTemp SpearmanFCAccName1],'pdf');
  2390. % saveas(gcf,[SaveTemp SpearmanFCAccName1],'png');
  2391. % % saveas(gcf,[SaveTemp SpearmanFCAccName1 '.eps'],'epsc');
  2392. % close all
  2393. %% Spearman Acc-FC - Topo plot (40 vs Random) All BANDS
  2394. % BOI=[1 4 8 13 30 39.5 43;4 8 13 30 37 41.5 100];
  2395. % BandName={'Delta','Theta','Alpha','Beta','Gamma-1','Gamma-E','Gamma-2'};
  2396. % % Gamma-E is 40
  2397. % BandHzName={'1-4 Hz','4-8 Hz','8-13Hz','13-30 Hz','30-37 Hz','39-41 Hz','43-100 Hz'};
  2398. % figure;
  2399. %
  2400. % Group1 = 1; % 1= 40Hz flicker group
  2401. % Group2 = 3; % 3 = Random group
  2402. % dataName = [TrialTypeName{iTrialType} FreqFuncNames{iFreqFunc} GroupName{Group1} GroupName{Group2}];
  2403. %
  2404. % IncludedSubj=union(SubjG{1},SubjG{3}); % 40 vs Random
  2405. % for iBOI=1:size(BOI,2)
  2406. % NeedI=find(Fplot>=BOI(1,iBOI)&Fplot<=BOI(2,iBOI));
  2407. % if iFreqFunc==3
  2408. % WPLItemp=squeeze(FreqFunc{iFreqFunc}(WPLIall(IncludedSubj,NeedI,EEGchInd,EEGchInd,iTrialType),[],2));
  2409. %
  2410. % else
  2411. % WPLItemp=squeeze(FreqFunc{iFreqFunc}(WPLIall(IncludedSubj,NeedI,EEGchInd,EEGchInd,iTrialType),2));
  2412. % end
  2413. % Acctemp=Acc(IncludedSubj);
  2414. %
  2415. % for iCh=1:length(ChanEEGLab)-1
  2416. % % for jCh=2:length(ChanEEGLab)
  2417. % [rSpear(iCh,iCh+1:end,iBOI,iTrialType),pSpear(iCh,iCh+1:end,iBOI,iTrialType)]=corr(squeeze(WPLItemp(:,iCh,iCh+1:end)),Acctemp,'type','spearman','rows','pairwise');
  2418. % % end
  2419. % end
  2420. % subplotLU(1,size(BOI,2),1,iBOI);
  2421. %
  2422. % Adj=rSpear(:,:,iBOI,iTrialType);
  2423. % AdjP=pSpear(:,:,iBOI,iTrialType);
  2424. % Adj(AdjP>pSpearEdgeTh)=0;
  2425. % FC_BrainEEGLu(ChanPosColin27,Adj,[],FCspear)
  2426. % axis off
  2427. % if iBOI==1
  2428. % yt=text(0,0.1,0.1,'Sig-Corr. FC-Acc','horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  2429. % end
  2430. % text(-0.1,0,0.1,[BandName{iBOI} ' (' BandHzName{iBOI} ')'],'horizontalalignment','center','verticalalignment','top','fontsize',10)
  2431. % if iBOI==4
  2432. % title(dataName, 'Units', 'normalized', 'Position', [0.5, 0.9, 0])
  2433. % end
  2434. % end
  2435. %
  2436. % LuFontStandard;
  2437. % papersizePX=[0 0 6*size(BOI,2) 6+2];
  2438. % set(gcf, 'PaperUnits', 'centimeters');
  2439. % set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  2440. %
  2441. % SpearmanFCAccfigName = ['SpearmanFCAcc' dataName];
  2442. % saveas(gcf,[SaveTemp SpearmanFCAccfigName],'pdf');
  2443. % saveas(gcf,[SaveTemp SpearmanFCAccfigName],'png');
  2444. % saveas(gcf,[SaveTemp SpearmanFCAccfigName 'eps'],'epsc');
  2445. % close all
  2446. %% Spearman Acc-FC - Topo plot (40+RANDOM) ALPHA
  2447. % BOI=[8 8 10; 13 10 13];
  2448. % BandName={'Alpha','LowAlpha','HighAlpha'};
  2449. % % Gamma-E is 40
  2450. % BandHzName={'8-13Hz','8-10Hz','10-13Hz',};
  2451. % figure;
  2452. %
  2453. % Group1 = 1; % 1= 40Hz flicker group
  2454. % Group2 = 3; % 2 = Random group
  2455. % dataName = [TrialTypeName{iTrialType} FreqFuncNames{iFreqFunc} GroupName{Group1} GroupName{Group2}];
  2456. %
  2457. % IncludedSubj=union(SubjG{Group1},SubjG{Group2}); %40 vs Random
  2458. % for iBOI=1:size(BOI,2)
  2459. % NeedI=find(Fplot>=BOI(1,iBOI)&Fplot<=BOI(2,iBOI));
  2460. % if iFreqFunc==3
  2461. % WPLItemp=squeeze(FreqFunc{iFreqFunc}(WPLIall(IncludedSubj,NeedI,EEGchInd,EEGchInd,iTrialType),[],2));
  2462. %
  2463. % else
  2464. % WPLItemp=squeeze(FreqFunc{iFreqFunc}(WPLIall(IncludedSubj,NeedI,EEGchInd,EEGchInd,iTrialType),2));
  2465. % end
  2466. % Acctemp=Acc(IncludedSubj);
  2467. %
  2468. % for iCh=1:length(ChanEEGLab)-1
  2469. % % for jCh=2:length(ChanEEGLab)
  2470. % [rSpear(iCh,iCh+1:end,iBOI,iTrialType),pSpear(iCh,iCh+1:end,iBOI,iTrialType)]=corr(squeeze(WPLItemp(:,iCh,iCh+1:end)),Acctemp,'type','spearman','rows','pairwise');
  2471. % % end
  2472. % end
  2473. % subplotLU(1,size(BOI,2),1,iBOI);
  2474. %
  2475. % Adj=rSpear(:,:,iBOI,iTrialType);
  2476. % AdjP=pSpear(:,:,iBOI,iTrialType);
  2477. % Adj(AdjP>pSpearEdgeTh)=0;
  2478. % FC_BrainEEGLu(ChanPosColin27,Adj,[],FCspear)
  2479. % axis off
  2480. % if iBOI==1
  2481. % yt=text(0,0.1,0.1,'Sig-Corr. FC-Acc','horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  2482. % end
  2483. % text(-0.1,0,0.1,[BandName{iBOI} ' (' BandHzName{iBOI} ')'],'horizontalalignment','center','verticalalignment','top','fontsize',10)
  2484. %
  2485. % if iBOI==2
  2486. % title(dataName, 'Units', 'normalized', 'Position', [0.5, 0.9, 0])
  2487. % end
  2488. % end
  2489. %
  2490. % LuFontStandard;
  2491. % papersizePX=[0 0 6*size(BOI,2) 6+2];
  2492. % set(gcf, 'PaperUnits', 'centimeters');
  2493. % set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  2494. %
  2495. %
  2496. % SpearmanFCAccName1 = ['SpearmanFCAccALPHA' dataName];
  2497. % % saveas(gcf,[SaveTemp SpearmanFCAccName1],'pdf');
  2498. % saveas(gcf,[SaveTemp SpearmanFCAccName1],'png');
  2499. % % saveas(gcf,[SaveTemp SpearmanFCAccName1 '.eps'],'epsc');
  2500. % close all
  2501. %% Spearman RT-FC - Topo plot (40+Light) ALL BANDS
  2502. % BOI=[1 4 8 13 30 39.5 43;4 8 13 30 37 41.5 100];
  2503. % BandName={'Delta','Theta','Alpha','Beta','Gamma-1','Gamma-E','Gamma-2'};
  2504. % % Gamma-E is 40
  2505. % BandHzName={'1-4 Hz','4-8 Hz','8-13Hz','13-30 Hz','30-37 Hz','39-41 Hz','43-100 Hz'};
  2506. %
  2507. % figure;
  2508. %
  2509. % Group1 = 1; % 1= 40Hz flicker group
  2510. % Group2 = 2; % 2 = Light group
  2511. % dataName = [TrialTypeName{iTrialType} FreqFuncNames{iFreqFunc} GroupName{Group1} GroupName{Group2}];
  2512. %
  2513. % IncludedSubj=union(SubjG{Group1},SubjG{Group2}); %40 vs Light
  2514. % for iBOI=1:size(BOI,2)
  2515. % NeedI=find(Fplot>=BOI(1,iBOI)&Fplot<=BOI(2,iBOI));
  2516. % if iFreqFunc==3
  2517. % WPLItemp=squeeze(FreqFunc{iFreqFunc}(WPLIall(IncludedSubj,NeedI,EEGchInd,EEGchInd,iTrialType),[],2));
  2518. %
  2519. % else
  2520. % WPLItemp=squeeze(FreqFunc{iFreqFunc}(WPLIall(IncludedSubj,NeedI,EEGchInd,EEGchInd,iTrialType),2));
  2521. % end
  2522. % RTtemp=SubjsAvgRT(IncludedSubj);
  2523. % %Acctemp=Acc(IncludedSubj);
  2524. %
  2525. % for iCh=1:length(ChanEEGLab)-1
  2526. % % for jCh=2:length(ChanEEGLab)
  2527. % [rSpear(iCh,iCh+1:end,iBOI,iTrialType),pSpear(iCh,iCh+1:end,iBOI,iTrialType)]=corr(squeeze(WPLItemp(:,iCh,iCh+1:end)),RTtemp,'type','spearman','rows','pairwise');
  2528. % % end
  2529. % end
  2530. % subplotLU(1,size(BOI,2),1,iBOI);
  2531. %
  2532. % Adj=rSpear(:,:,iBOI,iTrialType);
  2533. % AdjP=pSpear(:,:,iBOI,iTrialType);
  2534. % Adj(AdjP>pSpearEdgeTh)=0;
  2535. % FC_BrainEEGLu(ChanPosColin27,Adj,[],FCspear)
  2536. % axis off
  2537. % if iBOI==1
  2538. % yt=text(0,0.1,0.1,'Sig-Corr. FC-RT','horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  2539. % end
  2540. % text(-0.1,0,0.1,[BandName{iBOI} ' (' BandHzName{iBOI} ')'],'horizontalalignment','center','verticalalignment','top','fontsize',10)
  2541. %
  2542. % if iBOI==4
  2543. % title(dataName, 'Units', 'normalized', 'Position', [0.5, 0.9, 0])
  2544. % end
  2545. % end
  2546. %
  2547. % LuFontStandard;
  2548. % papersizePX=[0 0 6*size(BOI,2) 6+2];
  2549. % set(gcf, 'PaperUnits', 'centimeters');
  2550. % set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  2551. %
  2552. %
  2553. % SpearmanFCRTName1 = ['SpearmanFCRT' dataName];
  2554. % saveas(gcf,[SaveTemp SpearmanFCRTName1],'pdf');
  2555. % saveas(gcf,[SaveTemp SpearmanFCRTName1],'png');
  2556. % saveas(gcf,[SaveTemp SpearmanFCRTName1 '.eps'],'epsc');
  2557. % close all
  2558. %% Spearman RT-FC - Topo plot (40 vs Random) ALL BANDS
  2559. % BOI=[1 4 8 13 30 39.5 43;4 8 13 30 37 41.5 100];
  2560. % BandName={'Delta','Theta','Alpha','Beta','Gamma-1','Gamma-E','Gamma-2'};
  2561. % % Gamma-E is 40
  2562. % BandHzName={'1-4 Hz','4-8 Hz','8-13Hz','13-30 Hz','30-37 Hz','39-41 Hz','43-100 Hz'};
  2563. % figure;
  2564. %
  2565. % Group1 = 1; % 1= 40Hz flicker group
  2566. % Group2 = 3; % 3 = Random group
  2567. % dataName = [TrialTypeName{iTrialType} FreqFuncNames{iFreqFunc} GroupName{Group1} GroupName{Group2}];
  2568. %
  2569. % IncludedSubj=union(SubjG{Group1},SubjG{Group2}); %40 vs Light
  2570. % for iBOI=1:size(BOI,2)
  2571. % NeedI=find(Fplot>=BOI(1,iBOI)&Fplot<=BOI(2,iBOI));
  2572. % if iFreqFunc==3
  2573. % WPLItemp=squeeze(FreqFunc{iFreqFunc}(WPLIall(IncludedSubj,NeedI,EEGchInd,EEGchInd,iTrialType),[],2));
  2574. %
  2575. % else
  2576. % WPLItemp=squeeze(FreqFunc{iFreqFunc}(WPLIall(IncludedSubj,NeedI,EEGchInd,EEGchInd,iTrialType),2));
  2577. % end
  2578. % RTtemp=SubjsAvgRT(IncludedSubj);
  2579. % %Acctemp=Acc(IncludedSubj);
  2580. %
  2581. % for iCh=1:length(ChanEEGLab)-1
  2582. % % for jCh=2:length(ChanEEGLab)
  2583. % [rSpear(iCh,iCh+1:end,iBOI,iTrialType),pSpear(iCh,iCh+1:end,iBOI,iTrialType)]=corr(squeeze(WPLItemp(:,iCh,iCh+1:end)),RTtemp,'type','spearman','rows','pairwise');
  2584. % % end
  2585. % end
  2586. % subplotLU(1,size(BOI,2),1,iBOI);
  2587. %
  2588. % Adj=rSpear(:,:,iBOI,iTrialType);
  2589. % AdjP=pSpear(:,:,iBOI,iTrialType);
  2590. % Adj(AdjP>pSpearEdgeTh)=0;
  2591. % FC_BrainEEGLu(ChanPosColin27,Adj,[],FCspear)
  2592. % axis off
  2593. % if iBOI==1
  2594. % yt=text(0,0.1,0.1,'Sig-Corr. FC-RT','horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  2595. % end
  2596. % text(-0.1,0,0.1,[BandName{iBOI} ' (' BandHzName{iBOI} ')'],'horizontalalignment','center','verticalalignment','top','fontsize',10)
  2597. %
  2598. % if iBOI==4
  2599. % title(dataName, 'Units', 'normalized', 'Position', [0.5, 0.9, 0])
  2600. % end
  2601. % end
  2602. %
  2603. % LuFontStandard;
  2604. % papersizePX=[0 0 6*size(BOI,2) 6+2];
  2605. % set(gcf, 'PaperUnits', 'centimeters');
  2606. % set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  2607. %
  2608. %
  2609. % SpearmanFCRTName2 = ['SpearmanFCRT' dataName];
  2610. % saveas(gcf,[SaveTemp SpearmanFCRTName2],'pdf');
  2611. % saveas(gcf,[SaveTemp SpearmanFCRTName2],'png');
  2612. % saveas(gcf,[SaveTemp SpearmanFCRTName2 '.eps'],'epsc');
  2613. % close all
  2614. %% Spearman RT-FC - Topo plot (40+Light) ALPHA ONLY
  2615. % BOI=[8 8 10; 13 10 13];
  2616. % BandName={'Alpha','LowAlpha','HighAlpha'};
  2617. % % Gamma-E is 40
  2618. % BandHzName={'8-13Hz','8-10Hz','10-13Hz',};
  2619. %
  2620. % figure;
  2621. %
  2622. % Group1 = 1; % 1= 40Hz flicker group
  2623. % Group2 = 2; % 2 = Light group
  2624. % dataName = [TrialTypeName{iTrialType} FreqFuncNames{iFreqFunc} GroupName{Group1} GroupName{Group2}];
  2625. %
  2626. % IncludedSubj=union(SubjG{Group1},SubjG{Group2}); %40 vs Light
  2627. % for iBOI=1:size(BOI,2)
  2628. % NeedI=find(Fplot>=BOI(1,iBOI)&Fplot<=BOI(2,iBOI));
  2629. % if iFreqFunc==3
  2630. % WPLItemp=squeeze(FreqFunc{iFreqFunc}(WPLIall(IncludedSubj,NeedI,EEGchInd,EEGchInd,iTrialType),[],2));
  2631. %
  2632. % else
  2633. % WPLItemp=squeeze(FreqFunc{iFreqFunc}(WPLIall(IncludedSubj,NeedI,EEGchInd,EEGchInd,iTrialType),2));
  2634. % end
  2635. % RTtemp=SubjsAvgRT(IncludedSubj);
  2636. % %Acctemp=Acc(IncludedSubj);
  2637. %
  2638. % for iCh=1:length(ChanEEGLab)-1
  2639. % % for jCh=2:length(ChanEEGLab)
  2640. % [rSpear(iCh,iCh+1:end,iBOI,iTrialType),pSpear(iCh,iCh+1:end,iBOI,iTrialType)]=corr(squeeze(WPLItemp(:,iCh,iCh+1:end)),RTtemp,'type','spearman','rows','pairwise');
  2641. % % end
  2642. % end
  2643. % subplotLU(1,size(BOI,2),1,iBOI);
  2644. %
  2645. % Adj=rSpear(:,:,iBOI,iTrialType);
  2646. % AdjP=pSpear(:,:,iBOI,iTrialType);
  2647. % Adj(AdjP>pSpearEdgeTh)=0;
  2648. % FC_BrainEEGLu(ChanPosColin27,Adj,[],FCspear)
  2649. % axis off
  2650. % if iBOI==1
  2651. % yt=text(0,0.1,0.1,'Sig-Corr. FC-RT','horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  2652. % end
  2653. % text(-0.1,0,0.1,[BandName{iBOI} ' (' BandHzName{iBOI} ')'],'horizontalalignment','center','verticalalignment','top','fontsize',10)
  2654. %
  2655. % if iBOI==4
  2656. % title(dataName, 'Units', 'normalized', 'Position', [0.5, 0.9, 0])
  2657. % end
  2658. % end
  2659. %
  2660. % LuFontStandard;
  2661. % papersizePX=[0 0 6*size(BOI,2) 6+2];
  2662. % set(gcf, 'PaperUnits', 'centimeters');
  2663. % set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  2664. %
  2665. %
  2666. % SpearmanFCRTName1 = ['SpearmanFCRTALPHA' dataName];
  2667. % % saveas(gcf,[SaveTemp SpearmanFCRTName1],'pdf');
  2668. % saveas(gcf,[SaveTemp SpearmanFCRTName1],'png');
  2669. % saveas(gcf,[SaveTemp SpearmanFCRTName1 '.eps'],'epsc');
  2670. % close all
  2671. %% Spearman RT-FC - Topo plot (40 vs Random) ALPHA ONLY
  2672. % BOI=[8 8 10; 13 10 13];
  2673. % BandName={'Alpha','LowAlpha','HighAlpha'};
  2674. % % Gamma-E is 40
  2675. % BandHzName={'8-13Hz','8-10Hz','10-13Hz',};
  2676. %
  2677. % figure;
  2678. %
  2679. % Group1 = 1; % 1= 40Hz flicker group
  2680. % Group2 = 3; % 3 = Random group
  2681. % dataName = [TrialTypeName{iTrialType} FreqFuncNames{iFreqFunc} GroupName{Group1} GroupName{Group2}];
  2682. %
  2683. % IncludedSubj=union(SubjG{Group1},SubjG{Group2}); %40 vs Random
  2684. % for iBOI=1:size(BOI,2)
  2685. % NeedI=find(Fplot>=BOI(1,iBOI)&Fplot<=BOI(2,iBOI));
  2686. % if iFreqFunc==3
  2687. % WPLItemp=squeeze(FreqFunc{iFreqFunc}(WPLIall(IncludedSubj,NeedI,EEGchInd,EEGchInd,iTrialType),[],2));
  2688. %
  2689. % else
  2690. % WPLItemp=squeeze(FreqFunc{iFreqFunc}(WPLIall(IncludedSubj,NeedI,EEGchInd,EEGchInd,iTrialType),2));
  2691. % end
  2692. % RTtemp=SubjsAvgRT(IncludedSubj);
  2693. % %Acctemp=Acc(IncludedSubj);
  2694. %
  2695. % for iCh=1:length(ChanEEGLab)-1
  2696. % % for jCh=2:length(ChanEEGLab)
  2697. % [rSpear(iCh,iCh+1:end,iBOI,iTrialType),pSpear(iCh,iCh+1:end,iBOI,iTrialType)]=corr(squeeze(WPLItemp(:,iCh,iCh+1:end)),RTtemp,'type','spearman','rows','pairwise');
  2698. % % end
  2699. % end
  2700. % subplotLU(1,size(BOI,2),1,iBOI);
  2701. %
  2702. % Adj=rSpear(:,:,iBOI,iTrialType);
  2703. % AdjP=pSpear(:,:,iBOI,iTrialType);
  2704. % Adj(AdjP>pSpearEdgeTh)=0;
  2705. % FC_BrainEEGLu(ChanPosColin27,Adj,[],FCspear)
  2706. % axis off
  2707. % if iBOI==1
  2708. % yt=text(0,0.1,0.1,'Sig-Corr. FC-RT','horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  2709. % end
  2710. % text(-0.1,0,0.1,[BandName{iBOI} ' (' BandHzName{iBOI} ')'],'horizontalalignment','center','verticalalignment','top','fontsize',10)
  2711. %
  2712. % if iBOI==2
  2713. % title(dataName, 'Units', 'normalized', 'Position', [0.5, 0.9, 0])
  2714. % end
  2715. % end
  2716. %
  2717. % LuFontStandard;
  2718. % papersizePX=[0 0 6*size(BOI,2) 6+2];
  2719. % set(gcf, 'PaperUnits', 'centimeters');
  2720. % set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  2721. %
  2722. %
  2723. % SpearmanFCRTName2 = ['SpearmanFCRTALPHA' dataName];
  2724. % % saveas(gcf,[SaveTemp SpearmanFCRTName2],'pdf');
  2725. % saveas(gcf,[SaveTemp SpearmanFCRTName2],'png');
  2726. % % saveas(gcf,[SaveTemp SpearmanFCRTName2 '.eps'],'epsc');
  2727. % close all
  2728. end
  2729. end
  2730. close all
  2731. %% PSD fig2b
  2732. SubSavePSD=[SavePath 'PSD\'];
  2733. PowerLim=[-7 -2];
  2734. PowerLab={'-7' '-2'};
  2735. TLim=[-5 5];
  2736. TLab={'-5' '5'};
  2737. PLim=[-4 4];
  2738. PLab={'10e-4' '10e-0'};
  2739. for iFreqFunc=1:length(FreqFuncNames)
  2740. clear MapGroup diffTmap diffMap Tdata;
  2741. clear diffTPmap diffRPmap diffTmap
  2742. for iTrialType=1:length(TrialType)
  2743. % iCom=1;
  2744. SaveTemp=[SubSavePSD TrialTypeName{iTrialType} '\'];
  2745. SaveTemp=[SaveTemp FreqFuncNames{iFreqFunc} '\' ];
  2746. mkdir(SaveTemp)
  2747. % Band of interest
  2748. BOI=[5;15];
  2749. BOI=[1 4 8 13 30 39.5 43;4 8 13 30 37 41.5 100];
  2750. BandName={'Delta','Theta','Alpha','Beta','Gamma-1','Gamma-E','Gamma-2'};
  2751. BandHzName={'1-4 Hz','4-8 Hz','8-13Hz','13-30 Hz','30-37 Hz','39-41 Hz','43-100 Hz'};
  2752. %
  2753. figure;
  2754. IncludedSubj=union(SubjG{1},SubjG{2});
  2755. for iBOI=1:size(BOI,2)
  2756. NeedI=find(Fplot>=BOI(1,iBOI)&Fplot<=BOI(2,iBOI));
  2757. % for iCh=1:length(ChanEEGLab)
  2758. if iFreqFunc==3
  2759. PSDtemp=squeeze(FreqFunc{iFreqFunc}(LogPSD(IncludedSubj,NeedI,EEGchInd,iTrialType),[],2));
  2760. else
  2761. PSDtemp=squeeze(FreqFunc{iFreqFunc}(LogPSD(IncludedSubj,NeedI,EEGchInd,iTrialType),2));
  2762. end
  2763. Acctemp=Acc(IncludedSubj); % Behavior accuracy
  2764. % Power and acuracy correlation
  2765. [rSpear(:,iBOI,iTrialType),pSpear(:,iBOI,iTrialType)]=corr(PSDtemp,Acctemp,'type','spearman','rows','pairwise');
  2766. % end
  2767. %% Visualize the data with topoplot in eeglab
  2768. subplotLU(2,size(BOI,2),1,iBOI);
  2769. topoplot(rSpear(:,iBOI,iTrialType), ChanEEGLab,'colormap',colorMapPN,'maplimits',[-1 1]);
  2770. if iBOI==1
  2771. yt=text(-0.55,0,'EEG-Behavior R','horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  2772. end
  2773. if iBOI==size(BOI,2)
  2774. b=colorbar('southoutside');set(b,'position',[0.52 0.93 0.2 0.01],'xtick',[-1 1],'xticklabel',{'-1' '1'},'xlim',[-1 1]);
  2775. xlabel(b,'PSD-Acc Correlation','verticalalignment','top')
  2776. end
  2777. subplotLU(2,size(BOI,2),2,iBOI);
  2778. topoplot(log10(pSpear(:,iBOI,iTrialType)), ChanEEGLab,'colormap',colorMapPN,'maplimits',[-4 4]); % p-value
  2779. xlabel(BandName{iBOI});
  2780. text(0,-0.55,[BandName{iBOI} ' (' BandHzName{iBOI} ')'],'horizontalalignment','center','verticalalignment','top','fontsize',10)
  2781. if iBOI==1
  2782. a=ylabel('40Hz-BothControls');
  2783. % set(a,'position',[0.01 0.5 0.03 0.4],'verticalalignment','middle')
  2784. set(a,'verticalalignment','middle')
  2785. yt=text(-0.55,0,'P EEG-Behavior','horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  2786. end
  2787. if iBOI==size(BOI,2)
  2788. % c=colorbar('southoutside');set(c,'position',[0.52 0.46 0.2 0.03],'xtick',[-4 0],'xticklabel',{'10e-4' '10e-0'},'xlim',[-4 0]);
  2789. % xlabel(c,'P values','verticalalignment','top')
  2790. c=colorbar('southoutside');
  2791. set(gca,'xtick',[],'ytick',[])
  2792. set(c,'position',[0.52 0.51 0.2 0.01],'Limits',[PLim(1) 0],'ticks',[PLim(1) 0],'ticklabels',PLab);
  2793. xlabel(c,'P values','verticalalignment','top')
  2794. end
  2795. end
  2796. LuFontStandard;
  2797. papersizePX=[0 0 6*size(BOI,2) 6*2+3];
  2798. set(gcf, 'PaperUnits', 'centimeters');
  2799. set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  2800. saveas(gcf,[SaveTemp 'SpearmanPSDAcc'],'pdf');
  2801. saveas(gcf,[SaveTemp 'SpearmanPSDAcc'],'png');
  2802. saveas(gcf,[SaveTemp 'SpearmanPSDAcc.eps'],'epsc');
  2803. close all
  2804. %% Scatter plot with behavior - PSDAcc fig
  2805. IncludedSubj=[SubjG{1}(:);SubjG{2}(:)];
  2806. FlickerID=[zeros(size(SubjG{1}(:)))+1;zeros(size(SubjG{2}(:)))+2];
  2807. Param.Corr='Spearman'; %%%Type of correlation, see Matlab function corr for more details
  2808. Param.Pth=0.05; %%%threshold of Pvalue
  2809. Param.ColorMap=colorMapPN; %%%Color map for correlation link
  2810. Param.NodeColor=repmat([0.8 0.8 0.8],6,1); %%%Color of Nodes for correlation link plot
  2811. Param.Clim=[-0.6 0.6]; %%%Color Limit for Correlation
  2812. Param.Title='Pool All Sample'; %%Any title for label the figure
  2813. Param.MarkerSize=8; %%%MarkerSize of scatter
  2814. Param.SubjIDColor=FlickerColor;
  2815. Param.SubjID=FlickerID;
  2816. for iBOI=1:size(BOI,2)
  2817. NeedI=find(Fplot>=BOI(1,iBOI)&Fplot<=BOI(2,iBOI));
  2818. figure;
  2819. % for iCh=1:length(ChanEEGLab)
  2820. if iFreqFunc==3
  2821. PSDtemp=squeeze(FreqFunc{iFreqFunc}(LogPSD(IncludedSubj,NeedI,EEGchInd,iTrialType),[],2));
  2822. else
  2823. PSDtemp=squeeze(FreqFunc{iFreqFunc}(LogPSD(IncludedSubj,NeedI,EEGchInd,iTrialType),2));
  2824. end
  2825. Acctemp=Acc(IncludedSubj);
  2826. tempName=EEGch;
  2827. tempName{end+1}='Acc';
  2828. multiCorr2GroupSubplot(6,6,[PSDtemp Acctemp],size(PSDtemp,2)+1,tempName,Param)
  2829. LuFontStandard;
  2830. papersizePX=[0 0 6*6 6*6];
  2831. set(gcf, 'PaperUnits', 'centimeters');
  2832. set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  2833. saveas(gcf,[SaveTemp Param.Corr BandName{iBOI} BandHzName{iBOI} 'PSDAcc'],'pdf');
  2834. saveas(gcf,[SaveTemp Param.Corr BandName{iBOI} BandHzName{iBOI} 'PSDAcc'],'png');
  2835. saveas(gcf,[SaveTemp Param.Corr BandName{iBOI} BandHzName{iBOI} 'PSDAcc.eps'],'epsc');
  2836. close all
  2837. end
  2838. end
  2839. %
  2840. % close all
  2841. %
  2842. end
  2843. FreqFunc{1}=@nanmean;
  2844. FreqFunc{2}=@nanmedian;
  2845. FreqFunc{3}=@nanmax;
  2846. FreqFuncNames={'mean','median','peak'};
  2847. diffTPmap=zeros(length(ChanEEGLab),length(ChanEEGLab),size(BOI,2),length(TrialType));
  2848. diffRPmap=diffTPmap;
  2849. diffTmap=diffTPmap;
  2850. rSpear=diffTPmap;
  2851. pSpear=diffTPmap;
  2852. % % TNodeTh=10;
  2853. COHEdgeTh=0.2;
  2854. COHNodeTh=0.01;
  2855. TEdgeTh=3;
  2856. TNodeTh=0.001;
  2857. pTEdgeTh=0.05;
  2858. pSpearEdgeTh=0.05;
  2859. rSpearEdgeTh=0.1;
  2860. rSpearNodeTh=0.001;
  2861. ScaleCOH=0.5;
  2862. ScaleT=1;
  2863. ScaleSpear=0.2;
  2864. close all
  2865. %% PSD all three groups together plot- Ranktest fig2b
  2866. SubSavePSD=[SavePath 'PSD\'];
  2867. % SubSaveCOH=[SavePath 'COH\'];
  2868. ParamPSD.Ytick = [-8.1:3:-2.1];
  2869. ParamPSD.SigPlot='Ranktest';
  2870. saveDate = datestr(datetime, 'yy-mm-dd_HHMMSSFFF');
  2871. SubSavePSD=[SubSavePSD 'Ranktest_' saveDate '\'];
  2872. mkdir(SubSavePSD)
  2873. for iTrialType=3 %1:length(TrialType)
  2874. SaveTemp=[SubSavePSD TrialTypeName{iTrialType} '\'];
  2875. mkdir(SaveTemp)
  2876. SubSaveFig=[SaveTemp 'Chan\'];
  2877. mkdir(SubSaveFig)
  2878. FBand=[1 55];
  2879. %% HzAllChPSD40Random
  2880. figure;
  2881. for iCh=1:length(EEGchInd)
  2882. clear DataPlot
  2883. for iStimGroup=1:length(SubjG)
  2884. DataPlot{iStimGroup}= squeeze(LogPSD(SubjG{iStimGroup},:,EEGchInd(iCh),iTrialType));
  2885. Invalid=isnan(DataPlot{iStimGroup}(:,1));
  2886. DataPlot{iStimGroup}(Invalid,:)=[];
  2887. end
  2888. if isempty(DataPlot{1})||isempty(DataPlot{2})
  2889. continue;
  2890. end
  2891. % subplotLUpage(6,6,iCh);
  2892. ParamPSD.PathSave=[SaveTemp '40LightRandCh' EEGch{iCh}];
  2893. % figure;
  2894. % subplot('Position',[0.1 0.1 0.88 0.88])
  2895. [~,FlickComStatis{iTrialType,iCh}]=RateHist_GroupPlot(Fplot,DataPlot,FlickerColor,ParamPSD); %plot and stats,
  2896. text(27.5,-2,EEGch{iCh}); % text(27.5,0.3,EEGch{iCh}); %pre 6/14/24
  2897. set(gca,'xlim',FBand,'xtick',[1 4 8 13 30 37 39 41 43 55],'ylim',[-7 -2],'ytick',[-7 -6 -5 -4 -3 -2]);
  2898. xlabel('Frequency Hz')
  2899. ylabel('Normalized Power (Log)')
  2900. ax=gca;
  2901. ax.XGrid = 'on';
  2902. ax.YGrid = 'off';
  2903. LuFontStandard
  2904. papersizePX=[0 0 12 12];
  2905. set(gcf, 'PaperUnits', 'centimeters');
  2906. set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  2907. saveas(gcf,[SubSaveFig 'ThreeGroup_' EEGch{iCh}],'tiff');
  2908. saveas(gcf,[SubSaveFig 'ThreeGroup_' EEGch{iCh}],'png');
  2909. saveas(gcf,[SubSaveFig 'ThreeGroup_' EEGch{iCh}],'svg');
  2910. saveas(gcf,[SubSaveFig 'ThreeGroup_' EEGch{iCh} '.eps'],'epsc');
  2911. close all
  2912. end
  2913. %
  2914. papersizePX=[0 0 6*6 6*6];
  2915. set(gcf, 'PaperUnits', 'centimeters');
  2916. set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  2917. saveas(gcf,[SaveTemp num2str(FBand(1)) '-' num2str(FBand(2)) 'HzAllChPSD40LightRand' TrialTypeName{iTrialType}],'pdf');
  2918. saveas(gcf,[SaveTemp num2str(FBand(1)) '-' num2str(FBand(2)) 'HzAllChPSD40LightRand' TrialTypeName{iTrialType}],'png');
  2919. saveas(gcf,[SaveTemp num2str(FBand(1)) '-' num2str(FBand(2)) 'HzAllChPSD40LightRand' TrialTypeName{iTrialType}],'svg');
  2920. saveas(gcf,[SaveTemp num2str(FBand(1)) '-' num2str(FBand(2)) 'HzAllChPSD40LightRand' TrialTypeName{iTrialType} '.eps'],'epsc');
  2921. end
  2922. close all
  2923. ParamPSD.Ytick=[-8:4:0]; % reset this param to original value to not mess up other functions: ParamPSD.Ytick=[-8:4:0];
  2924. %% Set-up FOOOF save folder
  2925. fooof_starttime = datestr(datetime, 'yy-mm-dd_HHMMSSFFF');
  2926. fooof_save_folder=['FOOOF Results' '\' fooof_starttime '\'];
  2927. mkdir(fooof_save_folder)
  2928. %% Plot single FOOOF PSD
  2929. % set up fooof input variables:
  2930. allFreqs = Fplot; % = row vector of frequency values
  2931. iStimGroup = 3; % 1=40; 2=Light, 3=Random;
  2932. trialtype = 3; % 3=HIT&MISS
  2933. singleFOOOFEEGCh = 1;
  2934. subjectNumInStimGroup = 1;
  2935. input_power_spectrum = PSDall(SubjG{iStimGroup}(subjectNumInStimGroup),:,singleFOOOFEEGCh,trialtype); % size(PSDall) = 67 197 40 3
  2936. % PSDall(participants, frequencies, channels, trialtypes)
  2937. f_range = [2 55]; % f_range = fitting range - !!different upper ranges yield different results
  2938. % set default settings values - % settings = fooof model settings, in a struct
  2939. settings = struct(...
  2940. 'peak_width_limits', [0.5, 12], ...
  2941. 'max_n_peaks', Inf, ...
  2942. 'min_peak_height', 0.0, ...
  2943. 'peak_threshold', 2.0, ...
  2944. 'aperiodic_mode', 'fixed', ...
  2945. 'verbose', true);
  2946. return_model = 1; % return_model = boolean of whether to return the FOOOF model fit, optional
  2947. % run fooof:
  2948. firstsubjectfirstCh = fooof(Fplot, input_power_spectrum, f_range, settings, return_model);
  2949. % Plot FOOOF results:
  2950. % Extract data from the results structure
  2951. fooof_freqs = firstsubjectfirstCh.freqs; % Frequency values
  2952. power_spectrum = firstsubjectfirstCh.power_spectrum; % Original power spectrum
  2953. % fooofed_spectrum = firstsubjectfirstCh.fooofed_spectrum; % Full FOOOF fit
  2954. ap_fit = firstsubjectfirstCh.ap_fit; % Aperiodic fit
  2955. difference_spectrum = power_spectrum - ap_fit;
  2956. % fooofed_difference_spectrum = fooofed_spectrum - ap_fit;
  2957. % Plot settings
  2958. figure;
  2959. hold on;
  2960. % Plot the original power spectrum
  2961. plot(fooof_freqs, power_spectrum, 'k', 'LineWidth', 1.5, 'DisplayName', 'Power Spectrum');
  2962. % Plot the FOOOFed spectrum (from fooof results)
  2963. % plot(fooof_freqs, fooofed_spectrum, 'r', 'LineWidth', 1.5, 'DisplayName', 'FOOOFed Spectrum');
  2964. % Plot the aperiodic fit
  2965. plot(fooof_freqs, ap_fit, 'b--', 'LineWidth', 1.5, 'DisplayName', 'Aperiodic Fit');
  2966. % Plot the difference spectrum (subtracted in MATLAB)
  2967. plot(fooof_freqs, difference_spectrum, 'g--', 'LineWidth', 1.5, 'DisplayName', 'Adjusted Spectrum (PS-ApFit)');
  2968. % % Plot the FOOOFed difference spectrum (subtracted in MATLAB)
  2969. % plot(fooof_freqs, fooofed_difference_spectrum, '--', 'LineWidth', 1.5, 'DisplayName', 'fooofedSpect minus apfit');
  2970. % Set logarithmic scale for frequency
  2971. % set(gca, 'XScale', 'log'); % Logarithmic x-axis
  2972. % set(gca, 'YScale', 'log'); % Logarithmic y-axis
  2973. % Add labels, legend, and grid
  2974. xlabel('Frequency (Hz)');
  2975. ylabel('Power');
  2976. legend_handle = legend('show');
  2977. % legend('show');
  2978. grid on;
  2979. % Position the legend in the middle right of the plot
  2980. set(legend_handle, 'Location', 'east');
  2981. fooof_figname = ['S' num2str(SubjG{iStimGroup}(subjectNumInStimGroup)) '_FOOOF_' EEGch{singleFOOOFEEGCh} '_' num2str(f_range(1)) '-' num2str(f_range(2)) ' Hz'];
  2982. title(['S' num2str(SubjG{iStimGroup}(subjectNumInStimGroup)) ' (' GroupName{iStimGroup} ') - Ch:' EEGch{singleFOOOFEEGCh} ' - ' num2str(f_range(1)) '-' num2str(f_range(2)) ' Hz']);
  2983. hold off;
  2984. % Save the figure as a single file
  2985. saveas(gcf, [fooof_save_folder fooof_figname '.png']);
  2986. %% FOOOF PSD - prep for fig2c
  2987. % set up fooof input variables:
  2988. allFreqs = Fplot; % = row vector of frequency values
  2989. f_range = [2 45]; % f_range = fitting range
  2990. % set default settings values - % settings = fooof model settings, in a struct
  2991. settings = struct(...
  2992. 'peak_width_limits', [0.5, 12], ...
  2993. 'max_n_peaks', Inf, ...
  2994. 'min_peak_height', 0.0, ...
  2995. 'peak_threshold', 2.0, ...
  2996. 'aperiodic_mode', 'fixed', ...
  2997. 'verbose', true);
  2998. return_model = 1; % return_model = boolean of whether to return the FOOOF model fit, optional
  2999. trialType = 3;
  3000. input_power_spectrum = PSDall(SubjG{iStimGroup}(1),:,1,trialType);
  3001. dummy_fooof = fooof(Fplot, input_power_spectrum, f_range, settings, return_model);
  3002. dummy_fooof.difference_spectrum = dummy_fooof.power_spectrum - dummy_fooof.ap_fit;
  3003. nSubjPSDs = length(PSDall(:,1,1,3)); % effectively gets the total number of EEG subjects that have a PSD
  3004. nCh = length(ChanEEGLab);
  3005. emptyStruct = struct(); % Create an empty structure
  3006. clear allFooofResults allFoooFDiffPSD
  3007. allFooofResults = repmat(dummy_fooof, nSubjPSDs, nCh); % Replicate the empty structure 32 times
  3008. for iSub = 1:nSubjPSDs
  3009. for iCh=1:length(ChanEEGLab)
  3010. input_power_spectrum = PSDall(iSub,:,iCh,3); % power_spectrum = row vector of power values
  3011. % PSDall(participants, frequencies, channels, trialtypes)
  3012. % run fooof:
  3013. iSubiChFoofResults = fooof(Fplot, input_power_spectrum, f_range, settings, return_model);
  3014. % Extract data from the results structure
  3015. fooof_freqs = iSubiChFoofResults.freqs; % Frequency values
  3016. power_spectrum = iSubiChFoofResults.power_spectrum; % Original power spectrum
  3017. fooofed_spectrum = iSubiChFoofResults.fooofed_spectrum; % Full FOOOF fit
  3018. ap_fit = iSubiChFoofResults.ap_fit; % Aperiodic fit
  3019. iSubiChFoofResults.difference_spectrum = power_spectrum - ap_fit;
  3020. allFooofResults(iCh) = iSubiChFoofResults;
  3021. allFooofDiffPSD(iSub,:,iCh,trialType) = iSubiChFoofResults.difference_spectrum;
  3022. end
  3023. end
  3024. %% Plot average FOOOF results:
  3025. % % Extract data from the results structure
  3026. % fooof_freqs = iChFooofResults.freqs; % Frequency values
  3027. % power_spectrum = iChFooofResults.power_spectrum; % Original power spectrum
  3028. % fooofed_spectrum = iChFooofResults.fooofed_spectrum; % Full FOOOF fit
  3029. % ap_fit = iChFooofResults.ap_fit; % Aperiodic fit
  3030. % difference_spectrum = power_spectrum - ap_fit;
  3031. % % fooofed_difference_spectrum = fooofed_spectrum - ap_fit;
  3032. %
  3033. % % Plot settings
  3034. % figure;
  3035. % hold on;
  3036. %
  3037. % % Plot the original power spectrum
  3038. % plot(fooof_freqs, power_spectrum, 'k', 'LineWidth', 1.5, 'DisplayName', 'Power Spectrum');
  3039. %
  3040. % % Plot the FOOOFed spectrum (from fooof results)
  3041. % plot(fooof_freqs, fooofed_spectrum, 'r', 'LineWidth', 1.5, 'DisplayName', 'FOOOFed Spectrum');
  3042. %
  3043. % % Plot the aperiodic fit
  3044. % plot(fooof_freqs, ap_fit, 'b--', 'LineWidth', 1.5, 'DisplayName', 'Aperiodic Fit');
  3045. %
  3046. % % Plot the difference spectrum (subtracted in MATLAB)
  3047. % plot(fooof_freqs, difference_spectrum, 'g--', 'LineWidth', 1.5, 'DisplayName', 'powerSpect minus apfit');
  3048. %
  3049. % % % Plot the FOOOFed difference spectrum (subtracted in MATLAB)
  3050. % % plot(fooof_freqs, fooofed_difference_spectrum, '--', 'LineWidth', 1.5, 'DisplayName', 'fooofedSpect minus apfit');
  3051. %
  3052. % % Set logarithmic scale for frequency
  3053. % % set(gca, 'XScale', 'log'); % Logarithmic x-axis
  3054. % % set(gca, 'YScale', 'log'); % Logarithmic y-axis
  3055. %
  3056. % % Add labels, legend, and grid
  3057. % xlabel('Frequency (Hz)');
  3058. % ylabel('Power');
  3059. % legend('show');
  3060. % grid on;
  3061. % title(['FOOOF Analysis Results: ' num2str(f_range(1)) '-' num2str(f_range(2)) ' Hz']);
  3062. % hold off;
  3063. %%
  3064. % Predefine the number of stimulation groups, subjects, and channels
  3065. nStimGroups = 3;
  3066. nFreqs = length(fooof_freqs);
  3067. nCh = 32;
  3068. % Preallocate for average difference spectra per channel and stim group
  3069. avgDiffSpectraChannels = zeros(nStimGroups, nFreqs, nCh);
  3070. % Compute the average difference spectrum for each channel and stim group
  3071. for iStimGroup = 1:nStimGroups
  3072. for iCh = 1:nCh
  3073. % Extract difference spectra for the current stim group and channel
  3074. diffSpectraGroup = allFooofDiffPSD(SubjG{iStimGroup},:,iCh,trialtype); % Subj x Freq
  3075. % Average across subjects
  3076. avgDiffSpectraChannels(iStimGroup, :, iCh) = mean(diffSpectraGroup, 1, 'omitnan');
  3077. end
  3078. end
  3079. % Plot the average difference spectra for each channel, one plot per stim group
  3080. for iStimGroup = 1:nStimGroups
  3081. figure;
  3082. hold on;
  3083. colors = lines(nCh); % Generate distinct colors for each channel
  3084. for iCh = 1:nCh
  3085. plot(fooof_freqs, avgDiffSpectraChannels(iStimGroup, :, iCh), 'LineWidth', 1.5, ...
  3086. 'DisplayName', EEGch{iCh}, 'Color', colors(iCh, :));
  3087. end
  3088. % Customize the plot
  3089. xlabel('Frequency (Hz)');
  3090. ylabel('Difference Spectrum Power');
  3091. title(['Average Difference Spectrum by Channel - Stim Group ' GroupName{iStimGroup}]);
  3092. legend('show', 'Location', 'eastoutside');
  3093. grid on;
  3094. hold off;
  3095. % Save the figure
  3096. saveas(gcf, [fooof_save_folder 'Avg_Diff_Spectrum_by_Channel_StimGroup' GroupName{iStimGroup} '.png']);
  3097. end
  3098. %% FOOOFed Spectrum Plot + settings
  3099. figure;
  3100. hold on;
  3101. % Plot the original power spectrum
  3102. plot(fooof_freqs, power_spectrum, 'k', 'LineWidth', 1.5, 'DisplayName', 'Power Spectrum');
  3103. % Plot the FOOOFed spectrum (from fooof results)
  3104. % plot(fooof_freqs, fooofed_spectrum, 'r', 'LineWidth', 1.5, 'DisplayName', 'FOOOFed Spectrum');
  3105. % Plot the aperiodic fit
  3106. plot(fooof_freqs, ap_fit, 'b--', 'LineWidth', 1.5, 'DisplayName', 'Aperiodic Fit');
  3107. % Plot the difference spectrum (subtracted in MATLAB)
  3108. plot(fooof_freqs, difference_spectrum, 'g--', 'LineWidth', 1.5, 'DisplayName', 'Adjusted Spectrum (PS-ApFit)');
  3109. % % Plot the FOOOFed difference spectrum (subtracted in MATLAB)
  3110. % plot(fooof_freqs, fooofed_difference_spectrum, '--', 'LineWidth', 1.5, 'DisplayName', 'fooofedSpect minus apfit');
  3111. % Set logarithmic scale for frequency
  3112. % set(gca, 'XScale', 'log'); % Logarithmic x-axis
  3113. % set(gca, 'YScale', 'log'); % Logarithmic y-axis
  3114. % Add labels, legend, and grid
  3115. xlabel('Frequency (Hz)');
  3116. ylabel('Power');
  3117. legend_handle = legend('show');
  3118. % legend('show');
  3119. grid on;
  3120. % Position the legend in the middle right of the plot
  3121. set(legend_handle, 'Location', 'east');
  3122. fooof_figname = ['S' num2str(SubjG{iStimGroup}(subjectNumInStimGroup)) '_FOOOF_' EEGch{singleFOOOFEEGCh} '_' num2str(f_range(1)) '-' num2str(f_range(2)) ' Hz'];
  3123. title(['S' num2str(SubjG{iStimGroup}(subjectNumInStimGroup)) ' (' GroupName{iStimGroup} ') - Ch:' EEGch{singleFOOOFEEGCh} ' - ' num2str(f_range(1)) '-' num2str(f_range(2)) ' Hz']);
  3124. hold off;
  3125. % Save the figure as a single file
  3126. saveas(gcf, [fooof_save_folder fooof_figname '.png']);
  3127. %% PSD topographical Plots on Brain (heatmaps) - group difference, PSD-Acc correlation brain *** fig2c fig2d fig3
  3128. PowerLim=[-7 -2];
  3129. PowerLab={'-7' '-2'};
  3130. TLim=[-5 5];
  3131. TLab={'-5' '5'};
  3132. PLim=[-4 4];
  3133. PLab={'10e-4' '10e-0'};
  3134. fooof_freqs = firstsubjectfirstCh.freqs; % Frequency values
  3135. for iFreqFunc=3 %1:length(FreqFuncNames)
  3136. clear MapGroup diffTmap diffMap Tdata;
  3137. clear diffTPmap diffRPmap diffTmap
  3138. for iTrialType=3%1:length(TrialType)
  3139. % iCom=1;
  3140. PSDsaveDate = datestr(datetime, 'yy-mm-dd_HHMMSSFFF');
  3141. SaveTemp=['Fig2Results\' PSDsaveDate '\'];
  3142. SaveTemp=[SaveTemp FreqFuncNames{iFreqFunc} '\' ];
  3143. mkdir(SaveTemp)
  3144. BOI=[5;15];
  3145. BOI=[1 4 8 13 30 39.5 43;4 8 13 30 37 41.5 100];
  3146. BandName={'Delta','Theta','Alpha','Beta','Gamma-1','Gamma-E','Gamma-2'};
  3147. BandHzName={'1-4 Hz','4-8 Hz','8-13Hz','13-30 Hz','30-37 Hz','39-41 Hz','43-100 Hz'};
  3148. %
  3149. %% PSD on brain ***old Fig2C
  3150. % figure;
  3151. for iBOI=1:size(BOI,2)
  3152. NeedI=find(Fplot>=BOI(1,iBOI)&Fplot<=BOI(2,iBOI));
  3153. for iStimGroup=1:length(SubjG)
  3154. if iFreqFunc==3
  3155. Tdata{iStimGroup,iBOI,iTrialType}=squeeze(FreqFunc{iFreqFunc}(LogPSD(SubjG{iStimGroup},NeedI,EEGchInd,iTrialType),[],2));
  3156. MapGroup{iStimGroup,iBOI,iTrialType}=squeeze(nanmean(FreqFunc{iFreqFunc}(LogPSD(SubjG{iStimGroup},NeedI,EEGchInd,iTrialType),[],2),1));
  3157. else
  3158. Tdata{iStimGroup,iBOI,iTrialType}=squeeze(FreqFunc{iFreqFunc}(LogPSD(SubjG{iStimGroup},NeedI,EEGchInd,iTrialType),2));
  3159. MapGroup{iStimGroup,iBOI,iTrialType}=squeeze(nanmean(FreqFunc{iFreqFunc}(LogPSD(SubjG{iStimGroup},NeedI,EEGchInd,iTrialType),2),1));
  3160. end
  3161. subplotLU(2,size(BOI,2),iStimGroup,iBOI); %
  3162. [~,~,~,xmesh,ymesh]=topoplot(MapGroup{iStimGroup,iBOI,iTrialType}, ChanEEGLab,'colormap',jet,'maplimits',PowerLim);
  3163. if iStimGroup==2
  3164. xlabel(BandName{iBOI});
  3165. text(0,-0.55,[BandName{iBOI} ' (' BandHzName{iBOI} ')'],'horizontalalignment','center','verticalalignment','top','fontsize',10)
  3166. end
  3167. if iBOI==1
  3168. ylabel(GroupName{iStimGroup})
  3169. yt=text(-0.55,0,GroupName{iStimGroup},'horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  3170. end
  3171. end
  3172. end
  3173. % subplot('position',[0.5 0.51 0.3 0.01]);
  3174. % b=colorbar('southoutside');
  3175. % set(gca,'xtick',[],'ytick',[])
  3176. % set(b,'position',[0.5 0.5 0.3 0.03],'Limits',[0 1],'Ticks',[0 1],'Ticklabels',PowerLab);
  3177. % xlabel(b,'Log Normalized Power')
  3178. % LuFontStandard;
  3179. % papersizePX=[0 0 6*size(BOI,2) 6*2+3];
  3180. % set(gcf, 'PaperUnits', 'centimeters');
  3181. % set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  3182. % saveas(gcf,[SaveTemp 'GroupPSDonBrain'],'pdf');
  3183. % saveas(gcf,[SaveTemp 'GroupPSDonBrain'],'png');
  3184. % saveas(gcf,[SaveTemp 'GroupPSDonBrain.eps'],'epsc');
  3185. %% FOOOF PSD on brain *** old Fig2C
  3186. % % Ensure the desired colormap is loaded
  3187. % % batlow is part of the cmocean or scientific colormaps package
  3188. % if exist('batlow', 'file') == 2
  3189. % cmap = batlow; % Load batlow if available
  3190. % else
  3191. % cmap = parula; % Fallback to parula if batlow isn't available
  3192. % end
  3193. %
  3194. % figure;
  3195. % BOI=[1 4 8 13 30 39.5 43;4 8 13 30 37 41.5 100];
  3196. % BandName={'Delta','Theta','Alpha','Beta','Gamma-1','Gamma-E','Gamma-2'};
  3197. % BandHzName={'1-4 Hz','4-8 Hz','8-13Hz','13-30 Hz','30-37 Hz','39-41 Hz','43-100 Hz'};
  3198. % PowerLim = [-0.5 1];
  3199. % for iBOI=1:size(BOI,2)
  3200. % % NeedI=find(Fplot>=BOI(1,iBOI)&Fplot<=BOI(2,iBOI));
  3201. % NeedI=find(fooof_freqs>=BOI(1,iBOI)&fooof_freqs<=BOI(2,iBOI))';
  3202. %
  3203. % for iStimGroup=1:length(SubjG)
  3204. % if iFreqFunc==3
  3205. % FooofTdata{iStimGroup,iBOI,iTrialType}=squeeze(FreqFunc{iFreqFunc}(allFooofDiffPSD(SubjG{iStimGroup},NeedI,EEGchInd,iTrialType),[],2));
  3206. % FooofMapGroup{iStimGroup,iBOI,iTrialType}=squeeze(nanmean(FreqFunc{iFreqFunc}(allFooofDiffPSD(SubjG{iStimGroup},NeedI,EEGchInd,iTrialType),[],2),1));
  3207. %
  3208. % else
  3209. % FooofTdata{iStimGroup,iBOI,iTrialType}=squeeze(FreqFunc{iFreqFunc}(allFooofDiffPSD(SubjG{iStimGroup},NeedI,EEGchInd,iTrialType),2));
  3210. % FooofMapGroup{iStimGroup,iBOI,iTrialType}=squeeze(nanmean(FreqFunc{iFreqFunc}(allFooofDiffPSD(SubjG{iStimGroup},NeedI,EEGchInd,iTrialType),2),1));
  3211. % end
  3212. % subplotLU(3,size(BOI,2),iStimGroup,iBOI); %
  3213. % [~,~,~,xmesh,ymesh]=topoplot(FooofMapGroup{iStimGroup,iBOI,iTrialType}, ChanEEGLab,'colormap',parula,'maplimits',PowerLim);
  3214. % if iStimGroup==3
  3215. % xlabel(BandName{iBOI});
  3216. % text(0,-0.55,[BandName{iBOI} ' (' BandHzName{iBOI} ')'],'horizontalalignment','center','verticalalignment','top','fontsize',10)
  3217. %
  3218. % end
  3219. % if iBOI==1
  3220. % ylabel(GroupName{iStimGroup})
  3221. % yt=text(-0.55,0,GroupName{iStimGroup},'horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  3222. % end
  3223. % end
  3224. % end
  3225. % % subplot('position',[0.5 0.51 0.3 0.01]);
  3226. % b=colorbar('southoutside');
  3227. % set(gca,'xtick',[],'ytick',[])
  3228. % % set(b,'position',[0.5 0.5 0.3 0.03],'Limits',[0 1],'Ticks',[0 1],'Ticklabels',PowerLab);
  3229. % xlabel(b,'(not Log Normalized) Power')
  3230. % LuFontStandard;
  3231. % papersizePX=[0 0 6*size(BOI,2) 6*2+3];
  3232. % set(gcf, 'PaperUnits', 'centimeters');
  3233. % set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  3234. % % saveas(gcf,[SaveTemp 'FOOOFPSD'],'pdf');
  3235. % saveas(gcf,[SaveTemp 'FOOOFPSD_oldColorbarRange'],'png');
  3236. % % saveas(gcf,[SaveTemp 'PSD40BothControls.eps'],'epsc');
  3237. %% Updated FOOOF color map 1/22/25 *** Fig2C
  3238. % % Compute global 10% minimum and 90% maximum across all data
  3239. % allData = []; % Initialize an empty array to collect all data values
  3240. % for iBOI = 1:size(BOI, 2)
  3241. % for iStimGroup = 1:length(SubjG)
  3242. % % Collect all FooofMapGroup data into a single array
  3243. % allData = [allData; FooofMapGroup{iStimGroup, iBOI, iTrialType}(:)];
  3244. % end
  3245. % end
  3246. %
  3247. % % Compute 10% and 90% percentiles
  3248. % cbarMin = prctile(allData, 10);
  3249. % cbarMax = prctile(allData, 90);
  3250. %
  3251. % % Update PowerLim based on the computed values
  3252. % FOOOFPowerLim = [cbarMin, cbarMax];
  3253. %
  3254. % figure;
  3255. % for iBOI = 1:size(BOI, 2)
  3256. % NeedI = find(fooof_freqs >= BOI(1, iBOI) & fooof_freqs <= BOI(2, iBOI))';
  3257. %
  3258. % for iStimGroup = 1:3 %length(SubjG)
  3259. % subplotLU(3, size(BOI, 2), iStimGroup, iBOI);
  3260. % [~, ~, ~, xmesh, ymesh] = topoplot(FooofMapGroup{iStimGroup, iBOI, iTrialType}, ChanEEGLab, ...
  3261. % 'colormap', parula, 'maplimits', FOOOFPowerLim);
  3262. %
  3263. % if iStimGroup == 3
  3264. % xlabel(BandName{iBOI});
  3265. % text(0, -0.55, [BandName{iBOI} ' (' BandHzName{iBOI} ')'], ...
  3266. % 'horizontalalignment', 'center', 'verticalalignment', 'top', 'fontsize', 10);
  3267. % end
  3268. % if iBOI == 1
  3269. % ylabel(GroupName{iStimGroup});
  3270. % text(-0.55, 0, GroupName{iStimGroup}, ...
  3271. % 'horizontalalignment', 'center', 'verticalalignment', 'bottom', ...
  3272. % 'fontsize', 10, 'rotation', 90);
  3273. % end
  3274. % end
  3275. % end
  3276. %
  3277. % % b = colorbar('southoutside');
  3278. % % caxis(FOOOFPowerLim); % Set colorbar limits to match PowerLim
  3279. % % xlabel(b, '(not Log Normalized) Power');
  3280. %
  3281. % LuFontStandard;
  3282. % papersizePX = [0 0 6 * size(BOI, 2) 6 * 2 + 3];
  3283. % set(gcf, 'PaperUnits', 'centimeters');
  3284. % set(gcf, 'PaperPosition', papersizePX, 'PaperSize', papersizePX(3:4));
  3285. % saveas(gcf, [SaveTemp 'FOOOFPSD_updatedColorbarRange'], 'png');
  3286. %
  3287. %% Plot and save FOOOF vertical colorbar separately *** Fig2C
  3288. % figure;
  3289. %
  3290. % % Create a dummy image to generate a colorbar
  3291. % imagesc([0 1]); % Placeholder data
  3292. % colormap(parula); % Use the same colormap
  3293. % caxis(FOOOFPowerLim); % Apply the same color limits
  3294. %
  3295. % % Customize colorbar
  3296. % b = colorbar('eastoutside'); % Set colorbar orientation to vertical
  3297. %
  3298. % % Set ticks and format tick labels with 2 significant figures
  3299. % tickValues = linspace(FOOOFPowerLim(1), FOOOFPowerLim(2), 5); % Generate 5 evenly spaced tick values
  3300. % set(b, 'Ticks', tickValues); % Set tick positions
  3301. % set(b, 'TickLabels', arrayfun(@(x) sprintf('%.2g', x), tickValues, 'UniformOutput', false)); % Format tick labels
  3302. %
  3303. % % Add label to the colorbar
  3304. % ylabel(b, '(not Log Normalized) Power', 'fontsize', 12, ...
  3305. % 'rotation', 270, 'VerticalAlignment', 'bottom', ...
  3306. % 'HorizontalAlignment', 'center');
  3307. %
  3308. % % Adjust figure layout to fit the colorbar and label
  3309. % set(gca, 'Visible', 'off'); % Hide axes
  3310. % set(gcf, 'PaperUnits', 'centimeters');
  3311. % set(gcf, 'PaperPosition', [0 0 2 10]); % Adjust to fit vertical colorbar
  3312. % set(gcf, 'PaperSize', [2 10]);
  3313. %
  3314. % % Save colorbar as a PNG file
  3315. % saveas(gcf, [SaveTemp 'FOOOF_Colorbar'], 'png');
  3316. %% Group differences on brain *** old fig2d ? suppfig2
  3317. figure;
  3318. for iBOI=1:size(BOI,2)
  3319. for iCh=1:length(ChanEEGLab)
  3320. [~,diffTPmap(iCh,iBOI,iTrialType),~,stats]=ttest2(Tdata{2,iBOI,iTrialType}(:,iCh),Tdata{1,iBOI,iTrialType}(:,iCh));
  3321. [diffRPmap(iCh,iBOI,iTrialType),~,~]=ranksum(Tdata{2,iBOI,iTrialType}(:,iCh),Tdata{1,iBOI,iTrialType}(:,iCh));
  3322. diffTmap(iCh,iBOI,iTrialType)=stats.tstat;
  3323. end
  3324. subplotLU(2,size(BOI,2),1,iBOI);
  3325. topoplot(diffTmap(:,iBOI,iTrialType), ChanEEGLab,'colormap',colorMapPN,'maplimits',TLim);
  3326. if iBOI==1
  3327. yt=text(-0.55,0,['T, ' GroupName{Group1} '-' GroupName{Group2}],'horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  3328. end
  3329. if iBOI==size(BOI,2)
  3330. b=colorbar('southoutside');
  3331. % set(b,'position',[0.52 0.93 0.2 0.03],'xtick',[-6 6],'xticklabel',{'-6' '6'},'xlim',[-6 6]);
  3332. set(gca,'xtick',[],'ytick',[])
  3333. set(b,'position',[0.52 0.93 0.2 0.01],'ticks',TLim,'ticklabels',TLab);
  3334. xlabel(b,'T statistics','verticalalignment','top')
  3335. end
  3336. subplotLU(2,size(BOI,2),2,iBOI);
  3337. topoplot(log10(diffTPmap(:,iBOI,iTrialType)), ChanEEGLab,'colormap',colorMapPN,'maplimits',PLim);
  3338. xlabel(BandName{iBOI});
  3339. text(0,-0.55,[BandName{iBOI} ' (' BandHzName{iBOI} ')'],'horizontalalignment','center','verticalalignment','top','fontsize',10)
  3340. if iBOI==1
  3341. a=ylabel('40Hz-BothControls');
  3342. % a.Position=[0.01 0.5 0.03 0.4];
  3343. % a.verticalalignment='middle';
  3344. % set(a,'Position',[0.01 0.5 0.03 0.4],'Verticalalignment','middle')
  3345. set(a,'Verticalalignment','middle')
  3346. yt=text(-0.55,0,['P, ' GroupName{Group1} '-' GroupName{Group2}],'horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  3347. end
  3348. if iBOI==size(BOI,2)
  3349. c=colorbar('southoutside');
  3350. % set(c,'position',[0.52 0.46 0.2 0.03],'xtick',[-4 0],'xticklabel',{'10e-4' '10e-0'},'xlim',[-4 0]);
  3351. set(gca,'xtick',[],'ytick',[])
  3352. set(c,'position',[0.52 0.51 0.2 0.01],'Limits',[PLim(1) 0],'ticks',[PLim(1) 0],'ticklabels',PLab);
  3353. xlabel(c,'P values','verticalalignment','top')
  3354. end
  3355. end
  3356. LuFontStandard;
  3357. papersizePX=[0 0 6*size(BOI,2) 6*2+3];
  3358. set(gcf, 'PaperUnits', 'centimeters');
  3359. set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  3360. saveas(gcf,[SaveTemp 'PSDDiff40BothControls'],'pdf');
  3361. saveas(gcf,[SaveTemp 'PSDDiff40BothControls'],'png');
  3362. saveas(gcf,[SaveTemp 'PSDDiff40BothControls.eps'],'epsc');
  3363. %% Group differences on brain - Calculate t-test stat and plot *** fig2d - suppfig2
  3364. % Define group pairs to compare % see: open: GroupName
  3365. groupPairs = [
  3366. 1, 2; % First pair: Group1 = 1 (40 Hz), Group2 = 2 (Light)
  3367. 1, 3 % Second pair: Group1 = 1, Group2 = 3 (Random)
  3368. ];
  3369. % Channel subset Fp1, Cz, Oz
  3370. % Channels of interest (indices)
  3371. PSDttestChSubset = [1, 32, 16]; % [1, 32, 16] = [Fp1, Cz, Oz]
  3372. % Initialize a matrix to store p-values for FDR correction
  3373. allPValues = [];
  3374. % Loop through each group pair
  3375. for iPair = 1:size(groupPairs, 1)
  3376. Group1 = groupPairs(iPair, 1);
  3377. Group2 = groupPairs(iPair, 2);
  3378. dataName = [TrialTypeName{iTrialType} FreqFuncNames{iFreqFunc} GroupName{Group1} GroupName{Group2}];
  3379. figure;
  3380. for iBOI=1:size(BOI,2)
  3381. %% calculate stats test
  3382. for iCh=1:length(ChanEEGLab)
  3383. [~,diffTPmap(iCh,iBOI,iTrialType),~,stats]=ttest2(Tdata{Group1,iBOI,iTrialType}(:,iCh),Tdata{Group2,iBOI,iTrialType}(:,iCh)); % two-sided ttest
  3384. [diffRPmap(iCh,iBOI,iTrialType),~,~]=ranksum(Tdata{Group1,iBOI,iTrialType}(:,iCh),Tdata{Group2,iBOI,iTrialType}(:,iCh));
  3385. diffTmap(iCh,iBOI,iTrialType)=stats.tstat;
  3386. end
  3387. % Collect p-values for channels of interest
  3388. for chIdx = 1:length(PSDttestChSubset)
  3389. iCh = PSDttestChSubset(chIdx);
  3390. % allPValues = [p-value, t-stat, groupPair, BOI, Ch]
  3391. allPValues(end+1, :) = [diffTPmap(iCh, iBOI, iTrialType), diffTmap(iCh, iBOI, iTrialType), iPair, iBOI, iCh]; % Store p-value with metadata
  3392. end
  3393. %% Plot T-stat
  3394. subplotLU(2,size(BOI,2),1,iBOI);
  3395. topoplot(diffTmap(:,iBOI,iTrialType), ChanEEGLab,'colormap',colorMapPN,'maplimits',TLim);
  3396. if iBOI==1
  3397. yt=text(-0.55,0,['T, ' GroupName{Group1} '-' GroupName{Group2}],'horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  3398. end
  3399. if iBOI==size(BOI,2)
  3400. b=colorbar('southoutside');
  3401. % set(b,'position',[0.52 0.93 0.2 0.03],'xtick',[-6 6],'xticklabel',{'-6' '6'},'xlim',[-6 6]);
  3402. set(gca,'xtick',[],'ytick',[])
  3403. set(b,'position',[0.52 0.93 0.2 0.01],'ticks',TLim,'ticklabels',TLab);
  3404. xlabel(b,'T statistics','verticalalignment','top')
  3405. end
  3406. %% Plot p-value
  3407. subplotLU(2,size(BOI,2),2,iBOI);
  3408. topoplot(log10(diffTPmap(:,iBOI,iTrialType)), ChanEEGLab,'colormap',colorMapPN,'maplimits',PLim);
  3409. xlabel(BandName{iBOI});
  3410. text(0,-0.55,[BandName{iBOI} ' (' BandHzName{iBOI} ')'],'horizontalalignment','center','verticalalignment','top','fontsize',10)
  3411. if iBOI==1
  3412. a=ylabel([GroupName{Group1} '-' GroupName{Group2}]);
  3413. % a.Position=[0.01 0.5 0.03 0.4];
  3414. % a.verticalalignment='middle';
  3415. % set(a,'Position',[0.01 0.5 0.03 0.4],'Verticalalignment','middle')
  3416. set(a,'Verticalalignment','middle')
  3417. yt=text(-0.55,0,['P, ' GroupName{Group1} '-' GroupName{Group2}],'horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  3418. end
  3419. if iBOI==size(BOI,2)
  3420. c=colorbar('southoutside');
  3421. % set(c,'position',[0.52 0.46 0.2 0.03],'xtick',[-4 0],'xticklabel',{'10e-4' '10e-0'},'xlim',[-4 0]);
  3422. set(gca,'xtick',[],'ytick',[])
  3423. set(c,'position',[0.52 0.51 0.2 0.01],'Limits',[PLim(1) 0],'ticks',[PLim(1) 0],'ticklabels',PLab);
  3424. xlabel(c,'P values','verticalalignment','top')
  3425. end
  3426. end
  3427. LuFontStandard;
  3428. papersizePX=[0 0 6*size(BOI,2) 6*2+3];
  3429. set(gcf, 'PaperUnits', 'centimeters');
  3430. set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  3431. saveas(gcf,[SaveTemp 'PSDDiff' GroupName{Group1} '-' GroupName{Group2}],'pdf');
  3432. saveas(gcf,[SaveTemp 'PSDDiff' GroupName{Group1} '-' GroupName{Group2}],'png');
  3433. saveas(gcf,[SaveTemp 'PSDDiff' GroupName{Group1} '-' GroupName{Group2} '.eps'],'epsc');
  3434. end
  3435. %% PSD ttest FDR Correction fig2d stats
  3436. % Filter p-values for the first 5 bands
  3437. numBandsToInclude = 5; % Only include the first 5 bands
  3438. % allPValues = [p-value, t-stat, groupPiar, BOI, Ch]
  3439. filteredPValues = allPValues(allPValues(:, 4) <= numBandsToInclude, :);
  3440. % Extract only the raw p-values for FDR correction
  3441. pValuesForFDR = filteredPValues(:, 1);
  3442. % % Perform FDR Correction (three different methods below)
  3443. BHFDR = mafdr(pValuesForFDR, 'BHFDR', true); % Benjamini-Hochberg FDR correction
  3444. [FDRST,Q,aPrioriProb] = mafdr(pValuesForFDR); % Storey-Tibshirani method (2002)
  3445. % fdh_bh, parameters - gets the critical p value
  3446. q=0.05;
  3447. method='pdep';
  3448. report='yes';
  3449. [fdrbh_h, fdrbh_crit_p, fdrbh_adj_p]=fdr_bh(pValuesForFDR,q,method,report);
  3450. % Add FDR corrected p-values back to the filtered results
  3451. filteredPValuesPlusFDR = [filteredPValues, fdrAdjustedPValues];
  3452. % Display Results
  3453. fprintf('Channel\tGroup Pair\tBand\tT-stat\tRaw P-value\tFDR Corrected P-value\n');
  3454. for i = 1:size(filteredPValuesPlusFDR, 1)
  3455. % Get channel name from the struct
  3456. chName = ChanEEGLab(filteredPValuesPlusFDR(i, 5)).labels; % Access the 'labels' field of the struct
  3457. % Extract group indices
  3458. group1Idx = groupPairs(filteredPValuesPlusFDR(i, 3), 1);
  3459. group2Idx = groupPairs(filteredPValuesPlusFDR(i, 3), 2);
  3460. % Get group pair names
  3461. groupPair = sprintf('%s-%s', GroupName{group1Idx}, GroupName{group2Idx});
  3462. % Get band name
  3463. band = BandName{filteredPValuesPlusFDR(i, 4)}; % Assuming BandName is a cell array
  3464. % Get t-test stat
  3465. tstat = filteredPValuesPlusFDR(i, 2);
  3466. % Get raw and FDR-corrected p-values
  3467. rawP = filteredPValuesPlusFDR(i, 1);
  3468. fdrP = filteredPValuesPlusFDR(i, 6);
  3469. % Print result
  3470. fprintf('%s\t%s\t%s\t%.4f\t%.4f\t%.4f\n', chName, groupPair, band, tstat, rawP, fdrP);
  3471. end
  3472. % Display total number of comparisons
  3473. numComparisons = size(filteredPValuesPlusFDR, 1);
  3474. fprintf('Total number of comparisons used for FDR correction: %d\n', numComparisons);
  3475. %% NEW PSD-ACC and PSD_RT 2025-02-04 MKA
  3476. %% % Set-up for PSD-Beh - Spearman PSD-Acc and PSD-RT ***fig3b & fig3c
  3477. SaveTemp = ['Fig3Results\' PSDsaveDate '\'];
  3478. mkdir(SaveTemp)
  3479. % PSD groups: 1=40Hz, 3=Random, 6=LightRT
  3480. PSDBehaviorGroup1 = 1;
  3481. PSDBehaviorGroup2 = 3;
  3482. PSDBehaviorGroup3 = 6;
  3483. allGroupNames = [GroupName{PSDBehaviorGroup1} GroupName{PSDBehaviorGroup2} GroupName{PSDBehaviorGroup3}];
  3484. dataName = [TrialTypeName{iTrialType} FreqFuncNames{iFreqFunc} allGroupNames];
  3485. IncludedSubj = union(SubjG{PSDBehaviorGroup1}, union(SubjG{PSDBehaviorGroup2}, SubjG{PSDBehaviorGroup3}));
  3486. figstatsChSubset = PSDttestChSubset;
  3487. % Create separate figures for PSD-Acc and PSD-RT
  3488. figureAcc = figure;
  3489. figureRT = figure;
  3490. Acctemp = Acc(IncludedSubj);
  3491. RTtemp = SubjsAvgRT(IncludedSubj);
  3492. %%% Variables for scatterplot:
  3493. FlickerID=[zeros(size(SubjG{PSDBehaviorGroup1}(:)))+1;zeros(size(SubjG{PSDBehaviorGroup2}(:)))+2;zeros(size(SubjG{PSDBehaviorGroup3}(:)))+3];
  3494. Param.Corr='Spearman'; %%%Type of correlation, see Matlab function corr for more details
  3495. Param.Pth=0.05; %%%threshold of Pvalue
  3496. Param.EdgeColor=colorMapPN; %%%Color map for correlation link
  3497. Param.NodeColor=repmat([0.8 0.8 0.8],6,1); %%%Color of Nodes for correlation link plot
  3498. Param.Clim=[-1 1]; %%%Color Limit for Correlation
  3499. Param.Title='Pool All Sample'; %%Any title for label the figure
  3500. Param.MarkerSize=8; %%%MarkerSize of scatter
  3501. Param.SubjIDColor=FlickerColor;
  3502. Param.SubjID=FlickerID;
  3503. %%% Loop thru BOI for Fig3b and Fig3c heatmaps & plot
  3504. for iBOI = 1:size(BOI, 2)
  3505. NeedI = find(Fplot >= BOI(1, iBOI) & Fplot <= BOI(2, iBOI));
  3506. if iFreqFunc == 3
  3507. PSDtemp = squeeze(FreqFunc{iFreqFunc}(LogPSD(IncludedSubj, NeedI, EEGchInd, iTrialType), [], 2));
  3508. else
  3509. PSDtemp = squeeze(FreqFunc{iFreqFunc}(LogPSD(IncludedSubj, NeedI, EEGchInd, iTrialType), 2));
  3510. end
  3511. %% Scatterplot of correlation - indiv behavior and PSD fig3a+
  3512. %% % PSD_Acc scatter
  3513. PSDAcc_scatterplot = figure;
  3514. tempName=EEGch;
  3515. tempName{end+1}='Acc';
  3516. multiCorr2GroupSubplot(6,6,[PSDtemp Acctemp],size(PSDtemp,2)+1,tempName,Param)
  3517. LuFontStandard;
  3518. papersizePX=[0 0 6*6 6*6];
  3519. set(gcf, 'PaperUnits', 'centimeters');
  3520. set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  3521. % saveas(gcf,[SaveTemp Param.Corr BName{iBOI} BName2{iBOI} 'PSDAcc_' allGroupNames],'pdf');
  3522. saveas(gcf,[SaveTemp Param.Corr BName{iBOI} BName2{iBOI} 'PSDAcc_' allGroupNames],'png');
  3523. saveas(gcf,[SaveTemp Param.Corr BName{iBOI} BName2{iBOI} 'PSDAcc_' allGroupNames '.eps'],'epsc');
  3524. %% % PSD_RT scatter
  3525. PSDRT_scatterplot = figure;
  3526. tempName=EEGch;
  3527. tempName{end+1}='RT';
  3528. multiCorr2GroupSubplot(6,6,[PSDtemp RTtemp],size(PSDtemp,2)+1,tempName,Param)
  3529. LuFontStandard;
  3530. papersizePX=[0 0 6*6 6*6];
  3531. set(gcf, 'PaperUnits', 'centimeters');
  3532. set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  3533. % saveas(gcf,[SaveTemp Param.Corr BName{iBOI} BName2{iBOI} 'PSDAcc_' allGroupNames],'pdf');
  3534. saveas(gcf,[SaveTemp Param.Corr BName{iBOI} BName2{iBOI} 'PSDRT_' allGroupNames],'png');
  3535. saveas(gcf,[SaveTemp Param.Corr BName{iBOI} BName2{iBOI} 'PSDRT_' allGroupNames '.eps'],'epsc');
  3536. %% Plot PSD-Acc and PSD-RT: R values and P Values topoplots
  3537. %%% Calculate correlations & pvalues
  3538. [PSDAcc_rSpear(:, iBOI, iTrialType), PSDAcc_pSpear(:, iBOI, iTrialType)] = corr(PSDtemp, Acctemp, 'type', 'spearman', 'rows', 'pairwise');
  3539. [PSDRT_rSpear(:, iBOI, iTrialType), PSDRT_pSpear(:, iBOI, iTrialType)] = corr(PSDtemp, RTtemp, 'type', 'spearman', 'rows', 'pairwise');
  3540. %% % Plot PSD-Acc Correlation
  3541. figure(figureAcc);
  3542. %%%% Plot R value correlation topoplot PSD-Acc(fig3b)
  3543. subplotLU(2,size(BOI,2),1,iBOI);
  3544. topoplot(PSDAcc_rSpear(:,iBOI,iTrialType), ChanEEGLab,'colormap',colorMapPN,'maplimits',[-1 1]); % Changed from 'maplimits',[-1 1]
  3545. if iBOI==1
  3546. yt=text(-0.55,0,'EEG-Behavior R','horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  3547. end
  3548. if iBOI==size(BOI,2)
  3549. b=colorbar('southoutside');set(b,'position',[0.52 0.93 0.2 0.01],'xtick',[-1 1],'xticklabel',{'-1' '1'},'xlim',[-1 1]);
  3550. xlabel(b,'PSD-Acc Correlation','verticalalignment','top')
  3551. end
  3552. %%%% Plot PSD-Acc P-value topoplot (fig3b supplement)
  3553. subplotLU(2,size(BOI,2),2,iBOI);
  3554. topoplot(log10(PSDAcc_pSpear(:,iBOI,iTrialType)), ChanEEGLab,'colormap',colorMapPN,'maplimits',[-4 4]);
  3555. xlabel(BandName{iBOI});
  3556. text(0,-0.55,[BandName{iBOI} ' (' BandHzName{iBOI} ')'],'horizontalalignment','center','verticalalignment','top','fontsize',10)
  3557. if iBOI==1
  3558. a=ylabel(allGroupNames);
  3559. % set(a,'position',[0.01 0.5 0.03 0.4],'verticalalignment','middle')
  3560. set(a,'verticalalignment','middle')
  3561. yt=text(-0.55,0,'P EEG-Behavior','horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  3562. end
  3563. if iBOI==size(BOI,2)
  3564. % c=colorbar('southoutside');set(c,'position',[0.52 0.46 0.2 0.03],'xtick',[-4 0],'xticklabel',{'10e-4' '10e-0'},'xlim',[-4 0]);
  3565. % xlabel(c,'P values','verticalalignment','top')
  3566. c=colorbar('southoutside');
  3567. set(gca,'xtick',[],'ytick',[])
  3568. set(c,'position',[0.52 0.51 0.2 0.01],'Limits',[PLim(1) 0],'ticks',[PLim(1) 0],'ticklabels',PLab);
  3569. xlabel(c,'P values','verticalalignment','top')
  3570. end
  3571. %% % Plot PSD-RT Correlation
  3572. figure(figureRT);
  3573. %%%% Plot PSD-RT R value (correlation) heatmap (fig3c)
  3574. subplotLU(2,size(BOI,2),1,iBOI);
  3575. topoplot(PSDRT_rSpear(:,iBOI,iTrialType), ChanEEGLab,'colormap',colorMapPN,'maplimits',[-1 1]); % Changed from 'maplimits',[-1 1]
  3576. if iBOI==1
  3577. yt=text(-0.55,0,'EEG-Behavior R','horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  3578. end
  3579. if iBOI==size(BOI,2)
  3580. b=colorbar('southoutside');set(b,'position',[0.52 0.93 0.2 0.01],'xtick',[-1 1],'xticklabel',{'-1' '1'},'xlim',[-1 1]);
  3581. xlabel(b,'PSD-RT Correlation','verticalalignment','top')
  3582. end
  3583. %%%% Plot PSD-RT P-value fig3c supp
  3584. subplotLU(2,size(BOI,2),2,iBOI);
  3585. topoplot(log10(PSDRT_pSpear(:,iBOI,iTrialType)), ChanEEGLab,'colormap',colorMapPN,'maplimits',[-4 4]);
  3586. xlabel(BandName{iBOI});
  3587. text(0,-0.55,[BandName{iBOI} ' (' BandHzName{iBOI} ')'],'horizontalalignment','center','verticalalignment','top','fontsize',10)
  3588. if iBOI==1
  3589. a=ylabel(allGroupNames);
  3590. % set(a,'position',[0.01 0.5 0.03 0.4],'verticalalignment','middle')
  3591. set(a,'verticalalignment','middle')
  3592. yt=text(-0.55,0,'P EEG-Behavior','horizontalalignment','center','verticalalignment','bottom','fontsize',10,'rotation',90);
  3593. end
  3594. %%%%% add colorbar
  3595. if iBOI==size(BOI,2)
  3596. % c=colorbar('southoutside');set(c,'position',[0.52 0.46 0.2 0.03],'xtick',[-4 0],'xticklabel',{'10e-4' '10e-0'},'xlim',[-4 0]);
  3597. % xlabel(c,'P values','verticalalignment','top')
  3598. c=colorbar('southoutside');
  3599. set(gca,'xtick',[],'ytick',[])
  3600. set(c,'position',[0.52 0.51 0.2 0.01],'Limits',[PLim(1) 0],'ticks',[PLim(1) 0],'ticklabels',PLab);
  3601. xlabel(c,'P values','verticalalignment','top')
  3602. end
  3603. end
  3604. %%% Apply formatting and save figures for heatmaps
  3605. LuFontStandard;
  3606. papersizePX = [0 0 6*size(BOI,2) 6*2+3];
  3607. set(figureAcc, 'PaperUnits', 'centimeters');
  3608. set(figureAcc, 'PaperPosition', papersizePX, 'PaperSize', papersizePX(3:4));
  3609. saveas(figureAcc, [SaveTemp 'SpearmanPSDAcc_' allGroupNames], 'png');
  3610. set(figureRT, 'PaperUnits', 'centimeters');
  3611. set(figureRT, 'PaperPosition', papersizePX, 'PaperSize', papersizePX(3:4));
  3612. saveas(figureRT, [SaveTemp 'SpearmanPSDRT_' allGroupNames], 'png');
  3613. close all;
  3614. %% Perform FDR correction for PSD_Acc and PSD_RT *** stats for fig3b and fig3c
  3615. % Define the subset of channels and bands
  3616. figstatsChSubset = [1, 32, 16]; % Indices for Fp1, Cz, Oz
  3617. numBandsToInclude = 5; % First 5 bands
  3618. % Initialize container for p-values
  3619. PSDAcc_allPValues = [];
  3620. PSDRT_allPValues = [];
  3621. %%% Extract p-values and r-values for specified channels and bands
  3622. for iBOI = 1:numBandsToInclude
  3623. for iCh = figstatsChSubset
  3624. % Get the p-value and r-value for the current channel and band
  3625. PSDAcc_pValue = PSDAcc_pSpear(iCh, iBOI, iTrialType);
  3626. PSDAcc_rValue = PSDAcc_rSpear(iCh, iBOI, iTrialType);
  3627. PSDRT_pValue = PSDRT_pSpear(iCh, iBOI, iTrialType);
  3628. PSDRT_rValue = PSDRT_rSpear(iCh, iBOI, iTrialType);
  3629. % Append to the list of all p-values with associated metadata
  3630. PSDAcc_allPValues = [PSDAcc_allPValues; PSDAcc_rValue, PSDAcc_pValue, iBOI, iCh];
  3631. PSDRT_allPValues = [PSDRT_allPValues; PSDRT_rValue, PSDRT_pValue, iBOI, iCh];
  3632. end
  3633. end
  3634. %%% Perform FDR correction
  3635. % Extract raw p-values for FDR correction
  3636. PSDAcc_pValuesForFDR = PSDAcc_allPValues(:, 2);
  3637. PSDRT_pValuesForFDR = PSDRT_allPValues(:, 2);
  3638. % Perform FDR correction
  3639. [PSDAcc_fdrAdjustedPValues,Q_Acc,aPrioriProb_Acc,R_squared_Acc] = mafdr(PSDAcc_pValuesForFDR, 'BHFDR', true);
  3640. [PSDRT_fdrAdjustedPValues,Q_RT,aPrioriProb_RT,R_squared_RT] = mafdr(PSDRT_pValuesForFDR, 'BHFDR', true);
  3641. % Append FDR-adjusted p-values
  3642. PSDAcc_allPValues = [PSDAcc_allPValues, PSDAcc_fdrAdjustedPValues];
  3643. PSDRT_allPValues = [PSDRT_allPValues, PSDRT_fdrAdjustedPValues];
  3644. %%% Display and Save Results for PSD_Acc
  3645. outputFile_Acc = fullfile(SaveTemp, 'PSD_Acc_FDR_corrected_stats.txt');
  3646. fid_Acc = fopen(outputFile_Acc, 'w');
  3647. fprintf(fid_Acc, 'Channel\tBand\tR-value\tRaw_P-value\tFDR Corrected P-value\n');
  3648. fprintf('FDR-corrected p-values for PSD-Acc:\n');
  3649. % Display total number of comparisons
  3650. numComparisons = size(PSDAcc_allPValues, 1);
  3651. fprintf('Total number of comparisons used for FDR correction: %d\n', numComparisons);
  3652. fprintf('Channel\tBand\tR-value\tRaw_P-value\tFDR Corrected P-value\n');
  3653. for i = 1:size(PSDAcc_allPValues, 1)
  3654. chName = ChanEEGLab(PSDAcc_allPValues(i, 4)).labels;
  3655. band = BandName{PSDAcc_allPValues(i, 3)};
  3656. rValue = PSDAcc_allPValues(i, 1);
  3657. rawP = PSDAcc_allPValues(i, 2);
  3658. fdrP = PSDAcc_allPValues(i, 5);
  3659. fprintf('%s\t%s\t%.4f\t%.4f\t%.4f\n', chName, band, rValue, rawP, fdrP);
  3660. fprintf(fid_Acc, '%s\t%s\t%.4f\t%.4f\t%.4f\n', chName, band, rValue, rawP, fdrP);
  3661. end
  3662. fclose(fid_Acc);
  3663. fprintf('PSD-Acc results saved to %s\n', outputFile_Acc);
  3664. %%% Display and Save Results for PSD_RT
  3665. outputFile_RT = fullfile(SaveTemp, 'PSD_RT_FDR_corrected_stats.txt');
  3666. fid_RT = fopen(outputFile_RT, 'w');
  3667. fprintf(fid_RT, 'Ch\tBand\tR-value\tRaw_P-value\tFDR Corrected P-value\n');
  3668. fprintf('FDR-corrected p-values for PSD-RT:\n');
  3669. % Display total number of comparisons
  3670. numComparisons = size(PSDRT_allPValues, 1);
  3671. fprintf('Total number of comparisons used for FDR correction: %d\n', numComparisons);
  3672. fprintf('Ch\tBand\tR-value\tRaw_P-value\tFDR Corrected P-value\n');
  3673. for i = 1:size(PSDRT_allPValues, 1)
  3674. chName = ChanEEGLab(PSDRT_allPValues(i, 4)).labels;
  3675. band = BandName{PSDRT_allPValues(i, 3)};
  3676. rValue = PSDRT_allPValues(i, 1);
  3677. rawP = PSDRT_allPValues(i, 2);
  3678. fdrP = PSDRT_allPValues(i, 5);
  3679. fprintf('%s\t%s\t%.4f\t%.4f\t%.4f\n', chName, band, rValue, rawP, fdrP);
  3680. fprintf(fid_RT, '%s\t%s\t%.4f\t%.4f\t%.4f\n', chName, band, rValue, rawP, fdrP);
  3681. end
  3682. fclose(fid_RT);
  3683. fprintf('PSD-RT results saved to %s\n', outputFile_RT);
  3684. end
  3685. end
  3686. %% PSD-Acc - Loop of correlation scatter plot? (future fig3a+b?)
  3687. FlickerColor=[31 125 184; 150 27 27; 219 129 50]/255; %40, Random, Light % blue, red, gold
  3688. for iFreqFunc=3%1:length(FunGroupName)
  3689. clear MapGroup diffTmap diffMap Tdata;
  3690. clear diffTPmap diffRPmap diffTmap
  3691. for iTrialType=3%1:length(TrialType)
  3692. % iCom=1;
  3693. SaveTemp=[SubSavePSD TrialTypeName{iTrialType} '\'];
  3694. SaveTemp=[SaveTemp FunGroupName{iFreqFunc} '\' ];
  3695. mkdir(SaveTemp)
  3696. BOI=[5;15];
  3697. BOI=[1 4 8 13 30 39.5 43;4 8 13 30 37 41.5 100];
  3698. BName={'Delta','Theta','Alpha','Beta','Gamma-1','Gamma-E','Gamma-2'};
  3699. BName2={'1-4 Hz','4-8 Hz','8-13Hz','13-30 Hz','30-37 Hz','39-41 Hz','43-100 Hz'};
  3700. PSDBehaviorGroup1 = 1;
  3701. PSDBehaviorGroup2 = 3;
  3702. PSDBehaviorGroup3 = 6; % 6 = LightRT
  3703. allGroupNames = [GroupName{PSDBehaviorGroup1} GroupName{PSDBehaviorGroup2} GroupName{PSDBehaviorGroup3}];
  3704. dataName = [TrialTypeName{iTrialType} FreqFuncNames{iFreqFunc} allGroupNames];
  3705. IncludedSubj=[SubjG{PSDBehaviorGroup1}(:);SubjG{PSDBehaviorGroup2}(:);SubjG{PSDBehaviorGroup3}(:)]; % 1 3 6 = 40, Random, Light
  3706. FlickerID=[zeros(size(SubjG{PSDBehaviorGroup1}(:)))+1;zeros(size(SubjG{PSDBehaviorGroup2}(:)))+2;zeros(size(SubjG{PSDBehaviorGroup3}(:)))+3];
  3707. Param.Corr='Spearman'; %%%Type of correlation, see Matlab function corr for more details
  3708. Param.Pth=0.05; %%%threshold of Pvalue
  3709. Param.EdgeColor=colorMapPN; %%%Color map for correlation link
  3710. Param.NodeColor=repmat([0.8 0.8 0.8],6,1); %%%Color of Nodes for correlation link plot
  3711. Param.Clim=[-1 1]; %%%Color Limit for Correlation
  3712. Param.Title='Pool All Sample'; %%Any title for label the figure
  3713. Param.MarkerSize=8; %%%MarkerSize of scatter
  3714. Param.SubjIDColor=FlickerColor;
  3715. Param.SubjID=FlickerID;
  3716. %% Scatterplot of correlation - indiv behavior and PSD
  3717. for iBOI=1:size(BOI,2)
  3718. NeedI=find(Fplot>=BOI(1,iBOI)&Fplot<=BOI(2,iBOI));
  3719. figure;
  3720. % for iCh=1:length(ChanEEGLab)
  3721. if iFreqFunc==3
  3722. PSDtemp=squeeze(FunGroup{iFreqFunc}(LogPSD(IncludedSubj,NeedI,EEGchInd,iTrialType),[],2));
  3723. else
  3724. PSDtemp=squeeze(FunGroup{iFreqFunc}(LogPSD(IncludedSubj,NeedI,EEGchInd,iTrialType),2));
  3725. end
  3726. Acctemp=Acc(IncludedSubj);
  3727. tempName=EEGch;
  3728. tempName{end+1}='Acc';
  3729. multiCorr2GroupSubplot(6,6,[PSDtemp Acctemp],size(PSDtemp,2)+1,tempName,Param)
  3730. LuFontStandard;
  3731. papersizePX=[0 0 6*6 6*6];
  3732. set(gcf, 'PaperUnits', 'centimeters');
  3733. set(gcf,'PaperPosition',papersizePX,'PaperSize',papersizePX(3:4));
  3734. % saveas(gcf,[SaveTemp Param.Corr BName{iBOI} BName2{iBOI} 'PSDAcc_' allGroupNames],'pdf');
  3735. saveas(gcf,[SaveTemp Param.Corr BName{iBOI} BName2{iBOI} 'PSDAcc_' allGroupNames],'png');
  3736. saveas(gcf,[SaveTemp Param.Corr BName{iBOI} BName2{iBOI} 'PSDAcc_' allGroupNames '.eps'],'epsc');
  3737. close all
  3738. end
  3739. end
  3740. %
  3741. % cl

S4_ALLSubjects_NoCut_NoNotch_2_100hz.m at commit 8afead0, no license · at the source

Overview

Authors: Matthew K. Attokaren1, Lu Zhang1,2, Sindhura Mettupalli1, Annabelle C. Singer1
  1. Coulter Department of Biomedical Engineering, Georgia Institute of Technology & Emory University, Atlanta, GA, United States
  2. National Institute of Mental Health, National Institutes of Health, Bethesda, MD, United States
Journal: Imaging neuroscience (Cambridge, Mass.), volume 4, article IMAG.a.1229
Dates: received 26 May 2025; accepted 8 April 2026; published online 13 May 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1162/imag.a.1229 · PMID 42146314 · PMCID PMC13175507 · OpenAlex W4413651284
Open access: diamond, a free copy (OpenAlex)
Status: code verified
Categories: EEG (modality), human (organism)
Methods: Spectral & time-frequency, Statistics, Smoothing, state filtering, decompositions, Connectivity, Preprocessing, Physiology & signal measures
Keywords: gamma sensory stimulation, vigilance, attention, delta oscillations, alpha oscillations, functional connectivity
Topic: Neural dynamics and brain function (Cognitive Neuroscience, Neuroscience), according to OpenAlex
Funding: David and Lucile Packard Foundation (David & Lucile Packard Foundation); NINDS (R01 NS109226, 2RF1NS109226); National Institute on Aging (RF1AG078736); McCamish Foundation; Friends and Alumni of Georgia Tech
Citations: not cited yet (Europe PMC); 60 references in the paper

Abstract

Gamma oscillations (30–100 Hz) have long been theorized to play a key role in sensory processing and attention by coordinating neural firing across distributed neurons. Gamma oscillations can be generated internally by neural circuits during attention or exogenously by stimuli that turn on and off at gamma frequencies. However, it remains unknown if driving gamma activity via exogenous sensory stimulation affects attention. We tested the hypothesis that non-invasive audiovisual stimulation in the form of flashing lights and sounds (flicker) at 40 Hz improves attention in an attentional vigilance task and affects neural oscillations associated with attention. We recorded scalp EEG activity of healthy adults (n = 62) during 1 hour of either 40 Hz audiovisual flicker, no flicker as control, or randomized flicker as sham stimulation, while subjects performed a psychomotor vigilance task. Participants exposed to 40 Hz flicker stimulation had better accuracy and faster reaction times than participants in the control groups. The 40 Hz group showed increased 40 Hz activity compared to the control groups in agreement with previous studies. Surprisingly, 40 Hz subjects had significantly lower delta power (2–4 Hz), which is associated with arousal, and higher functional connectivity in lower alpha (8–10 Hz), which is associated with attention processes. Furthermore, decreased delta power and increased lower alpha functional connectivity were correlated with better attention task performance. This study reveals how 40 Hz audiovisual stimulation improves attention performance with potential implications for therapeutic interventions for attention disorders and attention improvement.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 13 matches between paragraphs and lines of code.

singerlabgt/FlickerEEGAttention_HealthyAdults_Manuscript

License: none: the authors keep all their rights
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: 8afead01e123de2e52cdaf8a3ce5694f63251062, 20 February 2026
Languages: MATLAB (851), R (11), Python (6)
Size: 1,025 files, 868 scripts
Software Heritage: not archived
Found in: “Data and Code Availability”
Holds: README, environment (Behavior Analysis/requirements.txt), tests
Not found: license file, CITATION.cff, continuous integration, documentation
Tools: Statistics and Machine Learning Toolbox (158 files), Signal Processing Toolbox (64 files), EEGLAB (27 files), SPM (14 files), CircStat (10 files), fdr_bh (Benjamini-Hochberg FDR) (9 files), Wavelet Toolbox (9 files), Image Processing Toolbox (8 files), Parallel Computing Toolbox (5 files), NumPy (4 files), pandas (4 files), FieldTrip (3 files), SciPy (3 files), Chronux (2 files), ERPLAB (2 files), ICLabel (2 files), Matplotlib (2 files), seaborn (2 files), Violinplot-Matlab (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
869 files

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 868 scripts, each with its path and the digest of its content;
  • 13 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data and Code Availability

Data from this study is available on OpenNeuro https://doi.org/10.18112/openneuro.ds006222.v1.0.0. Code is available on GitHub: https://github.com/singerlabgt/FlickerEEGAttention_HealthyAdults_Manuscript.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 4 authors, 6 keywords, 5 funders, 60 references.

Cite

This paper

Attokaren, M. K., Zhang, L., Mettupalli, S., & Singer, A. C. (2026). 40 Hz audiovisual stimulation improves sustained attention and related brain oscillations. Imaging neuroscience (Cambridge, Mass.), 4, IMAG.a.1229. https://doi.org/10.1162/imag.a.1229

BibTeX

@article{attokaren202640,
author = {Attokaren, Matthew K. and Zhang, Lu and Mettupalli, Sindhura and Singer, Annabelle C.},
title = {{40 Hz audiovisual stimulation improves sustained attention and related brain oscillations}},
journal = {Imaging neuroscience (Cambridge, Mass.)},
year = {2026},
month = may,
volume = {4},
pages = {IMAG.a.1229},
publisher = {MIT Press},
issn = {2837-6056},
doi = {10.1162/imag.a.1229},
url = {https://doi.org/10.1162/imag.a.1229},
pmid = {42146314},
pmcid = {PMC13175507}
}

RIS

TY - JOUR
AU - Attokaren, Matthew K.
AU - Zhang, Lu
AU - Mettupalli, Sindhura
AU - Singer, Annabelle C.
TI - 40 Hz audiovisual stimulation improves sustained attention and related brain oscillations
T2 - Imaging neuroscience (Cambridge, Mass.)
J2 - Imaging Neurosci (Camb)
PY - 2026
DA - 2026/05/13
VL - 4
SP - IMAG.a.1229
SN - 2837-6056
PB - MIT Press
DO - 10.1162/imag.a.1229
UR - https://doi.org/10.1162/imag.a.1229
LA - en
ER -

CSL-JSON

{
"id": "10.1162/imag.a.1229",
"type": "article-journal",
"title": "40 Hz audiovisual stimulation improves sustained attention and related brain oscillations",
"container-title": "Imaging neuroscience (Cambridge, Mass.)",
"author": [
{
"family": "Attokaren",
"given": "Matthew K."
},
{
"family": "Zhang",
"given": "Lu"
},
{
"family": "Mettupalli",
"given": "Sindhura"
},
{
"family": "Singer",
"given": "Annabelle C."
}
],
"container-title-short": "Imaging Neurosci (Camb)",
"volume": "4",
"page": "IMAG.a.1229",
"DOI": "10.1162/imag.a.1229",
"PMID": "42146314",
"PMCID": "PMC13175507",
"ISSN": "2837-6056",
"publisher": "MIT Press",
"URL": "https://doi.org/10.1162/imag.a.1229",
"language": "en",
"issued": {
"date-parts": [
[
2026,
5,
13
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41598-026-49900-6 [code]
Global neural oscillations underlie performance variability and attentional state fluctuations in humans.
Journal: Scientific reports
In common: ERPLAB, Chronux, Wavelet Toolbox, 7 other tools, 1 reference
[2] doi:10.1038/s41593-026-02357-2 [code]
Experience reorganizes content-specific memory traces in macaques.
Journal: Nature neuroscience
In common: Chronux, fdr_bh (Benjamini-Hochberg FDR), CircStat, 11 other tools
[3] doi:10.64898/2026.03.12.710517 [code]
Cortical excitability inversely modulates fMRI connectivity via low-frequency neuronal coupling
Journal: bioRxiv (preprint)
In common: ERPLAB, Chronux, CircStat, 10 other tools, 1 reference
[4] doi:10.1038/s41467-026-71151-2 [code]
Common and distinct neural correlates of social interaction processing and theory of mind in narratives.
Journal: Nature communications
In common: Violinplot-Matlab, fdr_bh (Benjamini-Hochberg FDR), FieldTrip, 10 other tools
[5] doi:10.1038/s41467-026-75347-4 [code]
Sleep reveals dynamics integrating and segregating movement and stimulus representations in V1.
Journal: Nature communications
In common: Chronux, CircStat, EEGLAB, 9 other tools
[6] doi:10.1016/j.celrep.2026.117646 [code]
Medial entorhinal-hippocampal desynchronization parallels the emergence of memory impairment in a mouse model of Alzheimer's disease pathology.
Journal: Cell reports
In common: Chronux, CircStat, EEGLAB, 9 other tools
[7] doi:10.1038/s41467-026-73106-z [code]
Respiratory pauses highlight sleep architecture in mice.
Journal: Nature communications
In common: Chronux, CircStat, EEGLAB, 8 other tools, EEG
[8] doi:10.1016/j.neuron.2026.03.034 [code]
Dentate gyrus interneurons modulate winner-take-all network dynamics in freely behaving mice.
Journal: Neuron
In common: Chronux, CircStat, EEGLAB, 9 other tools
[9] doi:10.1038/s41467-026-74565-0 [code]
The functional neurobiology of dispositions towards negative emotions.
Journal: Nature communications
In common: Violinplot-Matlab, fdr_bh (Benjamini-Hochberg FDR), FieldTrip, 8 other tools
[10] doi:10.1002/hbm.70577 [code]
Disgust Propensity, Not Disgust Sensitivity, Shapes the Reactivity of a Subjective Disgust Circuit in Humans.
Journal: Human brain mapping
In common: Violinplot-Matlab, fdr_bh (Benjamini-Hochberg FDR), FieldTrip, 8 other tools

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.