diff --git a/Figures/ERPs_phase_reset_nm_alpha.asv b/Figures/ERPs_phase_reset_nm_alpha.asv
deleted file mode 100644
index dcaac0c03abcd53767942426e6f79d1d82ad9eb4..0000000000000000000000000000000000000000
--- a/Figures/ERPs_phase_reset_nm_alpha.asv
+++ /dev/null
@@ -1,434 +0,0 @@
-clear all;
-close all;
-
-addpath(genpath('/users/nemo/software/eeglab'));
-addpath(genpath('/users/nemo/projects/RSN'));
-addpath(genpath('/users/nemo/software/Henry/useful_functions'));
-% addpath(genpath('/user/HS301/m17462/matlab/colorGradient'));
-
-% Savefolder = '/vol/research/nemo/datasets/RSN/data/analysis/Figures/';
-
-% incl_sub = setdiff(1:19,[12]); % 14 excluded because no phasic trials, 9 because no wake eve trials
-incl_sub = setdiff(1:19,[12]); % 14 excluded because no phasic trials, 9 because no wake eve trials
-
-
-%% ERP phase bins - REM
-
-% load('/parallel_scratch/nemo/RSN/analysis/analysis/erp_allsub/ERP_nm_allsub_REM_mICA_avref09-Feb-2024.mat');
-load('/parallel_scratch/nemo/RSN/analysis/analysis/erp_allsub/ERP_nm_broadband_allsub_REM_mICA_avref24-May-2024');
-
-%% plot alpha ERP - averaged across first 10 triggers
-
-t = (-dt*fs:6*dt*fs-1)/fs*1000; % time in ms
-
-colors = linspecer(4);
-
-fig = figure('Renderer','painters','units','normalized','outerposition',[0 0 1 1])
-
-% tileplot = tiledlayout(1,4);
-% 
-% nexttile(1)
-% for trig = 1:10
-for s = 1:size(
-last_trig = ERP_nm_all.trial_data(s,1,:,1)
-trig = 1:40;
-
-% subplot(2,5,trig)
-for con = 1:4
-plot(t,squeeze(nanmean(nanmean(ERP_nm_all.trial_data(incl_sub,con,trig,:),1),3)),'Color',colors(con,:),'LineWidth',2)
-hold on
-end
-
-% hold on
-% plot(t,squeeze(nanmean(nanmean(nanmean(ERP_nm_all.trial_data(incl_sub,1:4,:,:),3),2),1)),'Color','k','LineWidth',3)
-
-xline(0,'LineStyle','--','LineWidth',2);
-
-xlim([-200 1000]);
-xlabel('Time (ms)')
-ylabel('Amplitude (\muV)')
-set(gca,'Fontsize',15,'TickDir','out','LineWidth',2);
-box off
-% axis square
-% legend({'Peak' 'Falling' 'Trough' 'Rising'});
-legend off
-ylim([-3 3]);
-xticks(-200:200:1000);
-title(['Stimulus ',num2str(trig)]);
-
-% end
-
-
-%% plot alpha ERP - averaged across first 10 triggers - alphafilt
-
-t = (-dt*fs:dt*fs-1)/fs*1000; % time in ms
-
-colors = linspecer(4);
-
-fig = figure('Renderer','painters','units','normalized','outerposition',[0 0 1 1])
-
-% tileplot = tiledlayout(1,4);
-% 
-% nexttile(1)
-% for trig = 1:10
-trig = 1:5;
-% subplot(2,5,trig)
-for con = 1:4
-plot(t,squeeze(nanmean(nanmean(ERP_nm_all.trial_data_alphafilt_real(incl_sub,con,trig,:),1),3)),'Color',colors(con,:),'LineWidth',2)
-hold on
-end
-
-% hold on
-% plot(t,squeeze(nanmean(nanmean(nanmean(ERP_nm_all.trial_data_alphafilt_real(incl_sub,1:4,:,:),3),2),1)),'Color','k','LineWidth',3)
-
-xline(0,'LineStyle','--','LineWidth',2);
-
-xlim([-200 1000]);
-xlabel('Time (ms)')
-ylabel('Amplitude (\muV)')
-set(gca,'Fontsize',15,'TickDir','out','LineWidth',2);
-box off
-axis square
-% legend({'Peak' 'Falling' 'Trough' 'Rising'});
-legend off
-ylim([-3 3]);
-xticks(-200:200:1000);
-title(['Stimulus ',num2str(trig)]);
-
-% end
-
-
-%% check non-uniformity across circle - alphafilt
-
-fig = figure('Renderer','painters','units','normalized','outerposition',[0 0 1 1])
-
-% for trig = 1:10
-trig = 1:5;
-% 
-% for samp = 1:size(ERP_nm_all.trial_data_alphafilt_phase,4)
-%         
-%     mean_phase_allcon_REM = reshape(ERP_nm_all.trial_data_alphafilt_phase(incl_sub,1:4,trig,samp),[length(incl_sub)*4 1]);
-%     resultant_phase_allcon_REM = reshape(ERP_nm_all.trial_data_alphafilt_phase_r(incl_sub,1:4,trig,samp),[length(incl_sub)*4 1]);
-%     
-%     mean_phase_allcon_REM(isnan(mean_phase_allcon_REM))= [];
-%     resultant_phase_allcon_REM(resultant_phase_allcon_REM == 0)= [];
-%     
-%     [p_val_REM(samp) z_val_REM(samp)] = circ_rtest(mean_phase_allcon_REM,resultant_phase_allcon_REM);
-%     
-% end
-
-% plot alpha ERP phase
-
-t = (-dt*fs:dt*fs-1)/fs*1000; % time in ms
-
-colors = linspecer(4);
-
-
-% tileplot = tiledlayout(1,4);
-% 
-% nexttile(1)
-
-% subplot(2,5,trig)    
-    
-for con = 1:4
-% plot(t,squeeze(nanmean(ERP_nm_all.trial_data(:,con,trig,:),1)),'Color',colors(con,:))
-incl_sub2 = find(~isnan(squeeze(ERP_nm_all.trial_data_alphafilt_phase(:,con,trig,1))) == 1);
-incl_sub3 = intersect(incl_sub,incl_sub2);
-plot(t,squeeze(circ_mean(ERP_nm_all.trial_data_alphafilt_phase(incl_sub3,con,trig,:))),'Color',colors(con,:))
-hold on
-% sig_samps = find(p_val_REM < 0.05);
-% plot(t(sig_samps),ones(length(sig_samps),1)*3.5,'*','Color','k');
-% hold on
-end
-
-hold on
-% plot(t,squeeze(circ_mean(nanmean(ERP_nm_all.trial_data_alphafilt_phase(incl_sub,1:4,trig,:),2))),'Color','k','LineWidth',3)
-xline(0,'LineStyle','--','LineWidth',2);
-
-xlim([-500 500]);
-xlabel('Time (ms)')
-ylabel('Phase (radians)')
-set(gca,'Fontsize',15,'TickDir','out','LineWidth',2);
-box off
-axis square
-% legend({'Peak' 'Falling' 'Trough' 'Rising'});
-legend off
-ylim([-4 4]);
-yticks([-pi:pi:pi]);
-yticklabels({'-\pi' '0' '\pi'})
-xticks(-200:200:1000);
-title(['Stimulus ',num2str(trig)]);
-
-clear p_val_REM z_val_REM
-% end
-
-% saveas(fig,[Savefolder,'Figure6_ERP_alpha_nm_alphafilt_phase.svg']);
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-%% plot alpha ERP
-
-t = (-dt*fs:dt*fs-1)/fs*1000; % time in ms
-
-colors = linspecer(4);
-
-fig = figure('Renderer','painters','units','normalized','outerposition',[0 0 1 1])
-
-% tileplot = tiledlayout(1,4);
-% 
-% nexttile(1)
-for trig = 1:10
-
-subplot(2,5,trig)
-for con = 1:4
-plot(t,squeeze(nanmean(ERP_nm_all.trial_data(incl_sub,con,trig,:),1)),'Color',colors(con,:),'LineWidth',2)
-hold on
-end
-
-hold on
-plot(t,squeeze(nanmean(nanmean(ERP_nm_all.trial_data(incl_sub,1:4,trig,:),2),1)),'Color','k','LineWidth',3)
-
-xline(0,'LineStyle','--','LineWidth',2);
-
-xlim([-500 500]);
-xlabel('Time (ms)')
-ylabel('Amplitude (\muV)')
-set(gca,'Fontsize',15,'TickDir','out','LineWidth',2);
-box off
-axis square
-% legend({'Peak' 'Falling' 'Trough' 'Rising'});
-legend off
-ylim([-3 3]);
-xticks(-200:200:1000);
-title(['Stimulus ',num2str(trig)]);
-
-end
-
-% saveas(fig,[Savefolder,'Figure6_ERP_alpha_nm_broadband.svg']);
-
-%% plot alpha ERP - alpha filtered
-
-t = (-dt*fs:dt*fs-1)/fs*1000; % time in ms
-
-colors = linspecer(4);
-
-fig = figure('Renderer','painters','units','normalized','outerposition',[0 0 1 1])
-
-% tileplot = tiledlayout(1,4);
-% 
-% nexttile(1)
-for trig = 1:10
-
-subplot(2,5,trig)
-
-for con = 1:4
-plot(t,squeeze(nanmean(ERP_nm_all.trial_data_alphafilt_real(incl_sub,con,trig,:),1)),'Color',colors(con,:))
-hold on
-end
-
-hold on
-plot(t,squeeze(nanmean(nanmean(ERP_nm_all.trial_data_alphafilt_real(incl_sub,1:4,trig,:),2),1)),'Color','k','LineWidth',3)
-
-xline(0,'LineStyle','--','LineWidth',2);
-
-xlim([-500 500]);
-xlabel('Time (ms)')
-ylabel('Amplitude (\muV)')
-set(gca,'Fontsize',15,'TickDir','out','LineWidth',2);
-box off
-axis square
-% legend({'Peak' 'Falling' 'Trough' 'Rising'});
-legend off
-ylim([-2 2]);
-xticks(-200:200:1000);
-title(['Stimulus ',num2str(trig)]);
-
-end
-
-% saveas(fig,[Savefolder,'Figure6_ERP_alpha_nm_alphafilt.svg']);
-
-
-%% plot alpha ERP - theta filtered
-
-t = (-dt*fs:dt*fs-1)/fs*1000; % time in ms
-
-colors = linspecer(4);
-
-fig = figure('Renderer','painters','units','normalized','outerposition',[0 0 1 1])
-
-% tileplot = tiledlayout(1,4);
-% 
-% nexttile(1)
-for trig = 1:10
-
-subplot(2,5,trig)
-
-for con = 1:4
-plot(t,squeeze(nanmean(ERP_nm_all.trial_data_thetafilt_real(incl_sub,con,trig,:),1)),'Color',colors(con,:))
-hold on
-end
-
-hold on
-plot(t,squeeze(nanmean(nanmean(ERP_nm_all.trial_data_thetafilt_real(incl_sub,1:4,trig,:),2),1)),'Color','k','LineWidth',3)
-
-xline(0,'LineStyle','--','LineWidth',2);
-
-xlim([-500 500]);
-xlabel('Time (ms)')
-ylabel('Amplitude (\muV)')
-set(gca,'Fontsize',15,'TickDir','out','LineWidth',2);
-box off
-axis square
-% legend({'Peak' 'Falling' 'Trough' 'Rising'});
-legend off
-ylim([-2 2]);
-xticks(-200:200:1000);
-title(['Stimulus ',num2str(trig)]);
-
-end
-
-% saveas(fig,[Savefolder,'Figure6_ERP_alpha_nm_thetafilt.svg']);
-
-
-%% check non-uniformity across circle - alphafilt
-
-fig = figure('Renderer','painters','units','normalized','outerposition',[0 0 1 1])
-
-for trig = 1:10
-
-for samp = 1:size(ERP_nm_all.trial_data_alphafilt_phase,4)
-        
-    mean_phase_allcon_REM = reshape(ERP_nm_all.trial_data_alphafilt_phase(incl_sub,1:4,trig,samp),[length(incl_sub)*4 1]);
-    resultant_phase_allcon_REM = reshape(ERP_nm_all.trial_data_alphafilt_phase_r(incl_sub,1:4,trig,samp),[length(incl_sub)*4 1]);
-    
-    mean_phase_allcon_REM(isnan(mean_phase_allcon_REM))= [];
-    resultant_phase_allcon_REM(resultant_phase_allcon_REM == 0)= [];
-    
-    [p_val_REM(samp) z_val_REM(samp)] = circ_rtest(mean_phase_allcon_REM,resultant_phase_allcon_REM);
-    
-end
-
-% plot alpha ERP phase
-
-t = (-dt*fs:dt*fs-1)/fs*1000; % time in ms
-
-colors = linspecer(4);
-
-
-% tileplot = tiledlayout(1,4);
-% 
-% nexttile(1)
-
-subplot(2,5,trig)    
-    
-for con = 1:4
-% plot(t,squeeze(nanmean(ERP_nm_all.trial_data(:,con,trig,:),1)),'Color',colors(con,:))
-incl_sub2 = find(~isnan(squeeze(ERP_nm_all.trial_data_alphafilt_phase(:,con,trig,1))) == 1);
-incl_sub3 = intersect(incl_sub,incl_sub2);
-plot(t,squeeze(circ_mean(ERP_nm_all.trial_data_alphafilt_phase(incl_sub3,con,trig,:))),'Color',colors(con,:))
-hold on
-sig_samps = find(p_val_REM < 0.05);
-plot(t(sig_samps),ones(length(sig_samps),1)*3.5,'*','Color','k');
-hold on
-end
-
-hold on
-% plot(t,squeeze(circ_mean(nanmean(ERP_nm_all.trial_data_alphafilt_phase(incl_sub,1:4,trig,:),2))),'Color','k','LineWidth',3)
-xline(0,'LineStyle','--','LineWidth',2);
-
-xlim([-500 500]);
-xlabel('Time (ms)')
-ylabel('Phase (radians)')
-set(gca,'Fontsize',15,'TickDir','out','LineWidth',2);
-box off
-axis square
-% legend({'Peak' 'Falling' 'Trough' 'Rising'});
-legend off
-ylim([-4 4]);
-yticks([-pi:pi:pi]);
-yticklabels({'-\pi' '0' '\pi'})
-xticks(-200:200:1000);
-title(['Stimulus ',num2str(trig)]);
-
-clear p_val_REM z_val_REM
-end
-
-% saveas(fig,[Savefolder,'Figure6_ERP_alpha_nm_alphafilt_phase.svg']);
-
-
-%% check non-uniformity across circle - thetafilt
-
-fig = figure('Renderer','painters','units','normalized','outerposition',[0 0 1 1])
-
-for trig = 1:10
-
-for samp = 1:size(ERP_nm_all.trial_data_alphafilt_phase,4)
-        
-    mean_phase_allcon_REM = reshape(ERP_nm_all.trial_data_thetafilt_phase(incl_sub,1:4,trig,samp),[length(incl_sub)*4 1]);
-    resultant_phase_allcon_REM = reshape(ERP_nm_all.trial_data_thetafilt_phase_r(incl_sub,1:4,trig,samp),[length(incl_sub)*4 1]);
-    
-    mean_phase_allcon_REM(isnan(mean_phase_allcon_REM))= [];
-    resultant_phase_allcon_REM(resultant_phase_allcon_REM == 0)= [];
-    
-    [p_val_REM(samp) z_val_REM(samp)] = circ_rtest(mean_phase_allcon_REM,resultant_phase_allcon_REM);
-    
-end
-
-% plot alpha ERP phase
-
-t = (-dt*fs:dt*fs-1)/fs*1000; % time in ms
-
-colors = linspecer(4);
-
-
-% tileplot = tiledlayout(1,4);
-% 
-% nexttile(1)
-
-subplot(2,5,trig)    
-    
-for con = 1:4
-% plot(t,squeeze(nanmean(ERP_nm_all.trial_data(:,con,trig,:),1)),'Color',colors(con,:))
-incl_sub2 = find(~isnan(squeeze(ERP_nm_all.trial_data_thetafilt_phase(:,con,trig,1))) == 1);
-incl_sub3 = intersect(incl_sub,incl_sub2);
-plot(t,squeeze(circ_mean(ERP_nm_all.trial_data_thetafilt_phase(incl_sub3,con,trig,:))),'Color',colors(con,:))
-hold on
-sig_samps = find(p_val_REM < 0.05);
-plot(t(sig_samps),ones(length(sig_samps),1)*3.5,'*','Color','k');
-hold on
-end
-
-hold on
-% plot(t,squeeze(circ_mean(nanmean(ERP_nm_all.trial_data_alphafilt_phase(incl_sub,1:4,trig,:),2))),'Color','k','LineWidth',3)
-xline(0,'LineStyle','--','LineWidth',2);
-
-xlim([-500 500]);
-xlabel('Time (ms)')
-ylabel('Phase (radians)')
-set(gca,'Fontsize',15,'TickDir','out','LineWidth',2);
-box off
-axis square
-% legend({'Peak' 'Falling' 'Trough' 'Rising'});
-legend off
-ylim([-4 4]);
-yticks([-pi:pi:pi]);
-yticklabels({'-\pi' '0' '\pi'})
-xticks(-200:200:1000);
-title(['Stimulus ',num2str(trig)]);
-
-clear p_val_REM z_val_REM
-end
-
-saveas(fig,[Savefolder,'Figure6_ERP_alpha_nm_thetafilt_phase.svg']);
-
diff --git a/Figures/Figure3C_3H_eBOSC_frequency_corr_nwaves.asv b/Figures/Figure3C_3H_eBOSC_frequency_corr_nwaves.asv
deleted file mode 100644
index bc4aa13901263aa7a56e8819b0d1d12bf75aaee3..0000000000000000000000000000000000000000
--- a/Figures/Figure3C_3H_eBOSC_frequency_corr_nwaves.asv
+++ /dev/null
@@ -1,160 +0,0 @@
-clear all;
-close all;
-
-%%
-
-load('/parallel_scratch/nemo/RSN/analysis/analysis/oscillation_detection/Oscillation_Freq_ISI_allsub_d03_08-Feb-2024.mat');
-close
-
-Savefolder = '/parallel_scratch/nemo/RSN/analysis/analysis/Figures/';
-
-%%
-incl_sub = setdiff(1:19,[12]);
-
-nonan_sub_alpha = intersect(find(~isnan(alpha_peak_findpeaks)==1),incl_sub);
-nonan_sub_alphastim_alpha = intersect(find(~isnan(alphastim_alpha_peak_findpeaks)==1),incl_sub);
-nonan_sub_alpha = intersect(nonan_sub_alpha,nonan_sub_alphastim_alpha);
-length(nonan_sub_alpha)
-
-[r p] = corr(alpha_peak_findpeaks(nonan_sub_alpha)',alphastim_alpha_peak_findpeaks(nonan_sub_alpha)')
-
-fig = figure('Renderer','painters','units','normalized','outerposition',[0 0 1 1])
-
-mdl_alpha = fitlm(alpha_peak_findpeaks(nonan_sub_alpha)',alphastim_alpha_peak_findpeaks(nonan_sub_alpha)');
-pl = plot(mdl_alpha);
-pl(2,1).Color = [0.5 0.5 0.5];
-pl(3,1).Color = [0.5 0.5 0.5];
-pl(4,1).Color = [0.5 0.5 0.5];
-pl(2,1).LineWidth = 3;
-pl(3,1).LineWidth = 3;
-pl(4,1).LineWidth = 3;
-
-hold on
-plot(alpha_peak_findpeaks(nonan_sub_alpha),alphastim_alpha_peak_findpeaks(nonan_sub_alpha),'o','MarkerFaceColor','w','MarkerEdgeColor','k','MarkerSize',15,'LineWidth',3);  
-p_reg_alpha = mdl_alpha.Coefficients(2,4)
-r_squared = mdl_alpha.Rsquared.Ordinary
-intercept_alpha = mdl_alpha.Coefficients(1,1);
-slope_alpha = mdl_alpha.Coefficients(2,1);
-slope_SE_alpha = mdl_alpha.Coefficients(2,2)
-
-%    fig = figure('Renderer','painters','units','normalized','outerposition',[0 0 1 1])
-%     scatter(mch_allcon_ifq_alpha_off(incl_sub),nanmean(m_ISI_theta(incl_sub,5:8),2),100,'k','LineWidth',3);
-%     scatter(alpha_peak_findpeaks(nonan_sub_alpha),alphastim_alpha_peak_findpeaks(nonan_sub_alpha),300,'MarkerEdgeColor','k','LineWidth',3);
-%     l = lsline;
-%     l.LineWidth = 3;
-%     l.Color = 'r';
-    xlabel('REM Alpha Freq (Hz)');
-    ylabel('Alpha ISI Freq (Hz)');
-%     title(['OFF, r = ', num2str(r), ', p = ',num2str(p)]);
-    set(gca,'Fontsize',35,'TickDir','out','LineWidth',3);
-    title('');    
-    box off
-    axis square
-    xlim([7 11]);
-    ylim([7 11]);
-    xticks(7.0:1.0:11.0);
-    xticklabels(7.0:1.0:11.0);
-    xtickformat('%.1f')
-    yticks(7.0:1.0:11.0);
-    yticklabels(7.0:1.0:11.0);
-    ytickformat('%.1f')
-    legend off
-
-% saveas(fig,[Savefolder,'Figure3C_corr_REM_alpha_oscillations_ISI_Alpha_',num2str(lower_freq_alpha),'_',num2str(higher_freq_alpha),'.svg']);
-
-%%
-
-
-incl_sub = setdiff(1:19,[12]);
-
-nonan_sub_theta = intersect(find(~isnan(theta_peak_findpeaks)==1),incl_sub);
-nonan_sub_thetastim_theta = intersect(find(~isnan(thetastim_theta_peak_findpeaks)==1),incl_sub);
-nonan_sub_theta = intersect(nonan_sub_theta,nonan_sub_thetastim_theta);
-length(nonan_sub_theta)
-
-[r p] = corr(theta_peak_findpeaks(nonan_sub_theta)',thetastim_theta_peak_findpeaks(nonan_sub_theta)')
-
-fig = figure('Renderer','painters','units','normalized','outerposition',[0 0 1 1])
-
-mdl_theta = fitlm(theta_peak_findpeaks(nonan_sub_theta)',thetastim_theta_peak_findpeaks(nonan_sub_theta)');
-pl = plot(mdl_theta);
-pl(2,1).Color = [0.5 0.5 0.5];
-pl(3,1).Color = [0.5 0.5 0.5];
-pl(4,1).Color = [0.5 0.5 0.5];
-pl(2,1).LineWidth = 3;
-pl(3,1).LineWidth = 3;
-pl(4,1).LineWidth = 3;
-
-hold on
-plot(theta_peak_findpeaks(nonan_sub_theta),thetastim_theta_peak_findpeaks(nonan_sub_theta),'o','MarkerFaceColor','w','MarkerEdgeColor','k','MarkerSize',15,'LineWidth',3);  
-p_reg_theta = mdl_theta.Coefficients(2,4)
-r_squared_theta = mdl_theta.Rsquared.Ordinary
-intercept_theta = mdl_theta.Coefficients(1,1);
-slope_theta = mdl_theta.Coefficients(2,1);
-slope_SE_alpha = mdl_alpha.Coefficients(2,2)
-
-
-%    fig = figure('Renderer','painters','units','normalized','outerposition',[0 0 1 1])
-%     scatter(mch_allcon_ifq_alpha_off(incl_sub),nanmean(m_ISI_theta(incl_sub,5:8),2),100,'k','LineWidth',3);
-%     scatter(theta_peak_findpeaks(nonan_sub_theta)',thetastim_theta_peak_findpeaks(nonan_sub_theta)',300,'MarkerEdgeColor','k','LineWidth',3);
-%     l = lsline;
-%     l.LineWidth = 3;
-%     l.Color = 'k';
-    xlabel('REM Theta Freq (Hz)');
-    ylabel('Theta ISI Freq (Hz)');
-%     title(['OFF, r = ', num2str(r), ', p = ',num2str(p)]);
-    set(gca,'Fontsize',35,'TickDir','out','LineWidth',3);
-    title('');    
-    box off
-    axis square
-    xlim([4.5 8]);
-    ylim([4.5 8]);
-    xticks(4.5:1:8);
-    xticklabels(4.5:1:8);
-    xtickformat('%.1f')
-    yticks(4.5:1:8);
-    yticklabels(4.5:1:8);
-    ytickformat('%.1f')
-    legend off
-
-
-    saveas(fig,[Savefolder,'Figure3H_corr_REM_theta_oscillations_ISI_Theta_',num2str(lower_freq_theta),'_',num2str(higher_freq_theta),'.svg']);
-
-%% number of waves
-
-m_n_rem_alpha = round(nanmean(n_rem_alpha(incl_sub)))
-% sem_n_rem_alpha = round(nanstd(n_rem_alpha(incl_sub))/sqrt(length(incl_sub)))
-sd_n_rem_alpha = round(nanstd(n_rem_alpha(incl_sub)))
-
-m_n_rem_theta = round(nanmean(n_rem_theta(incl_sub)))
-% sem_n_rem_theta = round(nanstd(n_rem_theta(incl_sub))/sqrt(length(incl_sub)))
-sd_n_rem_theta = round(nanstd(n_rem_theta(incl_sub)))
-
-
-m_abundance_rem_alpha = nanmean(abundance_rem_alpha(incl_sub))
-% sem_abundance_rem_alpha = nanstd(abundance_rem_alpha(incl_sub))/sqrt(length(incl_sub))
-sd_abundance_rem_alpha = nanstd(abundance_rem_alpha(incl_sub))
-
-m_abundance_rem_theta = nanmean(abundance_rem_theta(incl_sub))
-% sem_abundance_rem_theta = nanstd(abundance_rem_theta(incl_sub))/sqrt(length(incl_sub))
-sd_abundance_rem_theta = nanstd(abundance_rem_theta(incl_sub))
-
-
-% m_density_rem_alpha = nanmean(density_rem_alpha(incl_sub))
-% sem_density_rem_alpha = nanstd(density_rem_alpha(incl_sub))/sqrt(length(incl_sub))
-% 
-% m_density_rem_theta = nanmean(density_rem_theta(incl_sub))
-% sem_density_rem_theta = nanstd(density_rem_theta(incl_sub))/sqrt(length(incl_sub))
-% 
-%% FWHM
-
-m_fwhm_rem_osc = nanmean(fwhmx_osc(incl_sub))
-sd_fwhm_rem_osc = nanstd(fwhmx_osc(incl_sub))
-
-m_fwhm_rem_alphastim = nanmean(fwhmx_alphastim(incl_sub))
-sd_fwhm_rem_alphastim = nanstd(fwhmx_alphastim(incl_sub))
-
-m_fwhm_rem_thetastim = nanmean(fwhmx_thetastim(incl_sub))
-sd_fwhm_rem_thetastim = nanstd(fwhmx_thetastim(incl_sub))
-
-
diff --git a/Figures/not_used/ERPs_phase_response_alpha.m b/Figures/not_used/ERPs_phase_response_alpha.m
deleted file mode 100644
index 8ac90c1adcb63373c854ebd0094e0ea26a57bc8e..0000000000000000000000000000000000000000
--- a/Figures/not_used/ERPs_phase_response_alpha.m
+++ /dev/null
@@ -1,1002 +0,0 @@
-clear all;
-close all;
-
-addpath(genpath('/user/HS301/m17462/matlab/eeglab'));
-addpath(genpath('/user/HS301/m17462/matlab/Scripts/RSN'));
-addpath(genpath('/user/HS301/m17462/matlab/Henry/useful_functions'));
-addpath(genpath('/user/HS301/m17462/matlab/colorGradient'));
-
-Savefolder = '/vol/research/nemo/datasets/RSN/data/analysis/Figures/';
-
-incl_sub = setdiff(1:19,12);
-
-Folderpath = '/vol/research/nemo/datasets/RSN/data/hdEEG/';
-sub_Folderpath = dir([Folderpath,'RSN*']);
-
-%% load all trial data - REM
-
-for s = 1:length(sub_Folderpath)
-    
-    erp_good_file = dir([Folderpath,sub_Folderpath(s).name,filesep,'*_sleep*_erp_good.mat']);
-    
-    load([Folderpath,sub_Folderpath(s).name,filesep,erp_good_file(1).name]);
-    
-    ERP_trials_allsub_REM{s} = ERP; 
-    
-end
-
-
-%% load all trial data - wake eve
-
-for s = 1:length(sub_Folderpath)
-    
-    erp_good_file = dir([Folderpath,sub_Folderpath(s).name,filesep,'*_wake_e*_erp_good.mat']);
-    
-    load([Folderpath,sub_Folderpath(s).name,filesep,erp_good_file(1).name]);
-    
-    ERP_trials_allsub_wake_e{s} = ERP; 
-    
-end
-
-%% load all trial data - wake mor
-
-for s = 1:length(sub_Folderpath)
-    
-    erp_good_file = dir([Folderpath,sub_Folderpath(s).name,filesep,'*_wake_m*_erp_good.mat']);
-    
-    load([Folderpath,sub_Folderpath(s).name,filesep,erp_good_file(1).name]);
-    
-    ERP_trials_allsub_wake_m{s} = ERP; 
-    
-end
-
-%% phase response curve - tonic REM
-
-colors = linspecer(4);
-             
-%compute for many onset phases speeding up/slowing down
-      
-state = 1; % 1 = eyes open, 2 = eyes closed
-chan = 2;
-             
-clear angle_start_end phase_after R_after angs clear sem_end pred_angle R_end
-
-trigsamp = 501;
-nbins = 10;
-fs = 500;
-osc_freq = 9;
-
-% angs(1,:) = wrapToPi(deg2rad([0:360/10:360]));
-angs(1,:) = wrapToPi(deg2rad([0:360/nbins:360]));
-angs(2,:) = angs-pi/nbins;
-angs(3,:) = angs(1,:)+pi/nbins;
-             
-fig = figure('Renderer','painters','units','normalized','outerposition',[0 0 .4 1])
-hold on
-% for chan = [12]
-%                  clear phase_after phase_start_end
-G = 1;
-
-cycs_all = [20/2 100/2 150/2 200/2 250/2];
-
-for angle = 1:length(cycs_all) % [10 75 85 95 105];%75 100 125 150]%(25:25:200)
-    cycs = cycs_all(angle);
-    A = 1;
-    for ang = 1:size(angs,2) %0:18:360;%[36:10:360-36]
-             
-         for cond = 1
-                 clear y
-                 for s = 1:length(incl_sub)
-                     
-                     sub = incl_sub(s);
-                     
-%                         tmp = [];
-%                         for condo = 1:4
-%                          tmp = [tmp;group_phase.sub{sub}.chan{chan}.cond{condo}];
-                        tmp = ERP_trials_allsub_REM{sub}.trial_data_alphafilt_phase(ERP_trials_allsub_REM{sub}.tonic_ndx,:); % trials x samps
-%                         tmp = ERP_trials_allsub_REM{sub}.trial_data_alphafilt_phase(ERP_trials_allsub_REM{sub}.goodtrial_ndx,:); % trials x samps
-%                           tmp = ERP_trials_allsub_wake_m{sub}.trial_data_alphafilt_phase{state}(ERP_trials_allsub_wake_m{sub}.good_trial_ndx{state},:);
-%                           tmp = ERP_trials_allsub_wake_e{sub}.trial_data_alphafilt_phase{state}(ERP_trials_allsub_wake_e{sub}.good_trial_ndx{state},:);
-%                         end
-                                              
-                        tmp = tmp(tmp(:,trigsamp)>angs(2,ang)&tmp(:,trigsamp)<angs(3,ang),:);
-                         
-                        if ~isempty(tmp)
-                                      
-                        Mn = circ_mean(tmp); % average across all trials in that bin
-                        R = circ_r(tmp);
-                        to_find = Mn(trigsamp);
-                     
-                        cosed = cos(circ_dist(Mn,to_find));
-                                          
-                        [pks,locs] = findpeaks(cosed);
-                                       
-                     
-                        phase_before(s,cond) =  Mn(trigsamp);
-                        cycle_length(s,cond) = abs(trigsamp-locs(find(locs==trigsamp)-1)); % duration from previous peak to stim peak in samps
-                        phase_after(s,cond,1) = Mn(trigsamp);
-                        R_after(s,cond,1) = R(trigsamp);
-                                        
-                    
-                        D = (diff(locs(locs>=trigsamp-100 & locs<=trigsamp)));
-                        %D = D(D>30);
-                        cycles(s,cond,chan) = mean(D); % mean duration of all cycles between -200 ms and 0
-                    
-                    
-                            for cyc = 1:round(cycs)
-                                   
-                                phase_after(s,cond,cyc+1) = Mn(trigsamp+cyc); % mean phase for all samples from trigsamp to endsamp
-                                phase_change(s,cond,cyc) = circ_dist(phase_after(s,cond,cyc),phase_before(s,cond)); % phase distance from each samp after to trigsamp
-                                R_after(s,cond,cyc+1) = R(trigsamp+cyc);%  resultant for all samples from trigsamp to endsamp      
-                        
-                            end
-                    
-                                                               
-                        else
-                 
-                            phase_after(s,cond,:) = NaN(1,1,cycs+1);
-                            R_after(s,cond,:) = NaN(1,1,cycs+1);
-                        end
-            
-                 end
-           
-         end
-             
-             
-%              if cycs == trigsamp
-%              phase_series(ang,:) = squeeze(circ_mean(phase_after));
-%              end
-%              clear cumul_phase_change
-             
-             cumul_phase_change = unwrap(squeeze(phase_after(:,1,:))')';
-
-             angle_start_end(1,ang) = (wrapToPi(angs(1,ang))); % angle at center of bin
-             tmp = phase_after(:,cond,end); % phase of endsamp
-             tmp2 = R_after(:,cond,end); %  resultant of endsamp
-             tmp2 = tmp2(~isnan(tmp2));
-             tmp = tmp(~isnan(tmp));
-                 
-             angle_start_end(2,ang) = (wrapToPi(circ_mean(tmp))); % mean of endsamp phase across participants
-             yy = num2str(cycs);
-             yy = str2num(yy(end-1:end))/(fs/osc_freq); % calculate how much of a cycle this amount of time is
-             pred_angle(ang) = wrapToPi(angs(1,ang)+ 2*yy*pi);
-             R_end(1,ang) = circ_r(tmp,tmp2); % resultant of endsamp across participants
-                
-             angle_start_end(3,ang) = (circ_dist(circ_mean(tmp),pred_angle(ang))); % phase distance from endsamp to predicted angle
-
-             sem_end(1,ang) = circ_std(tmp)./sqrt(length(incl_sub)); % sem of end phase
-             sem_end(2,ang) = circ_std(circ_dist(angs(1,ang),tmp))./sqrt(length(incl_sub)); % sem of difference between angle at center of bin and end phase
-
-             A = A+1;
-             
-             end
-             
-             [B I] = sort((angle_start_end(2,:)));
-             
-             colorz = linspecer(length(angle_start_end(3,:)));
-
-             subplot(5,4,(angle-1)*4+1)
-             hold on
-             errorbar(angle_start_end(1,:),angle_start_end(2,:),sem_end(1,:),'Color','k','LineStyle','none') % angle at center of bin vs mean phase at end, errorbar
-             scatter(angle_start_end(1,:),pred_angle,2.5,'k','filled') % angle at center of bin vs predicted angle
-             
-             for i = 1:length(angle_start_end(3,:))
-
-                if chan == 2
-
-                    scatter(angle_start_end(1,i),angle_start_end(2,i),25,colorz(i,:),'filled','MarkerEdgeColor','k') % angle at center of bins vs mean at end 
-                else
-                    scatter(angle_start_end(1,:),angle_start_end(2,:),25,colors(2,:))
-                    errorbar(angle_start_end(1,:),angle_start_end(2,:),sem_end(1,:),'Color',colors(2,:),'LineStyle','none') 
-                end
-             end
-             
-%              G = G+1;
-
-             ylim([-pi pi])
-             xlim([-pi pi])
-             title(['Delay: ',num2str(cycs/fs*1000),' ms'])     
-
-             xlabel('Start Phase')
-             ylabel('End Phase')
-             
-             subplot(5,4,(angle-1)*4+2)
-             hold on
-             errorbar(angle_start_end(1,:),angle_start_end(3,:),sem_end(2,:),'Color','k','LineStyle','none') % angle at center of bin vs phase distance from endsamp to predicted angle, errorbar
-             
-             for i = 1:length(angle_start_end(3,:))
-                if chan == 2
-                    scatter(angle_start_end(1,i),angle_start_end(3,i),25,colorz(i,:),'filled','MarkerEdgeColor','k')                         
-                else
-                    scatter(angle_start_end(1,:),angle_start_end(3,:),25,colors(2,:))
-                    errorbar(angle_start_end(1,:),angle_start_end(3,:),sem_end(2,:),'Color',colors(2,:),'LineStyle','none')    
-                end
-             end
-             
-%             G = G+1;
-            a = 1;
-            for A = [0 90 180 270]
-%             for A = [340 60 160 240]
-                xline(wrapToPi(deg2rad(A)),'Color',colors(a,:)) 
-                a = a+1;
-            end
-             %xlim([-180 180])
-             xlabel('Start Phase')
-             ylabel('Phase Change')
-             yline(0)
-%              text(.95*-pi,-.9*pi,'slower')
-             %ylim([-200 200])
-             %ylim([.1 .14])
-             ylim([-pi pi])
-             xlim([-pi pi])
-%              text(.95*-pi,.9*pi,'faster')
-                         
-                          
-%              subplot(5,4,G)
-%              polarplot([0 circ_mean(angle_start_end(2,:)')],[0 (circ_r(angle_start_end(2,:)',R_end'))],'k','LineWidth',2)
-%              hold on
-% 
-%             if chan == 2
-%                 B = 0.5;
-%                 for i = 1:length(angle_start_end(3,:))
-%                     polarscatter(angle_start_end(2,i),R_end(i),35,colorz(i,:),'filled','MarkerEdgeColor','k')
-%                     B = B+.025;
-%                     hold on
-%                     if G == 3
-%                         before(i) = angle_start_end(2,i);
-%                     else
-%   
-%                         hold on
-%                     end
-%                 end
-%                 
-%             else
-%                 polarplot(angle_start_end(1,:),(circ_dist(circ_mean(angle_start_end(3,:)'),angle_start_end(3,:)))+pi,'Color',colors(2,:))
-% 
-%             end
-%              hold on
-% 
-%              rlim([0 1.1])
-%              G = G+1;
-%              
-%              subplot(5,4,G)
-%              polarplot([0 circ_mean(angle_start_end(3,:)')],[0 (circ_r(angle_start_end(3,:)',R_end'))],'k','LineWidth',2)
-%              hold on
-% 
-%              if chan == 2
-%                 B = .5;
-%                 
-%                 for i = 1:length(angle_start_end(3,:))
-%                     polarscatter(angle_start_end(3,i),R_end(i),35,colorz(i,:),'filled','MarkerEdgeColor','k')
-%                     hold on
-%                     B = B+.025;
-%     
-%                     if G == 3
-%                         before(i) = angle_start_end(2,i);
-%                     else
-%                         hold on
-%                     end
-%                 end
-%              end
-%              
-%              hold on
-%              rlim([0 1.1])
-%              G = G+1;
-             
-end
-             
-             
-saveas(fig,[Savefolder,'Figure6_phase_response_curve_freq',num2str(osc_freq),'_alphafilt_tonic_REM.svg']);
-
-%% phase response curve - phasic REM
-
-% colors = linspecer(4);
-%              
-% %compute for many onset phases speeding up/slowing down
-%       
-% state = 1; % 1 = eyes open, 2 = eyes closed
-% chan = 2;
-%              
-% clear angle_start_end phase_after R_after angs clear sem_end pred_angle R_end
-% 
-% trigsamp = 501;
-% nbins = 10;
-% fs = 500;
-% osc_freq = 9;
-% 
-% % angs(1,:) = wrapToPi(deg2rad([0:360/10:360]));
-% angs(1,:) = wrapToPi(deg2rad([0:360/nbins:360]));
-% angs(2,:) = angs-pi/nbins;
-% angs(3,:) = angs(1,:)+pi/nbins;
-%              
-% fig = figure('Renderer','painters','units','normalized','outerposition',[0 0 .4 1])
-% hold on
-% % for chan = [12]
-% %                  clear phase_after phase_start_end
-% G = 1;
-% 
-% for cycs = [20/2 100/2 150/2 200/2 250/2]% [10 75 85 95 105];%75 100 125 150]%(25:25:200)
-%     A = 1;
-%     for ang = 1:size(angs,2) %0:18:360;%[36:10:360-36]
-%              
-%          for cond = 1
-%                  clear y
-%                  for s = 1:length(incl_sub)
-%                      
-%                      sub = incl_sub(s);
-%                      
-% %                         tmp = [];
-% %                         for condo = 1:4
-% %                          tmp = [tmp;group_phase.sub{sub}.chan{chan}.cond{condo}];
-%                         tmp = ERP_trials_allsub_REM{sub}.trial_data_alphafilt_phase(ERP_trials_allsub_REM{sub}.phasic_ndx,:); % trials x samps
-% %                         tmp = ERP_trials_allsub_REM{sub}.trial_data_alphafilt_phase(ERP_trials_allsub_REM{sub}.goodtrial_ndx,:); % trials x samps
-% %                           tmp = ERP_trials_allsub_wake_m{sub}.trial_data_alphafilt_phase{state}(ERP_trials_allsub_wake_m{sub}.good_trial_ndx{state},:);
-% %                           tmp = ERP_trials_allsub_wake_e{sub}.trial_data_alphafilt_phase{state}(ERP_trials_allsub_wake_e{sub}.good_trial_ndx{state},:);
-% %                         end
-%                                               
-%                         tmp = tmp(tmp(:,trigsamp)>angs(2,ang)&tmp(:,trigsamp)<angs(3,ang),:);
-%                          
-%                         if ~isempty(tmp)
-%                                       
-%                         Mn = circ_mean(tmp); % average across all trials in that bin
-%                         R = circ_r(tmp);
-%                         to_find = Mn(trigsamp);
-%                      
-%                         cosed = cos(circ_dist(Mn,to_find));
-%                                           
-%                         [pks,locs] = findpeaks(cosed);
-%                                        
-%                      
-%                         phase_before(s,cond) =  Mn(trigsamp);
-%                         cycle_length(s,cond) = abs(trigsamp-locs(find(locs==trigsamp)-1)); % duration from previous peak to stim peak in samps
-%                         phase_after(s,cond,1) = Mn(trigsamp);
-%                         R_after(s,cond,1) = R(trigsamp);
-%                                         
-%                     
-%                         D = (diff(locs(locs>=trigsamp-100 & locs<=trigsamp)));
-%                         %D = D(D>30);
-%                         cycles(s,cond,chan) = mean(D); % mean duration of all cycles between -200 ms and 0
-%                     
-%                     
-%                             for cyc = 1:round(cycs)
-%                                    
-%                                 phase_after(s,cond,cyc+1) = Mn(trigsamp+cyc); % mean phase for all samples from trigsamp to endsamp
-%                                 phase_change(s,cond,cyc) = circ_dist(phase_after(s,cond,cyc),phase_before(s,cond)); % phase distance from each samp after to trigsamp
-%                                 R_after(s,cond,cyc+1) = R(trigsamp+cyc);%  resultant for all samples from trigsamp to endsamp      
-%                         
-%                             end
-%                     
-%                                                                
-%                         else
-%                  
-%                             phase_after(s,cond,:) = NaN(1,1,cycs+1);
-%                             R_after(s,cond,:) = NaN(1,1,cycs+1);
-%                         end
-%             
-%                  end
-%            
-%          end
-%              
-%              
-% %              if cycs == trigsamp
-% %              phase_series(ang,:) = squeeze(circ_mean(phase_after));
-% %              end
-% %              clear cumul_phase_change
-%              
-%              cumul_phase_change = unwrap(squeeze(phase_after(:,1,:))')';
-% 
-%              angle_start_end(1,ang) = (wrapToPi(angs(1,ang))); % angle at center of bin
-%              tmp = phase_after(:,cond,end); % phase of endsamp
-%              tmp2 = R_after(:,cond,end); %  resultant of endsamp
-%              tmp2 = tmp2(~isnan(tmp2));
-%              tmp = tmp(~isnan(tmp));
-%                  
-%              angle_start_end(2,ang) = (wrapToPi(circ_mean(tmp))); % mean of endsamp phase across participants
-%              yy = num2str(cycs);
-%              yy = str2num(yy(end-1:end))/(fs/osc_freq); % calculate how much of a cycle this amount of time is
-%              pred_angle(ang) = wrapToPi(angs(1,ang)+ 2*yy*pi);
-%              R_end(1,ang) = circ_r(tmp,tmp2); % resultant of endsamp across participants
-%                 
-%              angle_start_end(3,ang) = (circ_dist(circ_mean(tmp),pred_angle(ang))); % phase distance from endsamp to predicted angle
-% 
-%              sem_end(1,ang) = circ_std(tmp)./sqrt(length(incl_sub)); % sem of end phase
-%              sem_end(2,ang) = circ_std(circ_dist(angs(1,ang),tmp))./sqrt(length(incl_sub)); % sem of difference between angle at center of bin and end phase
-% 
-%              A = A+1;
-%              
-%              end
-%              
-%              [B I] = sort((angle_start_end(2,:)));
-%              
-%              colorz = linspecer(length(angle_start_end(3,:)));
-% 
-%              subplot(5,4,G)
-%              hold on
-%              errorbar(angle_start_end(1,:),angle_start_end(2,:),sem_end(1,:),'Color','k','LineStyle','none') % angle at center of bin vs mean phase at end, errorbar
-%              scatter(angle_start_end(1,:),pred_angle,2.5,'k','filled') % angle at center of bin vs predicted angle
-%              
-%              for i = 1:length(angle_start_end(3,:))
-% 
-%                 if chan == 2
-% 
-%                     scatter(angle_start_end(1,i),angle_start_end(2,i),25,colorz(i,:),'filled','MarkerEdgeColor','k') % angle at center of bins vs mean at end 
-%                 else
-%                     scatter(angle_start_end(1,:),angle_start_end(2,:),25,colors(2,:))
-%                     errorbar(angle_start_end(1,:),angle_start_end(2,:),sem_end(1,:),'Color',colors(2,:),'LineStyle','none') 
-%                 end
-%              end
-%              
-%              G = G+1;
-% 
-%              ylim([-pi pi])
-%              xlim([-pi pi])
-%              title(['Delay: ',num2str(cycs/fs*1000),' ms'])     
-% 
-%              xlabel('Start Phase')
-%              ylabel('End Phase')
-%              
-%              subplot(5,4,G)
-%              hold on
-%              errorbar(angle_start_end(1,:),angle_start_end(3,:),sem_end(2,:),'Color','k','LineStyle','none') % angle at center of bin vs phase distance from endsamp to predicted angle, errorbar
-%              
-%              for i = 1:length(angle_start_end(3,:))
-%                 if chan == 2
-%                     scatter(angle_start_end(1,i),angle_start_end(3,i),25,colorz(i,:),'filled','MarkerEdgeColor','k')                         
-%                 else
-%                     scatter(angle_start_end(1,:),angle_start_end(3,:),25,colors(2,:))
-%                     errorbar(angle_start_end(1,:),angle_start_end(3,:),sem_end(2,:),'Color',colors(2,:),'LineStyle','none')    
-%                 end
-%              end
-%              
-%             G = G+1;
-%             a = 1;
-%             for A = [0 90 180 270]
-% %             for A = [340 60 160 240]
-%                 xline(wrapToPi(deg2rad(A)),'Color',colors(a,:)) 
-%                 a = a+1;
-%             end
-%              %xlim([-180 180])
-%              xlabel('Start Phase')
-%              ylabel('Phase Change')
-%              yline(0)
-%              text(.95*-pi,-.9*pi,'slower')
-%              %ylim([-200 200])
-%              %ylim([.1 .14])
-%              ylim([-pi pi])
-%              xlim([-pi pi])
-%              text(.95*-pi,.9*pi,'faster')
-%                          
-%                           
-%              subplot(5,4,G)
-%              polarplot([0 circ_mean(angle_start_end(2,:)')],[0 (circ_r(angle_start_end(2,:)',R_end'))],'k','LineWidth',2)
-%              hold on
-% 
-%             if chan == 2
-%                 B = 0.5;
-%                 for i = 1:length(angle_start_end(3,:))
-%                     polarscatter(angle_start_end(2,i),R_end(i),35,colorz(i,:),'filled','MarkerEdgeColor','k')
-%                     B = B+.025;
-%                     hold on
-%                     if G == 3
-%                         before(i) = angle_start_end(2,i);
-%                     else
-%   
-%                         hold on
-%                     end
-%                 end
-%                 
-%             else
-%                 polarplot(angle_start_end(1,:),(circ_dist(circ_mean(angle_start_end(3,:)'),angle_start_end(3,:)))+pi,'Color',colors(2,:))
-% 
-%             end
-%              hold on
-% 
-%              rlim([0 1.1])
-%              G = G+1;
-%              
-%              subplot(5,4,G)
-%              polarplot([0 circ_mean(angle_start_end(3,:)')],[0 (circ_r(angle_start_end(3,:)',R_end'))],'k','LineWidth',2)
-%              hold on
-% 
-%              if chan == 2
-%                 B = .5;
-%                 
-%                 for i = 1:length(angle_start_end(3,:))
-%                     polarscatter(angle_start_end(3,i),R_end(i),35,colorz(i,:),'filled','MarkerEdgeColor','k')
-%                     hold on
-%                     B = B+.025;
-%     
-%                     if G == 3
-%                         before(i) = angle_start_end(2,i);
-%                     else
-%                         hold on
-%                     end
-%                 end
-%              end
-%              
-%              hold on
-%              rlim([0 1.1])
-%              G = G+1;
-%              
-% end
-%              
-             
-% saveas(fig,[Savefolder,'Figure6_phase_response_curve_freq',num2str(osc_freq),'_phasic_REM.svg']);
-
-%% phase response curve - wake eve (eyes closed)
-
-% colors = linspecer(4);
-%              
-% %compute for many onset phases speeding up/slowing down
-%       
-% state = 2; % 1 = eyes open, 2 = eyes closed
-% chan = 2;
-%              
-% clear angle_start_end phase_after R_after angs clear sem_end pred_angle R_end
-% 
-% trigsamp = 501;
-% nbins = 10;
-% fs = 500;
-% osc_freq = 10;
-% 
-% % angs(1,:) = wrapToPi(deg2rad([0:360/10:360]));
-% angs(1,:) = wrapToPi(deg2rad([0:360/nbins:360]));
-% angs(2,:) = angs-pi/nbins;
-% angs(3,:) = angs(1,:)+pi/nbins;
-%              
-% fig = figure('Renderer','painters','units','normalized','outerposition',[0 0 .4 1])
-% hold on
-% % for chan = [12]
-% %                  clear phase_after phase_start_end
-% G = 1;
-% 
-% for cycs = [20/2 100/2 150/2 200/2 250/2]% [10 75 85 95 105];%75 100 125 150]%(25:25:200)
-%     A = 1;
-%     for ang = 1:size(angs,2) %0:18:360;%[36:10:360-36]
-%              
-%          for cond = 1
-%                  clear y
-%                  for s = 1:length(incl_sub)
-%                      
-%                      sub = incl_sub(s);
-%                      
-% %                         tmp = [];
-% %                         for condo = 1:4
-% %                          tmp = [tmp;group_phase.sub{sub}.chan{chan}.cond{condo}];
-% %                         tmp = ERP_trials_allsub_REM{sub}.trial_data_alphafilt_phase(ERP_trials_allsub_REM{sub}.phasic_ndx,:); % trials x samps
-% %                         tmp = ERP_trials_allsub_REM{sub}.trial_data_alphafilt_phase(ERP_trials_allsub_REM{sub}.goodtrial_ndx,:); % trials x samps
-% %                           tmp = ERP_trials_allsub_wake_m{sub}.trial_data_alphafilt_phase{state}(ERP_trials_allsub_wake_m{sub}.good_trial_ndx{state},:);
-%                           tmp = ERP_trials_allsub_wake_e{sub}.trial_data_alphafilt_phase{state}(ERP_trials_allsub_wake_e{sub}.good_trial_ndx{state},:);
-% %                         end
-%                                               
-%                         tmp = tmp(tmp(:,trigsamp)>angs(2,ang)&tmp(:,trigsamp)<angs(3,ang),:);
-%                          
-%                         if ~isempty(tmp)
-%                                       
-%                         Mn = circ_mean(tmp); % average across all trials in that bin
-%                         R = circ_r(tmp);
-%                         to_find = Mn(trigsamp);
-%                      
-%                         cosed = cos(circ_dist(Mn,to_find));
-%                                           
-%                         [pks,locs] = findpeaks(cosed);
-%                                        
-%                      
-%                         phase_before(s,cond) =  Mn(trigsamp);
-%                         cycle_length(s,cond) = abs(trigsamp-locs(find(locs==trigsamp)-1)); % duration from previous peak to stim peak in samps
-%                         phase_after(s,cond,1) = Mn(trigsamp);
-%                         R_after(s,cond,1) = R(trigsamp);
-%                                         
-%                     
-%                         D = (diff(locs(locs>=trigsamp-100 & locs<=trigsamp)));
-%                         %D = D(D>30);
-%                         cycles(s,cond,chan) = mean(D); % mean duration of all cycles between -200 ms and 0
-%                     
-%                     
-%                             for cyc = 1:round(cycs)
-%                                    
-%                                 phase_after(s,cond,cyc+1) = Mn(trigsamp+cyc); % mean phase for all samples from trigsamp to endsamp
-%                                 phase_change(s,cond,cyc) = circ_dist(phase_after(s,cond,cyc),phase_before(s,cond)); % phase distance from each samp after to trigsamp
-%                                 R_after(s,cond,cyc+1) = R(trigsamp+cyc);%  resultant for all samples from trigsamp to endsamp      
-%                         
-%                             end
-%                     
-%                                                                
-%                         else
-%                  
-%                             phase_after(s,cond,:) = NaN(1,1,cycs+1);
-%                             R_after(s,cond,:) = NaN(1,1,cycs+1);
-%                         end
-%             
-%                  end
-%            
-%          end
-%              
-%              
-% %              if cycs == trigsamp
-% %              phase_series(ang,:) = squeeze(circ_mean(phase_after));
-% %              end
-% %              clear cumul_phase_change
-%              
-%              cumul_phase_change = unwrap(squeeze(phase_after(:,1,:))')';
-% 
-%              angle_start_end(1,ang) = (wrapToPi(angs(1,ang))); % angle at center of bin
-%              tmp = phase_after(:,cond,end); % phase of endsamp
-%              tmp2 = R_after(:,cond,end); %  resultant of endsamp
-%              tmp2 = tmp2(~isnan(tmp2));
-%              tmp = tmp(~isnan(tmp));
-%                  
-%              angle_start_end(2,ang) = (wrapToPi(circ_mean(tmp))); % mean of endsamp phase across participants
-%              yy = num2str(cycs);
-%              yy = str2num(yy(end-1:end))/(fs/osc_freq); % calculate how much of a cycle this amount of time is
-%              pred_angle(ang) = wrapToPi(angs(1,ang)+ 2*yy*pi);
-%              R_end(1,ang) = circ_r(tmp,tmp2); % resultant of endsamp across participants
-%                 
-%              angle_start_end(3,ang) = (circ_dist(circ_mean(tmp),pred_angle(ang))); % phase distance from endsamp to predicted angle
-% 
-%              sem_end(1,ang) = circ_std(tmp)./sqrt(length(incl_sub)); % sem of end phase
-%              sem_end(2,ang) = circ_std(circ_dist(angs(1,ang),tmp))./sqrt(length(incl_sub)); % sem of difference between angle at center of bin and end phase
-% 
-%              A = A+1;
-%              
-%              end
-%              
-%              [B I] = sort((angle_start_end(2,:)));
-%              
-%              colorz = linspecer(length(angle_start_end(3,:)));
-% 
-%              subplot(5,4,G)
-%              hold on
-%              errorbar(angle_start_end(1,:),angle_start_end(2,:),sem_end(1,:),'Color','k','LineStyle','none') % angle at center of bin vs mean phase at end, errorbar
-%              scatter(angle_start_end(1,:),pred_angle,2.5,'k','filled') % angle at center of bin vs predicted angle
-%              
-%              for i = 1:length(angle_start_end(3,:))
-% 
-%                 if chan == 2
-% 
-%                     scatter(angle_start_end(1,i),angle_start_end(2,i),25,colorz(i,:),'filled','MarkerEdgeColor','k') % angle at center of bins vs mean at end 
-%                 else
-%                     scatter(angle_start_end(1,:),angle_start_end(2,:),25,colors(2,:))
-%                     errorbar(angle_start_end(1,:),angle_start_end(2,:),sem_end(1,:),'Color',colors(2,:),'LineStyle','none') 
-%                 end
-%              end
-%              
-%              G = G+1;
-% 
-%              ylim([-pi pi])
-%              xlim([-pi pi])
-%              title(['Delay: ',num2str(cycs/fs*1000),' ms'])     
-% 
-%              xlabel('Start Phase')
-%              ylabel('End Phase')
-%              
-%              subplot(5,4,G)
-%              hold on
-%              errorbar(angle_start_end(1,:),angle_start_end(3,:),sem_end(2,:),'Color','k','LineStyle','none') % angle at center of bin vs phase distance from endsamp to predicted angle, errorbar
-%              
-%              for i = 1:length(angle_start_end(3,:))
-%                 if chan == 2
-%                     scatter(angle_start_end(1,i),angle_start_end(3,i),25,colorz(i,:),'filled','MarkerEdgeColor','k')                         
-%                 else
-%                     scatter(angle_start_end(1,:),angle_start_end(3,:),25,colors(2,:))
-%                     errorbar(angle_start_end(1,:),angle_start_end(3,:),sem_end(2,:),'Color',colors(2,:),'LineStyle','none')    
-%                 end
-%              end
-%              
-%             G = G+1;
-%             a = 1;
-%             for A = [0 90 180 270]
-% %             for A = [340 60 160 240]
-%                 xline(wrapToPi(deg2rad(A)),'Color',colors(a,:)) 
-%                 a = a+1;
-%             end
-%              %xlim([-180 180])
-%              xlabel('Start Phase')
-%              ylabel('Phase Change')
-%              yline(0)
-%              text(.95*-pi,-.9*pi,'slower')
-%              %ylim([-200 200])
-%              %ylim([.1 .14])
-%              ylim([-pi pi])
-%              xlim([-pi pi])
-%              text(.95*-pi,.9*pi,'faster')
-%                          
-%                           
-%              subplot(5,4,G)
-%              polarplot([0 circ_mean(angle_start_end(2,:)')],[0 (circ_r(angle_start_end(2,:)',R_end'))],'k','LineWidth',2)
-%              hold on
-% 
-%             if chan == 2
-%                 B = 0.5;
-%                 for i = 1:length(angle_start_end(3,:))
-%                     polarscatter(angle_start_end(2,i),R_end(i),35,colorz(i,:),'filled','MarkerEdgeColor','k')
-%                     B = B+.025;
-%                     hold on
-%                     if G == 3
-%                         before(i) = angle_start_end(2,i);
-%                     else
-%   
-%                         hold on
-%                     end
-%                 end
-%                 
-%             else
-%                 polarplot(angle_start_end(1,:),(circ_dist(circ_mean(angle_start_end(3,:)'),angle_start_end(3,:)))+pi,'Color',colors(2,:))
-% 
-%             end
-%              hold on
-% 
-%              rlim([0 1.1])
-%              G = G+1;
-%              
-%              subplot(5,4,G)
-%              polarplot([0 circ_mean(angle_start_end(3,:)')],[0 (circ_r(angle_start_end(3,:)',R_end'))],'k','LineWidth',2)
-%              hold on
-% 
-%              if chan == 2
-%                 B = .5;
-%                 
-%                 for i = 1:length(angle_start_end(3,:))
-%                     polarscatter(angle_start_end(3,i),R_end(i),35,colorz(i,:),'filled','MarkerEdgeColor','k')
-%                     hold on
-%                     B = B+.025;
-%     
-%                     if G == 3
-%                         before(i) = angle_start_end(2,i);
-%                     else
-%                         hold on
-%                     end
-%                 end
-%              end
-%              
-%              hold on
-%              rlim([0 1.1])
-%              G = G+1;
-%              
-% end
-%              
-%              
-% % saveas(fig,[Savefolder,'Figure6_phase_response_curve_freq',num2str(osc_freq),'_wake_eve_EC.svg']);
-% 
-% %% phase response curve - wake mor (eyes closed)
-% 
-% colors = linspecer(4);
-%              
-% %compute for many onset phases speeding up/slowing down
-%       
-% state = 2; % 1 = eyes open, 2 = eyes closed
-% chan = 2;
-%              
-% clear angle_start_end phase_after R_after angs clear sem_end pred_angle R_end
-% 
-% trigsamp = 501;
-% nbins = 10;
-% fs = 500;
-% osc_freq = 10;
-% 
-% % angs(1,:) = wrapToPi(deg2rad([0:360/10:360]));
-% angs(1,:) = wrapToPi(deg2rad([0:360/nbins:360]));
-% angs(2,:) = angs-pi/nbins;
-% angs(3,:) = angs(1,:)+pi/nbins;
-%              
-% fig = figure('Renderer','painters','units','normalized','outerposition',[0 0 .4 1])
-% hold on
-% % for chan = [12]
-% %                  clear phase_after phase_start_end
-% G = 1;
-% 
-% for cycs = [20/2 100/2 150/2 200/2 250/2]% [10 75 85 95 105];%75 100 125 150]%(25:25:200)
-%     A = 1;
-%     for ang = 1:size(angs,2) %0:18:360;%[36:10:360-36]
-%              
-%          for cond = 1
-%                  clear y
-%                  for s = 1:length(incl_sub)
-%                      
-%                      sub = incl_sub(s);
-%                      
-% %                         tmp = [];
-% %                         for condo = 1:4
-% %                          tmp = [tmp;group_phase.sub{sub}.chan{chan}.cond{condo}];
-% %                         tmp = ERP_trials_allsub_REM{sub}.trial_data_alphafilt_phase(ERP_trials_allsub_REM{sub}.phasic_ndx,:); % trials x samps
-% %                         tmp = ERP_trials_allsub_REM{sub}.trial_data_alphafilt_phase(ERP_trials_allsub_REM{sub}.goodtrial_ndx,:); % trials x samps
-%                           tmp = ERP_trials_allsub_wake_m{sub}.trial_data_alphafilt_phase{state}(ERP_trials_allsub_wake_m{sub}.good_trial_ndx{state},:);
-% %                           tmp = ERP_trials_allsub_wake_e{sub}.trial_data_alphafilt_phase{state}(ERP_trials_allsub_wake_e{sub}.good_trial_ndx{state},:);
-% %                         end
-%                                               
-%                         tmp = tmp(tmp(:,trigsamp)>angs(2,ang)&tmp(:,trigsamp)<angs(3,ang),:);
-%                          
-%                         if ~isempty(tmp)
-%                                       
-%                         Mn = circ_mean(tmp); % average across all trials in that bin
-%                         R = circ_r(tmp);
-%                         to_find = Mn(trigsamp);
-%                      
-%                         cosed = cos(circ_dist(Mn,to_find));
-%                                           
-%                         [pks,locs] = findpeaks(cosed);
-%                                        
-%                      
-%                         phase_before(s,cond) =  Mn(trigsamp);
-%                         cycle_length(s,cond) = abs(trigsamp-locs(find(locs==trigsamp)-1)); % duration from previous peak to stim peak in samps
-%                         phase_after(s,cond,1) = Mn(trigsamp);
-%                         R_after(s,cond,1) = R(trigsamp);
-%                                         
-%                     
-%                         D = (diff(locs(locs>=trigsamp-100 & locs<=trigsamp)));
-%                         %D = D(D>30);
-%                         cycles(s,cond,chan) = mean(D); % mean duration of all cycles between -200 ms and 0
-%                     
-%                     
-%                             for cyc = 1:round(cycs)
-%                                    
-%                                 phase_after(s,cond,cyc+1) = Mn(trigsamp+cyc); % mean phase for all samples from trigsamp to endsamp
-%                                 phase_change(s,cond,cyc) = circ_dist(phase_after(s,cond,cyc),phase_before(s,cond)); % phase distance from each samp after to trigsamp
-%                                 R_after(s,cond,cyc+1) = R(trigsamp+cyc);%  resultant for all samples from trigsamp to endsamp      
-%                         
-%                             end
-%                     
-%                                                                
-%                         else
-%                  
-%                             phase_after(s,cond,:) = NaN(1,1,cycs+1);
-%                             R_after(s,cond,:) = NaN(1,1,cycs+1);
-%                         end
-%             
-%                  end
-%            
-%          end
-%              
-%              
-% %              if cycs == trigsamp
-% %              phase_series(ang,:) = squeeze(circ_mean(phase_after));
-% %              end
-% %              clear cumul_phase_change
-%              
-%              cumul_phase_change = unwrap(squeeze(phase_after(:,1,:))')';
-% 
-%              angle_start_end(1,ang) = (wrapToPi(angs(1,ang))); % angle at center of bin
-%              tmp = phase_after(:,cond,end); % phase of endsamp
-%              tmp2 = R_after(:,cond,end); %  resultant of endsamp
-%              tmp2 = tmp2(~isnan(tmp2));
-%              tmp = tmp(~isnan(tmp));
-%                  
-%              angle_start_end(2,ang) = (wrapToPi(circ_mean(tmp))); % mean of endsamp phase across participants
-%              yy = num2str(cycs);
-%              yy = str2num(yy(end-1:end))/(fs/osc_freq); % calculate how much of a cycle this amount of time is
-%              pred_angle(ang) = wrapToPi(angs(1,ang)+ 2*yy*pi);
-%              R_end(1,ang) = circ_r(tmp,tmp2); % resultant of endsamp across participants
-%                 
-%              angle_start_end(3,ang) = (circ_dist(circ_mean(tmp),pred_angle(ang))); % phase distance from endsamp to predicted angle
-% 
-%              sem_end(1,ang) = circ_std(tmp)./sqrt(length(incl_sub)); % sem of end phase
-%              sem_end(2,ang) = circ_std(circ_dist(angs(1,ang),tmp))./sqrt(length(incl_sub)); % sem of difference between angle at center of bin and end phase
-% 
-%              A = A+1;
-%              
-%              end
-%              
-%              [B I] = sort((angle_start_end(2,:)));
-%              
-%              colorz = linspecer(length(angle_start_end(3,:)));
-% 
-%              subplot(5,4,G)
-%              hold on
-%              errorbar(angle_start_end(1,:),angle_start_end(2,:),sem_end(1,:),'Color','k','LineStyle','none') % angle at center of bin vs mean phase at end, errorbar
-%              scatter(angle_start_end(1,:),pred_angle,2.5,'k','filled') % angle at center of bin vs predicted angle
-%              
-%              for i = 1:length(angle_start_end(3,:))
-% 
-%                 if chan == 2
-% 
-%                     scatter(angle_start_end(1,i),angle_start_end(2,i),25,colorz(i,:),'filled','MarkerEdgeColor','k') % angle at center of bins vs mean at end 
-%                 else
-%                     scatter(angle_start_end(1,:),angle_start_end(2,:),25,colors(2,:))
-%                     errorbar(angle_start_end(1,:),angle_start_end(2,:),sem_end(1,:),'Color',colors(2,:),'LineStyle','none') 
-%                 end
-%              end
-%              
-%              G = G+1;
-% 
-%              ylim([-pi pi])
-%              xlim([-pi pi])
-%              title(['Delay: ',num2str(cycs/fs*1000),' ms'])     
-% 
-%              xlabel('Start Phase')
-%              ylabel('End Phase')
-%              
-%              subplot(5,4,G)
-%              hold on
-%              errorbar(angle_start_end(1,:),angle_start_end(3,:),sem_end(2,:),'Color','k','LineStyle','none') % angle at center of bin vs phase distance from endsamp to predicted angle, errorbar
-%              
-%              for i = 1:length(angle_start_end(3,:))
-%                 if chan == 2
-%                     scatter(angle_start_end(1,i),angle_start_end(3,i),25,colorz(i,:),'filled','MarkerEdgeColor','k')                         
-%                 else
-%                     scatter(angle_start_end(1,:),angle_start_end(3,:),25,colors(2,:))
-%                     errorbar(angle_start_end(1,:),angle_start_end(3,:),sem_end(2,:),'Color',colors(2,:),'LineStyle','none')    
-%                 end
-%              end
-%              
-%             G = G+1;
-%             a = 1;
-%             for A = [0 90 180 270]
-% %             for A = [340 60 160 240]
-%                 xline(wrapToPi(deg2rad(A)),'Color',colors(a,:)) 
-%                 a = a+1;
-%             end
-%              %xlim([-180 180])
-%              xlabel('Start Phase')
-%              ylabel('Phase Change')
-%              yline(0)
-%              text(.95*-pi,-.9*pi,'slower')
-%              %ylim([-200 200])
-%              %ylim([.1 .14])
-%              ylim([-pi pi])
-%              xlim([-pi pi])
-%              text(.95*-pi,.9*pi,'faster')
-%                          
-%                           
-%              subplot(5,4,G)
-%              polarplot([0 circ_mean(angle_start_end(2,:)')],[0 (circ_r(angle_start_end(2,:)',R_end'))],'k','LineWidth',2)
-%              hold on
-% 
-%             if chan == 2
-%                 B = 0.5;
-%                 for i = 1:length(angle_start_end(3,:))
-%                     polarscatter(angle_start_end(2,i),R_end(i),35,colorz(i,:),'filled','MarkerEdgeColor','k')
-%                     B = B+.025;
-%                     hold on
-%                     if G == 3
-%                         before(i) = angle_start_end(2,i);
-%                     else
-%   
-%                         hold on
-%                     end
-%                 end
-%                 
-%             else
-%                 polarplot(angle_start_end(1,:),(circ_dist(circ_mean(angle_start_end(3,:)'),angle_start_end(3,:)))+pi,'Color',colors(2,:))
-% 
-%             end
-%              hold on
-% 
-%              rlim([0 1.1])
-%              G = G+1;
-%              
-%              subplot(5,4,G)
-%              polarplot([0 circ_mean(angle_start_end(3,:)')],[0 (circ_r(angle_start_end(3,:)',R_end'))],'k','LineWidth',2)
-%              hold on
-% 
-%              if chan == 2
-%                 B = .5;
-%                 
-%                 for i = 1:length(angle_start_end(3,:))
-%                     polarscatter(angle_start_end(3,i),R_end(i),35,colorz(i,:),'filled','MarkerEdgeColor','k')
-%                     hold on
-%                     B = B+.025;
-%     
-%                     if G == 3
-%                         before(i) = angle_start_end(2,i);
-%                     else
-%                         hold on
-%                     end
-%                 end
-%              end
-%              
-%              hold on
-%              rlim([0 1.1])
-%              G = G+1;
-%              
-% end
-%              
-%              
-% % saveas(fig,[Savefolder,'Figure6_phase_response_curve_freq',num2str(osc_freq),'_wake_mor_EC.svg']);
-% 
-%              
-%              
-% 
-% 
-%              
-%              
-% 
diff --git a/Figures/not_used/ERPs_phase_response_theta.m b/Figures/not_used/ERPs_phase_response_theta.m
deleted file mode 100644
index 4a5c752ccbcfe75ab672fae51f986171a5dd4b5a..0000000000000000000000000000000000000000
--- a/Figures/not_used/ERPs_phase_response_theta.m
+++ /dev/null
@@ -1,1002 +0,0 @@
-clear all;
-close all;
-
-addpath(genpath('/user/HS301/m17462/matlab/eeglab'));
-addpath(genpath('/user/HS301/m17462/matlab/Scripts/RSN'));
-addpath(genpath('/user/HS301/m17462/matlab/Henry/useful_functions'));
-addpath(genpath('/user/HS301/m17462/matlab/colorGradient'));
-
-Savefolder = '/vol/research/nemo/datasets/RSN/data/analysis/Figures/';
-
-incl_sub = setdiff(1:19,12);
-
-Folderpath = '/vol/research/nemo/datasets/RSN/data/hdEEG/';
-sub_Folderpath = dir([Folderpath,'RSN*']);
-
-%% load all trial data - REM
-
-for s = 1:length(sub_Folderpath)
-    
-    erp_good_file = dir([Folderpath,sub_Folderpath(s).name,filesep,'*_sleep*_erp_good.mat']);
-    
-    load([Folderpath,sub_Folderpath(s).name,filesep,erp_good_file(1).name]);
-    
-    ERP_trials_allsub_REM{s} = ERP; 
-    
-end
-
-
-%% load all trial data - wake eve
-
-for s = 1:length(sub_Folderpath)
-    
-    erp_good_file = dir([Folderpath,sub_Folderpath(s).name,filesep,'*_wake_e*_erp_good.mat']);
-    
-    load([Folderpath,sub_Folderpath(s).name,filesep,erp_good_file(1).name]);
-    
-    ERP_trials_allsub_wake_e{s} = ERP; 
-    
-end
-
-%% load all trial data - wake mor
-
-for s = 1:length(sub_Folderpath)
-    
-    erp_good_file = dir([Folderpath,sub_Folderpath(s).name,filesep,'*_wake_m*_erp_good.mat']);
-    
-    load([Folderpath,sub_Folderpath(s).name,filesep,erp_good_file(1).name]);
-    
-    ERP_trials_allsub_wake_m{s} = ERP; 
-    
-end
-
-%% phase response curve - tonic REM
-
-colors = linspecer(4);
-             
-%compute for many onset phases speeding up/slowing down
-      
-state = 1; % 1 = eyes open, 2 = eyes closed
-chan = 2;
-             
-clear angle_start_end phase_after R_after angs clear sem_end pred_angle R_end
-
-trigsamp = 501;
-nbins = 10;
-fs = 500;
-osc_freq = 6;
-
-% angs(1,:) = wrapToPi(deg2rad([0:360/10:360]));
-angs(1,:) = wrapToPi(deg2rad([0:360/nbins:360]));
-angs(2,:) = angs-pi/nbins;
-angs(3,:) = angs(1,:)+pi/nbins;
-             
-fig = figure('Renderer','painters','units','normalized','outerposition',[0 0 .4 1])
-hold on
-% for chan = [12]
-%                  clear phase_after phase_start_end
-G = 1;
-
-cycs_all = [20/2 100/2 150/2 200/2 250/2];
-
-for angle = 1:length(cycs_all) % [10 75 85 95 105];%75 100 125 150]%(25:25:200)
-    cycs = cycs_all(angle);
-    A = 1;
-    for ang = 1:size(angs,2) %0:18:360;%[36:10:360-36]
-             
-         for cond = 1
-                 clear y
-                 for s = 1:length(incl_sub)
-                     
-                     sub = incl_sub(s);
-                     
-%                         tmp = [];
-%                         for condo = 1:4
-%                          tmp = [tmp;group_phase.sub{sub}.chan{chan}.cond{condo}];
-                        tmp = ERP_trials_allsub_REM{sub}.trial_data_thetafilt_phase(ERP_trials_allsub_REM{sub}.tonic_ndx,:); % trials x samps
-%                         tmp = ERP_trials_allsub_REM{sub}.trial_data_alphafilt_phase(ERP_trials_allsub_REM{sub}.goodtrial_ndx,:); % trials x samps
-%                           tmp = ERP_trials_allsub_wake_m{sub}.trial_data_alphafilt_phase{state}(ERP_trials_allsub_wake_m{sub}.good_trial_ndx{state},:);
-%                           tmp = ERP_trials_allsub_wake_e{sub}.trial_data_alphafilt_phase{state}(ERP_trials_allsub_wake_e{sub}.good_trial_ndx{state},:);
-%                         end
-                                              
-                        tmp = tmp(tmp(:,trigsamp)>angs(2,ang)&tmp(:,trigsamp)<angs(3,ang),:);
-                         
-                        if ~isempty(tmp)
-                                      
-                        Mn = circ_mean(tmp); % average across all trials in that bin
-                        R = circ_r(tmp);
-                        to_find = Mn(trigsamp);
-                     
-                        cosed = cos(circ_dist(Mn,to_find));
-                                          
-                        [pks,locs] = findpeaks(cosed);
-                                       
-                     
-                        phase_before(s,cond) =  Mn(trigsamp);
-                        cycle_length(s,cond) = abs(trigsamp-locs(find(locs==trigsamp)-1)); % duration from previous peak to stim peak in samps
-                        phase_after(s,cond,1) = Mn(trigsamp);
-                        R_after(s,cond,1) = R(trigsamp);
-                                        
-                    
-                        D = (diff(locs(locs>=trigsamp-100 & locs<=trigsamp)));
-                        %D = D(D>30);
-                        cycles(s,cond,chan) = mean(D); % mean duration of all cycles between -200 ms and 0
-                    
-                    
-                            for cyc = 1:round(cycs)
-                                   
-                                phase_after(s,cond,cyc+1) = Mn(trigsamp+cyc); % mean phase for all samples from trigsamp to endsamp
-                                phase_change(s,cond,cyc) = circ_dist(phase_after(s,cond,cyc),phase_before(s,cond)); % phase distance from each samp after to trigsamp
-                                R_after(s,cond,cyc+1) = R(trigsamp+cyc);%  resultant for all samples from trigsamp to endsamp      
-                        
-                            end
-                    
-                                                               
-                        else
-                 
-                            phase_after(s,cond,:) = NaN(1,1,cycs+1);
-                            R_after(s,cond,:) = NaN(1,1,cycs+1);
-                        end
-            
-                 end
-           
-         end
-             
-             
-%              if cycs == trigsamp
-%              phase_series(ang,:) = squeeze(circ_mean(phase_after));
-%              end
-%              clear cumul_phase_change
-             
-             cumul_phase_change = unwrap(squeeze(phase_after(:,1,:))')';
-
-             angle_start_end(1,ang) = (wrapToPi(angs(1,ang))); % angle at center of bin
-             tmp = phase_after(:,cond,end); % phase of endsamp
-             tmp2 = R_after(:,cond,end); %  resultant of endsamp
-             tmp2 = tmp2(~isnan(tmp2));
-             tmp = tmp(~isnan(tmp));
-                 
-             angle_start_end(2,ang) = (wrapToPi(circ_mean(tmp))); % mean of endsamp phase across participants
-             yy = num2str(cycs);
-             yy = str2num(yy(end-1:end))/(fs/osc_freq); % calculate how much of a cycle this amount of time is
-             pred_angle(ang) = wrapToPi(angs(1,ang)+ 2*yy*pi);
-             R_end(1,ang) = circ_r(tmp,tmp2); % resultant of endsamp across participants
-                
-             angle_start_end(3,ang) = (circ_dist(circ_mean(tmp),pred_angle(ang))); % phase distance from endsamp to predicted angle
-
-             sem_end(1,ang) = circ_std(tmp)./sqrt(length(incl_sub)); % sem of end phase
-             sem_end(2,ang) = circ_std(circ_dist(angs(1,ang),tmp))./sqrt(length(incl_sub)); % sem of difference between angle at center of bin and end phase
-
-             A = A+1;
-             
-             end
-             
-             [B I] = sort((angle_start_end(2,:)));
-             
-             colorz = linspecer(length(angle_start_end(3,:)));
-
-             subplot(5,4,(angle-1)*4+1)
-             hold on
-             errorbar(angle_start_end(1,:),angle_start_end(2,:),sem_end(1,:),'Color','k','LineStyle','none') % angle at center of bin vs mean phase at end, errorbar
-             scatter(angle_start_end(1,:),pred_angle,2.5,'k','filled') % angle at center of bin vs predicted angle
-             
-             for i = 1:length(angle_start_end(3,:))
-
-                if chan == 2
-
-                    scatter(angle_start_end(1,i),angle_start_end(2,i),25,colorz(i,:),'filled','MarkerEdgeColor','k') % angle at center of bins vs mean at end 
-                else
-                    scatter(angle_start_end(1,:),angle_start_end(2,:),25,colors(2,:))
-                    errorbar(angle_start_end(1,:),angle_start_end(2,:),sem_end(1,:),'Color',colors(2,:),'LineStyle','none') 
-                end
-             end
-             
-%              G = G+1;
-
-             ylim([-pi pi])
-             xlim([-pi pi])
-             title(['Delay: ',num2str(cycs/fs*1000),' ms'])     
-
-             xlabel('Start Phase')
-             ylabel('End Phase')
-             
-             subplot(5,4,(angle-1)*4+2)
-             hold on
-             errorbar(angle_start_end(1,:),angle_start_end(3,:),sem_end(2,:),'Color','k','LineStyle','none') % angle at center of bin vs phase distance from endsamp to predicted angle, errorbar
-             
-             for i = 1:length(angle_start_end(3,:))
-                if chan == 2
-                    scatter(angle_start_end(1,i),angle_start_end(3,i),25,colorz(i,:),'filled','MarkerEdgeColor','k')                         
-                else
-                    scatter(angle_start_end(1,:),angle_start_end(3,:),25,colors(2,:))
-                    errorbar(angle_start_end(1,:),angle_start_end(3,:),sem_end(2,:),'Color',colors(2,:),'LineStyle','none')    
-                end
-             end
-             
-%             G = G+1;
-            a = 1;
-            for A = [0 90 180 270]
-%             for A = [340 60 160 240]
-                xline(wrapToPi(deg2rad(A)),'Color',colors(a,:)) 
-                a = a+1;
-            end
-             %xlim([-180 180])
-             xlabel('Start Phase')
-             ylabel('Phase Change')
-             yline(0)
-%              text(.95*-pi,-.9*pi,'slower')
-             %ylim([-200 200])
-             %ylim([.1 .14])
-             ylim([-pi pi])
-             xlim([-pi pi])
-%              text(.95*-pi,.9*pi,'faster')
-                         
-                          
-%              subplot(5,4,G)
-%              polarplot([0 circ_mean(angle_start_end(2,:)')],[0 (circ_r(angle_start_end(2,:)',R_end'))],'k','LineWidth',2)
-%              hold on
-% 
-%             if chan == 2
-%                 B = 0.5;
-%                 for i = 1:length(angle_start_end(3,:))
-%                     polarscatter(angle_start_end(2,i),R_end(i),35,colorz(i,:),'filled','MarkerEdgeColor','k')
-%                     B = B+.025;
-%                     hold on
-%                     if G == 3
-%                         before(i) = angle_start_end(2,i);
-%                     else
-%   
-%                         hold on
-%                     end
-%                 end
-%                 
-%             else
-%                 polarplot(angle_start_end(1,:),(circ_dist(circ_mean(angle_start_end(3,:)'),angle_start_end(3,:)))+pi,'Color',colors(2,:))
-% 
-%             end
-%              hold on
-% 
-%              rlim([0 1.1])
-%              G = G+1;
-%              
-%              subplot(5,4,G)
-%              polarplot([0 circ_mean(angle_start_end(3,:)')],[0 (circ_r(angle_start_end(3,:)',R_end'))],'k','LineWidth',2)
-%              hold on
-% 
-%              if chan == 2
-%                 B = .5;
-%                 
-%                 for i = 1:length(angle_start_end(3,:))
-%                     polarscatter(angle_start_end(3,i),R_end(i),35,colorz(i,:),'filled','MarkerEdgeColor','k')
-%                     hold on
-%                     B = B+.025;
-%     
-%                     if G == 3
-%                         before(i) = angle_start_end(2,i);
-%                     else
-%                         hold on
-%                     end
-%                 end
-%              end
-%              
-%              hold on
-%              rlim([0 1.1])
-%              G = G+1;
-             
-end
-             
-             
-saveas(fig,[Savefolder,'Figure6_phase_response_curve_freq',num2str(osc_freq),'_thetafilt_tonic_REM.svg']);
-
-%% phase response curve - phasic REM
-
-% colors = linspecer(4);
-%              
-% %compute for many onset phases speeding up/slowing down
-%       
-% state = 1; % 1 = eyes open, 2 = eyes closed
-% chan = 2;
-%              
-% clear angle_start_end phase_after R_after angs clear sem_end pred_angle R_end
-% 
-% trigsamp = 501;
-% nbins = 10;
-% fs = 500;
-% osc_freq = 6;
-% 
-% % angs(1,:) = wrapToPi(deg2rad([0:360/10:360]));
-% angs(1,:) = wrapToPi(deg2rad([0:360/nbins:360]));
-% angs(2,:) = angs-pi/nbins;
-% angs(3,:) = angs(1,:)+pi/nbins;
-%              
-% fig = figure('Renderer','painters','units','normalized','outerposition',[0 0 .4 1])
-% hold on
-% % for chan = [12]
-% %                  clear phase_after phase_start_end
-% G = 1;
-% 
-% for cycs = [20/2 100/2 150/2 200/2 250/2]% [10 75 85 95 105];%75 100 125 150]%(25:25:200)
-%     A = 1;
-%     for ang = 1:size(angs,2) %0:18:360;%[36:10:360-36]
-%              
-%          for cond = 1
-%                  clear y
-%                  for s = 1:length(incl_sub)
-%                      
-%                      sub = incl_sub(s);
-%                      
-% %                         tmp = [];
-% %                         for condo = 1:4
-% %                          tmp = [tmp;group_phase.sub{sub}.chan{chan}.cond{condo}];
-%                         tmp = ERP_trials_allsub_REM{sub}.trial_data_thetafilt_phase(ERP_trials_allsub_REM{sub}.phasic_ndx,:); % trials x samps
-% %                         tmp = ERP_trials_allsub_REM{sub}.trial_data_alphafilt_phase(ERP_trials_allsub_REM{sub}.goodtrial_ndx,:); % trials x samps
-% %                           tmp = ERP_trials_allsub_wake_m{sub}.trial_data_alphafilt_phase{state}(ERP_trials_allsub_wake_m{sub}.good_trial_ndx{state},:);
-% %                           tmp = ERP_trials_allsub_wake_e{sub}.trial_data_alphafilt_phase{state}(ERP_trials_allsub_wake_e{sub}.good_trial_ndx{state},:);
-% %                         end
-%                                               
-%                         tmp = tmp(tmp(:,trigsamp)>angs(2,ang)&tmp(:,trigsamp)<angs(3,ang),:);
-%                          
-%                         if ~isempty(tmp)
-%                                       
-%                         Mn = circ_mean(tmp); % average across all trials in that bin
-%                         R = circ_r(tmp);
-%                         to_find = Mn(trigsamp);
-%                      
-%                         cosed = cos(circ_dist(Mn,to_find));
-%                                           
-%                         [pks,locs] = findpeaks(cosed);
-%                                        
-%                      
-%                         phase_before(s,cond) =  Mn(trigsamp);
-%                         cycle_length(s,cond) = abs(trigsamp-locs(find(locs==trigsamp)-1)); % duration from previous peak to stim peak in samps
-%                         phase_after(s,cond,1) = Mn(trigsamp);
-%                         R_after(s,cond,1) = R(trigsamp);
-%                                         
-%                     
-%                         D = (diff(locs(locs>=trigsamp-100 & locs<=trigsamp)));
-%                         %D = D(D>30);
-%                         cycles(s,cond,chan) = mean(D); % mean duration of all cycles between -200 ms and 0
-%                     
-%                     
-%                             for cyc = 1:round(cycs)
-%                                    
-%                                 phase_after(s,cond,cyc+1) = Mn(trigsamp+cyc); % mean phase for all samples from trigsamp to endsamp
-%                                 phase_change(s,cond,cyc) = circ_dist(phase_after(s,cond,cyc),phase_before(s,cond)); % phase distance from each samp after to trigsamp
-%                                 R_after(s,cond,cyc+1) = R(trigsamp+cyc);%  resultant for all samples from trigsamp to endsamp      
-%                         
-%                             end
-%                     
-%                                                                
-%                         else
-%                  
-%                             phase_after(s,cond,:) = NaN(1,1,cycs+1);
-%                             R_after(s,cond,:) = NaN(1,1,cycs+1);
-%                         end
-%             
-%                  end
-%            
-%          end
-%              
-%              
-% %              if cycs == trigsamp
-% %              phase_series(ang,:) = squeeze(circ_mean(phase_after));
-% %              end
-% %              clear cumul_phase_change
-%              
-%              cumul_phase_change = unwrap(squeeze(phase_after(:,1,:))')';
-% 
-%              angle_start_end(1,ang) = (wrapToPi(angs(1,ang))); % angle at center of bin
-%              tmp = phase_after(:,cond,end); % phase of endsamp
-%              tmp2 = R_after(:,cond,end); %  resultant of endsamp
-%              tmp2 = tmp2(~isnan(tmp2));
-%              tmp = tmp(~isnan(tmp));
-%                  
-%              angle_start_end(2,ang) = (wrapToPi(circ_mean(tmp))); % mean of endsamp phase across participants
-%              yy = num2str(cycs);
-%              yy = str2num(yy(end-1:end))/(fs/osc_freq); % calculate how much of a cycle this amount of time is
-%              pred_angle(ang) = wrapToPi(angs(1,ang)+ 2*yy*pi);
-%              R_end(1,ang) = circ_r(tmp,tmp2); % resultant of endsamp across participants
-%                 
-%              angle_start_end(3,ang) = (circ_dist(circ_mean(tmp),pred_angle(ang))); % phase distance from endsamp to predicted angle
-% 
-%              sem_end(1,ang) = circ_std(tmp)./sqrt(length(incl_sub)); % sem of end phase
-%              sem_end(2,ang) = circ_std(circ_dist(angs(1,ang),tmp))./sqrt(length(incl_sub)); % sem of difference between angle at center of bin and end phase
-% 
-%              A = A+1;
-%              
-%              end
-%              
-%              [B I] = sort((angle_start_end(2,:)));
-%              
-%              colorz = linspecer(length(angle_start_end(3,:)));
-% 
-%              subplot(5,4,G)
-%              hold on
-%              errorbar(angle_start_end(1,:),angle_start_end(2,:),sem_end(1,:),'Color','k','LineStyle','none') % angle at center of bin vs mean phase at end, errorbar
-%              scatter(angle_start_end(1,:),pred_angle,2.5,'k','filled') % angle at center of bin vs predicted angle
-%              
-%              for i = 1:length(angle_start_end(3,:))
-% 
-%                 if chan == 2
-% 
-%                     scatter(angle_start_end(1,i),angle_start_end(2,i),25,colorz(i,:),'filled','MarkerEdgeColor','k') % angle at center of bins vs mean at end 
-%                 else
-%                     scatter(angle_start_end(1,:),angle_start_end(2,:),25,colors(2,:))
-%                     errorbar(angle_start_end(1,:),angle_start_end(2,:),sem_end(1,:),'Color',colors(2,:),'LineStyle','none') 
-%                 end
-%              end
-%              
-%              G = G+1;
-% 
-%              ylim([-pi pi])
-%              xlim([-pi pi])
-%              title(['Delay: ',num2str(cycs/fs*1000),' ms'])     
-% 
-%              xlabel('Start Phase')
-%              ylabel('End Phase')
-%              
-%              subplot(5,4,G)
-%              hold on
-%              errorbar(angle_start_end(1,:),angle_start_end(3,:),sem_end(2,:),'Color','k','LineStyle','none') % angle at center of bin vs phase distance from endsamp to predicted angle, errorbar
-%              
-%              for i = 1:length(angle_start_end(3,:))
-%                 if chan == 2
-%                     scatter(angle_start_end(1,i),angle_start_end(3,i),25,colorz(i,:),'filled','MarkerEdgeColor','k')                         
-%                 else
-%                     scatter(angle_start_end(1,:),angle_start_end(3,:),25,colors(2,:))
-%                     errorbar(angle_start_end(1,:),angle_start_end(3,:),sem_end(2,:),'Color',colors(2,:),'LineStyle','none')    
-%                 end
-%              end
-%              
-%             G = G+1;
-%             a = 1;
-%             for A = [0 90 180 270]
-% %             for A = [340 60 160 240]
-%                 xline(wrapToPi(deg2rad(A)),'Color',colors(a,:)) 
-%                 a = a+1;
-%             end
-%              %xlim([-180 180])
-%              xlabel('Start Phase')
-%              ylabel('Phase Change')
-%              yline(0)
-%              text(.95*-pi,-.9*pi,'slower')
-%              %ylim([-200 200])
-%              %ylim([.1 .14])
-%              ylim([-pi pi])
-%              xlim([-pi pi])
-%              text(.95*-pi,.9*pi,'faster')
-%                          
-%                           
-%              subplot(5,4,G)
-%              polarplot([0 circ_mean(angle_start_end(2,:)')],[0 (circ_r(angle_start_end(2,:)',R_end'))],'k','LineWidth',2)
-%              hold on
-% 
-%             if chan == 2
-%                 B = 0.5;
-%                 for i = 1:length(angle_start_end(3,:))
-%                     polarscatter(angle_start_end(2,i),R_end(i),35,colorz(i,:),'filled','MarkerEdgeColor','k')
-%                     B = B+.025;
-%                     hold on
-%                     if G == 3
-%                         before(i) = angle_start_end(2,i);
-%                     else
-%   
-%                         hold on
-%                     end
-%                 end
-%                 
-%             else
-%                 polarplot(angle_start_end(1,:),(circ_dist(circ_mean(angle_start_end(3,:)'),angle_start_end(3,:)))+pi,'Color',colors(2,:))
-% 
-%             end
-%              hold on
-% 
-%              rlim([0 1.1])
-%              G = G+1;
-%              
-%              subplot(5,4,G)
-%              polarplot([0 circ_mean(angle_start_end(3,:)')],[0 (circ_r(angle_start_end(3,:)',R_end'))],'k','LineWidth',2)
-%              hold on
-% 
-%              if chan == 2
-%                 B = .5;
-%                 
-%                 for i = 1:length(angle_start_end(3,:))
-%                     polarscatter(angle_start_end(3,i),R_end(i),35,colorz(i,:),'filled','MarkerEdgeColor','k')
-%                     hold on
-%                     B = B+.025;
-%     
-%                     if G == 3
-%                         before(i) = angle_start_end(2,i);
-%                     else
-%                         hold on
-%                     end
-%                 end
-%              end
-%              
-%              hold on
-%              rlim([0 1.1])
-%              G = G+1;
-%              
-% end
-%              
-%              
-% saveas(fig,[Savefolder,'Figure6_phase_response_curve_freq',num2str(osc_freq),'_thetafilt_phasic_REM.svg']);
-
-%% phase response curve - wake eve (eyes closed)
-
-% colors = linspecer(4);
-%              
-% %compute for many onset phases speeding up/slowing down
-%       
-% state = 2; % 1 = eyes open, 2 = eyes closed
-% chan = 2;
-%              
-% clear angle_start_end phase_after R_after angs clear sem_end pred_angle R_end
-% 
-% trigsamp = 501;
-% nbins = 10;
-% fs = 500;
-% osc_freq = 10;
-% 
-% % angs(1,:) = wrapToPi(deg2rad([0:360/10:360]));
-% angs(1,:) = wrapToPi(deg2rad([0:360/nbins:360]));
-% angs(2,:) = angs-pi/nbins;
-% angs(3,:) = angs(1,:)+pi/nbins;
-%              
-% fig = figure('Renderer','painters','units','normalized','outerposition',[0 0 .4 1])
-% hold on
-% % for chan = [12]
-% %                  clear phase_after phase_start_end
-% G = 1;
-% 
-% for cycs = [20/2 100/2 150/2 200/2 250/2]% [10 75 85 95 105];%75 100 125 150]%(25:25:200)
-%     A = 1;
-%     for ang = 1:size(angs,2) %0:18:360;%[36:10:360-36]
-%              
-%          for cond = 1
-%                  clear y
-%                  for s = 1:length(incl_sub)
-%                      
-%                      sub = incl_sub(s);
-%                      
-% %                         tmp = [];
-% %                         for condo = 1:4
-% %                          tmp = [tmp;group_phase.sub{sub}.chan{chan}.cond{condo}];
-% %                         tmp = ERP_trials_allsub_REM{sub}.trial_data_alphafilt_phase(ERP_trials_allsub_REM{sub}.phasic_ndx,:); % trials x samps
-% %                         tmp = ERP_trials_allsub_REM{sub}.trial_data_alphafilt_phase(ERP_trials_allsub_REM{sub}.goodtrial_ndx,:); % trials x samps
-% %                           tmp = ERP_trials_allsub_wake_m{sub}.trial_data_alphafilt_phase{state}(ERP_trials_allsub_wake_m{sub}.good_trial_ndx{state},:);
-%                           tmp = ERP_trials_allsub_wake_e{sub}.trial_data_alphafilt_phase{state}(ERP_trials_allsub_wake_e{sub}.good_trial_ndx{state},:);
-% %                         end
-%                                               
-%                         tmp = tmp(tmp(:,trigsamp)>angs(2,ang)&tmp(:,trigsamp)<angs(3,ang),:);
-%                          
-%                         if ~isempty(tmp)
-%                                       
-%                         Mn = circ_mean(tmp); % average across all trials in that bin
-%                         R = circ_r(tmp);
-%                         to_find = Mn(trigsamp);
-%                      
-%                         cosed = cos(circ_dist(Mn,to_find));
-%                                           
-%                         [pks,locs] = findpeaks(cosed);
-%                                        
-%                      
-%                         phase_before(s,cond) =  Mn(trigsamp);
-%                         cycle_length(s,cond) = abs(trigsamp-locs(find(locs==trigsamp)-1)); % duration from previous peak to stim peak in samps
-%                         phase_after(s,cond,1) = Mn(trigsamp);
-%                         R_after(s,cond,1) = R(trigsamp);
-%                                         
-%                     
-%                         D = (diff(locs(locs>=trigsamp-100 & locs<=trigsamp)));
-%                         %D = D(D>30);
-%                         cycles(s,cond,chan) = mean(D); % mean duration of all cycles between -200 ms and 0
-%                     
-%                     
-%                             for cyc = 1:round(cycs)
-%                                    
-%                                 phase_after(s,cond,cyc+1) = Mn(trigsamp+cyc); % mean phase for all samples from trigsamp to endsamp
-%                                 phase_change(s,cond,cyc) = circ_dist(phase_after(s,cond,cyc),phase_before(s,cond)); % phase distance from each samp after to trigsamp
-%                                 R_after(s,cond,cyc+1) = R(trigsamp+cyc);%  resultant for all samples from trigsamp to endsamp      
-%                         
-%                             end
-%                     
-%                                                                
-%                         else
-%                  
-%                             phase_after(s,cond,:) = NaN(1,1,cycs+1);
-%                             R_after(s,cond,:) = NaN(1,1,cycs+1);
-%                         end
-%             
-%                  end
-%            
-%          end
-%              
-%              
-% %              if cycs == trigsamp
-% %              phase_series(ang,:) = squeeze(circ_mean(phase_after));
-% %              end
-% %              clear cumul_phase_change
-%              
-%              cumul_phase_change = unwrap(squeeze(phase_after(:,1,:))')';
-% 
-%              angle_start_end(1,ang) = (wrapToPi(angs(1,ang))); % angle at center of bin
-%              tmp = phase_after(:,cond,end); % phase of endsamp
-%              tmp2 = R_after(:,cond,end); %  resultant of endsamp
-%              tmp2 = tmp2(~isnan(tmp2));
-%              tmp = tmp(~isnan(tmp));
-%                  
-%              angle_start_end(2,ang) = (wrapToPi(circ_mean(tmp))); % mean of endsamp phase across participants
-%              yy = num2str(cycs);
-%              yy = str2num(yy(end-1:end))/(fs/osc_freq); % calculate how much of a cycle this amount of time is
-%              pred_angle(ang) = wrapToPi(angs(1,ang)+ 2*yy*pi);
-%              R_end(1,ang) = circ_r(tmp,tmp2); % resultant of endsamp across participants
-%                 
-%              angle_start_end(3,ang) = (circ_dist(circ_mean(tmp),pred_angle(ang))); % phase distance from endsamp to predicted angle
-% 
-%              sem_end(1,ang) = circ_std(tmp)./sqrt(length(incl_sub)); % sem of end phase
-%              sem_end(2,ang) = circ_std(circ_dist(angs(1,ang),tmp))./sqrt(length(incl_sub)); % sem of difference between angle at center of bin and end phase
-% 
-%              A = A+1;
-%              
-%              end
-%              
-%              [B I] = sort((angle_start_end(2,:)));
-%              
-%              colorz = linspecer(length(angle_start_end(3,:)));
-% 
-%              subplot(5,4,G)
-%              hold on
-%              errorbar(angle_start_end(1,:),angle_start_end(2,:),sem_end(1,:),'Color','k','LineStyle','none') % angle at center of bin vs mean phase at end, errorbar
-%              scatter(angle_start_end(1,:),pred_angle,2.5,'k','filled') % angle at center of bin vs predicted angle
-%              
-%              for i = 1:length(angle_start_end(3,:))
-% 
-%                 if chan == 2
-% 
-%                     scatter(angle_start_end(1,i),angle_start_end(2,i),25,colorz(i,:),'filled','MarkerEdgeColor','k') % angle at center of bins vs mean at end 
-%                 else
-%                     scatter(angle_start_end(1,:),angle_start_end(2,:),25,colors(2,:))
-%                     errorbar(angle_start_end(1,:),angle_start_end(2,:),sem_end(1,:),'Color',colors(2,:),'LineStyle','none') 
-%                 end
-%              end
-%              
-%              G = G+1;
-% 
-%              ylim([-pi pi])
-%              xlim([-pi pi])
-%              title(['Delay: ',num2str(cycs/fs*1000),' ms'])     
-% 
-%              xlabel('Start Phase')
-%              ylabel('End Phase')
-%              
-%              subplot(5,4,G)
-%              hold on
-%              errorbar(angle_start_end(1,:),angle_start_end(3,:),sem_end(2,:),'Color','k','LineStyle','none') % angle at center of bin vs phase distance from endsamp to predicted angle, errorbar
-%              
-%              for i = 1:length(angle_start_end(3,:))
-%                 if chan == 2
-%                     scatter(angle_start_end(1,i),angle_start_end(3,i),25,colorz(i,:),'filled','MarkerEdgeColor','k')                         
-%                 else
-%                     scatter(angle_start_end(1,:),angle_start_end(3,:),25,colors(2,:))
-%                     errorbar(angle_start_end(1,:),angle_start_end(3,:),sem_end(2,:),'Color',colors(2,:),'LineStyle','none')    
-%                 end
-%              end
-%              
-%             G = G+1;
-%             a = 1;
-%             for A = [0 90 180 270]
-% %             for A = [340 60 160 240]
-%                 xline(wrapToPi(deg2rad(A)),'Color',colors(a,:)) 
-%                 a = a+1;
-%             end
-%              %xlim([-180 180])
-%              xlabel('Start Phase')
-%              ylabel('Phase Change')
-%              yline(0)
-%              text(.95*-pi,-.9*pi,'slower')
-%              %ylim([-200 200])
-%              %ylim([.1 .14])
-%              ylim([-pi pi])
-%              xlim([-pi pi])
-%              text(.95*-pi,.9*pi,'faster')
-%                          
-%                           
-%              subplot(5,4,G)
-%              polarplot([0 circ_mean(angle_start_end(2,:)')],[0 (circ_r(angle_start_end(2,:)',R_end'))],'k','LineWidth',2)
-%              hold on
-% 
-%             if chan == 2
-%                 B = 0.5;
-%                 for i = 1:length(angle_start_end(3,:))
-%                     polarscatter(angle_start_end(2,i),R_end(i),35,colorz(i,:),'filled','MarkerEdgeColor','k')
-%                     B = B+.025;
-%                     hold on
-%                     if G == 3
-%                         before(i) = angle_start_end(2,i);
-%                     else
-%   
-%                         hold on
-%                     end
-%                 end
-%                 
-%             else
-%                 polarplot(angle_start_end(1,:),(circ_dist(circ_mean(angle_start_end(3,:)'),angle_start_end(3,:)))+pi,'Color',colors(2,:))
-% 
-%             end
-%              hold on
-% 
-%              rlim([0 1.1])
-%              G = G+1;
-%              
-%              subplot(5,4,G)
-%              polarplot([0 circ_mean(angle_start_end(3,:)')],[0 (circ_r(angle_start_end(3,:)',R_end'))],'k','LineWidth',2)
-%              hold on
-% 
-%              if chan == 2
-%                 B = .5;
-%                 
-%                 for i = 1:length(angle_start_end(3,:))
-%                     polarscatter(angle_start_end(3,i),R_end(i),35,colorz(i,:),'filled','MarkerEdgeColor','k')
-%                     hold on
-%                     B = B+.025;
-%     
-%                     if G == 3
-%                         before(i) = angle_start_end(2,i);
-%                     else
-%                         hold on
-%                     end
-%                 end
-%              end
-%              
-%              hold on
-%              rlim([0 1.1])
-%              G = G+1;
-%              
-% end
-%              
-%              
-% % saveas(fig,[Savefolder,'Figure6_phase_response_curve_freq',num2str(osc_freq),'_wake_eve_EC.svg']);
-% 
-% %% phase response curve - wake mor (eyes closed)
-% 
-% colors = linspecer(4);
-%              
-% %compute for many onset phases speeding up/slowing down
-%       
-% state = 2; % 1 = eyes open, 2 = eyes closed
-% chan = 2;
-%              
-% clear angle_start_end phase_after R_after angs clear sem_end pred_angle R_end
-% 
-% trigsamp = 501;
-% nbins = 10;
-% fs = 500;
-% osc_freq = 10;
-% 
-% % angs(1,:) = wrapToPi(deg2rad([0:360/10:360]));
-% angs(1,:) = wrapToPi(deg2rad([0:360/nbins:360]));
-% angs(2,:) = angs-pi/nbins;
-% angs(3,:) = angs(1,:)+pi/nbins;
-%              
-% fig = figure('Renderer','painters','units','normalized','outerposition',[0 0 .4 1])
-% hold on
-% % for chan = [12]
-% %                  clear phase_after phase_start_end
-% G = 1;
-% 
-% for cycs = [20/2 100/2 150/2 200/2 250/2]% [10 75 85 95 105];%75 100 125 150]%(25:25:200)
-%     A = 1;
-%     for ang = 1:size(angs,2) %0:18:360;%[36:10:360-36]
-%              
-%          for cond = 1
-%                  clear y
-%                  for s = 1:length(incl_sub)
-%                      
-%                      sub = incl_sub(s);
-%                      
-% %                         tmp = [];
-% %                         for condo = 1:4
-% %                          tmp = [tmp;group_phase.sub{sub}.chan{chan}.cond{condo}];
-% %                         tmp = ERP_trials_allsub_REM{sub}.trial_data_alphafilt_phase(ERP_trials_allsub_REM{sub}.phasic_ndx,:); % trials x samps
-% %                         tmp = ERP_trials_allsub_REM{sub}.trial_data_alphafilt_phase(ERP_trials_allsub_REM{sub}.goodtrial_ndx,:); % trials x samps
-%                           tmp = ERP_trials_allsub_wake_m{sub}.trial_data_alphafilt_phase{state}(ERP_trials_allsub_wake_m{sub}.good_trial_ndx{state},:);
-% %                           tmp = ERP_trials_allsub_wake_e{sub}.trial_data_alphafilt_phase{state}(ERP_trials_allsub_wake_e{sub}.good_trial_ndx{state},:);
-% %                         end
-%                                               
-%                         tmp = tmp(tmp(:,trigsamp)>angs(2,ang)&tmp(:,trigsamp)<angs(3,ang),:);
-%                          
-%                         if ~isempty(tmp)
-%                                       
-%                         Mn = circ_mean(tmp); % average across all trials in that bin
-%                         R = circ_r(tmp);
-%                         to_find = Mn(trigsamp);
-%                      
-%                         cosed = cos(circ_dist(Mn,to_find));
-%                                           
-%                         [pks,locs] = findpeaks(cosed);
-%                                        
-%                      
-%                         phase_before(s,cond) =  Mn(trigsamp);
-%                         cycle_length(s,cond) = abs(trigsamp-locs(find(locs==trigsamp)-1)); % duration from previous peak to stim peak in samps
-%                         phase_after(s,cond,1) = Mn(trigsamp);
-%                         R_after(s,cond,1) = R(trigsamp);
-%                                         
-%                     
-%                         D = (diff(locs(locs>=trigsamp-100 & locs<=trigsamp)));
-%                         %D = D(D>30);
-%                         cycles(s,cond,chan) = mean(D); % mean duration of all cycles between -200 ms and 0
-%                     
-%                     
-%                             for cyc = 1:round(cycs)
-%                                    
-%                                 phase_after(s,cond,cyc+1) = Mn(trigsamp+cyc); % mean phase for all samples from trigsamp to endsamp
-%                                 phase_change(s,cond,cyc) = circ_dist(phase_after(s,cond,cyc),phase_before(s,cond)); % phase distance from each samp after to trigsamp
-%                                 R_after(s,cond,cyc+1) = R(trigsamp+cyc);%  resultant for all samples from trigsamp to endsamp      
-%                         
-%                             end
-%                     
-%                                                                
-%                         else
-%                  
-%                             phase_after(s,cond,:) = NaN(1,1,cycs+1);
-%                             R_after(s,cond,:) = NaN(1,1,cycs+1);
-%                         end
-%             
-%                  end
-%            
-%          end
-%              
-%              
-% %              if cycs == trigsamp
-% %              phase_series(ang,:) = squeeze(circ_mean(phase_after));
-% %              end
-% %              clear cumul_phase_change
-%              
-%              cumul_phase_change = unwrap(squeeze(phase_after(:,1,:))')';
-% 
-%              angle_start_end(1,ang) = (wrapToPi(angs(1,ang))); % angle at center of bin
-%              tmp = phase_after(:,cond,end); % phase of endsamp
-%              tmp2 = R_after(:,cond,end); %  resultant of endsamp
-%              tmp2 = tmp2(~isnan(tmp2));
-%              tmp = tmp(~isnan(tmp));
-%                  
-%              angle_start_end(2,ang) = (wrapToPi(circ_mean(tmp))); % mean of endsamp phase across participants
-%              yy = num2str(cycs);
-%              yy = str2num(yy(end-1:end))/(fs/osc_freq); % calculate how much of a cycle this amount of time is
-%              pred_angle(ang) = wrapToPi(angs(1,ang)+ 2*yy*pi);
-%              R_end(1,ang) = circ_r(tmp,tmp2); % resultant of endsamp across participants
-%                 
-%              angle_start_end(3,ang) = (circ_dist(circ_mean(tmp),pred_angle(ang))); % phase distance from endsamp to predicted angle
-% 
-%              sem_end(1,ang) = circ_std(tmp)./sqrt(length(incl_sub)); % sem of end phase
-%              sem_end(2,ang) = circ_std(circ_dist(angs(1,ang),tmp))./sqrt(length(incl_sub)); % sem of difference between angle at center of bin and end phase
-% 
-%              A = A+1;
-%              
-%              end
-%              
-%              [B I] = sort((angle_start_end(2,:)));
-%              
-%              colorz = linspecer(length(angle_start_end(3,:)));
-% 
-%              subplot(5,4,G)
-%              hold on
-%              errorbar(angle_start_end(1,:),angle_start_end(2,:),sem_end(1,:),'Color','k','LineStyle','none') % angle at center of bin vs mean phase at end, errorbar
-%              scatter(angle_start_end(1,:),pred_angle,2.5,'k','filled') % angle at center of bin vs predicted angle
-%              
-%              for i = 1:length(angle_start_end(3,:))
-% 
-%                 if chan == 2
-% 
-%                     scatter(angle_start_end(1,i),angle_start_end(2,i),25,colorz(i,:),'filled','MarkerEdgeColor','k') % angle at center of bins vs mean at end 
-%                 else
-%                     scatter(angle_start_end(1,:),angle_start_end(2,:),25,colors(2,:))
-%                     errorbar(angle_start_end(1,:),angle_start_end(2,:),sem_end(1,:),'Color',colors(2,:),'LineStyle','none') 
-%                 end
-%              end
-%              
-%              G = G+1;
-% 
-%              ylim([-pi pi])
-%              xlim([-pi pi])
-%              title(['Delay: ',num2str(cycs/fs*1000),' ms'])     
-% 
-%              xlabel('Start Phase')
-%              ylabel('End Phase')
-%              
-%              subplot(5,4,G)
-%              hold on
-%              errorbar(angle_start_end(1,:),angle_start_end(3,:),sem_end(2,:),'Color','k','LineStyle','none') % angle at center of bin vs phase distance from endsamp to predicted angle, errorbar
-%              
-%              for i = 1:length(angle_start_end(3,:))
-%                 if chan == 2
-%                     scatter(angle_start_end(1,i),angle_start_end(3,i),25,colorz(i,:),'filled','MarkerEdgeColor','k')                         
-%                 else
-%                     scatter(angle_start_end(1,:),angle_start_end(3,:),25,colors(2,:))
-%                     errorbar(angle_start_end(1,:),angle_start_end(3,:),sem_end(2,:),'Color',colors(2,:),'LineStyle','none')    
-%                 end
-%              end
-%              
-%             G = G+1;
-%             a = 1;
-%             for A = [0 90 180 270]
-% %             for A = [340 60 160 240]
-%                 xline(wrapToPi(deg2rad(A)),'Color',colors(a,:)) 
-%                 a = a+1;
-%             end
-%              %xlim([-180 180])
-%              xlabel('Start Phase')
-%              ylabel('Phase Change')
-%              yline(0)
-%              text(.95*-pi,-.9*pi,'slower')
-%              %ylim([-200 200])
-%              %ylim([.1 .14])
-%              ylim([-pi pi])
-%              xlim([-pi pi])
-%              text(.95*-pi,.9*pi,'faster')
-%                          
-%                           
-%              subplot(5,4,G)
-%              polarplot([0 circ_mean(angle_start_end(2,:)')],[0 (circ_r(angle_start_end(2,:)',R_end'))],'k','LineWidth',2)
-%              hold on
-% 
-%             if chan == 2
-%                 B = 0.5;
-%                 for i = 1:length(angle_start_end(3,:))
-%                     polarscatter(angle_start_end(2,i),R_end(i),35,colorz(i,:),'filled','MarkerEdgeColor','k')
-%                     B = B+.025;
-%                     hold on
-%                     if G == 3
-%                         before(i) = angle_start_end(2,i);
-%                     else
-%   
-%                         hold on
-%                     end
-%                 end
-%                 
-%             else
-%                 polarplot(angle_start_end(1,:),(circ_dist(circ_mean(angle_start_end(3,:)'),angle_start_end(3,:)))+pi,'Color',colors(2,:))
-% 
-%             end
-%              hold on
-% 
-%              rlim([0 1.1])
-%              G = G+1;
-%              
-%              subplot(5,4,G)
-%              polarplot([0 circ_mean(angle_start_end(3,:)')],[0 (circ_r(angle_start_end(3,:)',R_end'))],'k','LineWidth',2)
-%              hold on
-% 
-%              if chan == 2
-%                 B = .5;
-%                 
-%                 for i = 1:length(angle_start_end(3,:))
-%                     polarscatter(angle_start_end(3,i),R_end(i),35,colorz(i,:),'filled','MarkerEdgeColor','k')
-%                     hold on
-%                     B = B+.025;
-%     
-%                     if G == 3
-%                         before(i) = angle_start_end(2,i);
-%                     else
-%                         hold on
-%                     end
-%                 end
-%              end
-%              
-%              hold on
-%              rlim([0 1.1])
-%              G = G+1;
-%              
-% end
-%              
-%              
-% % saveas(fig,[Savefolder,'Figure6_phase_response_curve_freq',num2str(osc_freq),'_wake_mor_EC.svg']);
-% 
-%              
-%              
-% 
-% 
-%              
-%              
-% 
diff --git a/README.md b/README.md
index 5e7a5cff00518a09aed71f8c8a2328c1c8c9a78d..85d65aa89b91e28baaf7577f17fdf4efbe60b14f 100644
--- a/README.md
+++ b/README.md
@@ -10,10 +10,10 @@ This repository contains scripts for the paper:
 by Valeria Jaramillo, Henry Hebron, Sara Wong, Giuseppe Atzori, Ullrich Bartsch, Derk-Jan Dijk*, Ines R. Violante* (* contributed equally).  
 
 
-Preprint can be found here:
-[biorxiv DOI](https://doi.org/10.1101/2024.03.03.582907)
+Journal article has been published in SLEEP and can be found here:
+[SLEEP DOI](https://doi.org/10.1093/sleep/zsae193)
 
-Raw data and data to create the figures can be found here:
+Raw data, sleep scoring, and data to create the figures can be found here:
 [Zenodo DOI](10.5281/zenodo.10663994)
 
 
diff --git a/analyses/frequency/B_freq_allsub.asv b/analyses/frequency/B_freq_allsub.asv
deleted file mode 100644
index d1bfdfe01cb198b65a3f3b537cc3522631f1f890..0000000000000000000000000000000000000000
--- a/analyses/frequency/B_freq_allsub.asv
+++ /dev/null
@@ -1,297 +0,0 @@
-clear all;
-close all;
-
-addpath(genpath('/users/nemo/software/eeglab'));
-addpath(genpath('/users/nemo/projects/RSN'));
-
-addpath /users/nemo/software/Henry/useful_functions
-
-Folderpath = '/parallel_scratch/nemo/RSN/hdEEG/';
-sub_Folderpath = dir([Folderpath,'RSN*']);
-
-Savefolder = '/parallel_scratch/nemo/RSN/analysis/analysis/frequency_allsub/';
-
-%% Load freq for all subjects and channels
-
-fs = 500;
-
-for s = 1:length(sub_Folderpath)
-        
-    display(sub_Folderpath(s).name);
-
-    freq_file = dir([Folderpath,sub_Folderpath(s).name,filesep,'*_freq.mat']);
-    load([Folderpath,sub_Folderpath(s).name,filesep,freq_file(1).name]);
-    
-    goodREM_file = dir([Folderpath,sub_Folderpath(s).name,filesep,'*_sleep*_fil_czref_goodREM.mat']);
-    load([Folderpath,sub_Folderpath(s).name,filesep,goodREM_file(1).name]);
-    
-    nm_good_file = dir([Folderpath,sub_Folderpath(s).name,filesep,'*_nm_good_psd_allch.mat']);
-    load([Folderpath,sub_Folderpath(s).name,filesep,nm_good_file(1).name]);
-    
-    
-    freq_allch = nanmean(freq.ifq_rem_alpha,1);
-        
-    %%
-    freq.ifq_rem_alpha(:,end:end+1) = NaN;
-    
-    if ~isequal(length(rem_goodsamp2),length(freq.ifq_rem_alpha))
-        error('wrong data length');
-    end
-
-    freq_allnight_alpha = NaN(128,nepochs*epochl*fs);
-    
-    for ch = 1:128
-        freq_allnight_alpha(ch,rem_goodsamp2) = freq.ifq_rem_alpha(ch,:); 
-    end
-    
-    %%
-    
-    freq.ifq_rem_theta(:,end:end+1) = NaN;
-    
-    if ~isequal(length(rem_goodsamp2),length(freq.ifq_rem_theta))
-        error('wrong data length');
-    end
-    
-    
-    freq_allnight_theta = NaN(128,nepochs*epochl*fs);
-    
-    for ch = 1:128
-        freq_allnight_theta(ch,rem_goodsamp2) = freq.ifq_rem_theta(ch,:); 
-    end
-    
-    
-    %% calculate mean frequency for rem, phasic and tonic, first half, second half
-        
-    ifq.rem_alpha_allsub(s,:) = nanmean(freq_allnight_alpha(:,rem_goodsamp2),2);
-    ifq.phasic_alpha_allsub(s,:) = nanmean(freq_allnight_alpha(:,phasic_goodsamp2),2);
-    ifq.tonic_alpha_allsub(s,:) = nanmean(freq_allnight_alpha(:,tonic_goodsamp2),2);
-    
-    ifq.rem_theta_allsub(s,:) = nanmean(freq_allnight_theta(:,rem_goodsamp2),2);
-    ifq.phasic_theta_allsub(s,:) = nanmean(freq_allnight_theta(:,phasic_goodsamp2),2);
-    ifq.tonic_theta_allsub(s,:) = nanmean(freq_allnight_theta(:,tonic_goodsamp2),2);
-    
-    firsthalf_samp = rem_goodsamp2(1:floor(length(rem_goodsamp2)/2)-1);
-    secondhalf_samp = rem_goodsamp2(floor(length(rem_goodsamp2)/2):end);
-    
-    firsthalf_phasic_samp = intersect(firsthalf_samp,phasic_goodsamp2);
-    secondhalf_phasic_samp = intersect(secondhalf_samp,phasic_goodsamp2);
-    
-    firsthalf_tonic_samp = intersect(firsthalf_samp,tonic_goodsamp2);
-    secondhalf_tonic_samp = intersect(secondhalf_samp,tonic_goodsamp2);
-    
-    ifq.rem_firsthalf_alpha_allsub(s,:) = nanmean(freq_allnight_alpha(:,firsthalf_samp),2);
-    ifq.rem_secondhalf_alpha_allsub(s,:) = nanmean(freq_allnight_alpha(:,secondhalf_samp),2);
-    
-    ifq.rem_firsthalf_phasic_alpha_allsub(s,:) = nanmean(freq_allnight_alpha(:,firsthalf_phasic_samp),2);
-    ifq.rem_secondhalf_phasic_alpha_allsub(s,:) = nanmean(freq_allnight_alpha(:,secondhalf_phasic_samp),2);
-    
-    ifq.rem_firsthalf_tonic_alpha_allsub(s,:) = nanmean(freq_allnight_alpha(:,firsthalf_tonic_samp),2);
-    ifq.rem_secondhalf_tonic_alpha_allsub(s,:) = nanmean(freq_allnight_alpha(:,secondhalf_tonic_samp),2);
-    
-    ifq.rem_firsthalf_theta_allsub(s,:) = nanmean(freq_allnight_theta(:,firsthalf_samp),2);
-    ifq.rem_secondhalf_theta_allsub(s,:) = nanmean(freq_allnight_theta(:,secondhalf_samp),2);
-    
-    ifq.rem_firsthalf_phasic_theta_allsub(s,:) = nanmean(freq_allnight_theta(:,firsthalf_phasic_samp),2);
-    ifq.rem_secondhalf_phasic_theta_allsub(s,:) = nanmean(freq_allnight_theta(:,secondhalf_phasic_samp),2);
-    
-    ifq.rem_firsthalf_tonic_theta_allsub(s,:) = nanmean(freq_allnight_theta(:,firsthalf_tonic_samp),2);
-    ifq.rem_secondhalf_tonic_theta_allsub(s,:) = nanmean(freq_allnight_theta(:,secondhalf_tonic_samp),2);
-    
-    
-      %% calculate phasic/tonic trial ndx
-      
-     for con = 1:8
-         
-     ON_start_all = nm.ON_start_good{con}; 
-          
-     art_trial = nm.art_trialndx{con};
-     art_trialndx = find(art_trial == 1);
-     good_trialndx = setdiff(1:length(ON_start_all),art_trialndx);
-     nm.good_trialndx{con} = good_trialndx;
-
-%         art_trialndx = find(nm.art_trialndx{con} == 1);
-%         art_trialndx_new = find(nm.art_trialndx{con} == 1);
-%         good_trialndx = nm.good_trialndx{con};
-%         good_trialndx_new = nm.good_trialndx_new{con};
-        
-     phasic_trial = zeros(length(ON_start_all),1);
-     tonic_trial = zeros(length(ON_start_all),1);
-%      
-%      
-     for t = 1:length(ON_start_all)
-         
-         if ~ismember(t,art_trialndx)
-         
-         dt = 6;
-         trial_samps = ON_start_all(t)-dt*fs:ON_start_all(t)+2*dt*fs-1;
-         ifq_trials_alpha(t,:,:) = freq_allnight_alpha(:,trial_samps);
-         ifq_trials_theta(t,:,:) = freq_allnight_theta(:,trial_samps);
-
-            if length(intersect(trial_samps,tonic_goodsamp2)) == 2*dt*fs
-                tonic_trial(t) = 1;
-            elseif length(intersect(trial_samps,phasic_goodsamp2)) > 0
-                phasic_trial(t) = 1;   
-            end
-                        
-         end
-          
-     end
-%      
-%      phasic_goodtrialndx = intersect(find(phasic_trial == 1),good_trialndx);
-%      tonic_goodtrialndx = intersect(find(tonic_trial == 1),good_trialndx);
-%      
-% %      nm.phasic_goodtrialndx{con} = phasic_goodtrialndx;
-% %      nm.tonic_goodtrialndx{con} = tonic_goodtrialndx;     
-% %      nm.ntrials_phasic{con} = length(phasic_goodtrialndx);
-% %      nm.ntrials_tonic{con} = length(tonic_goodtrialndx);
-% %      
-% %      if ~ nm.ntrials_phasic{con} + nm.ntrials_tonic{con} + length(art_trialndx) == length(ON_start_all)
-% %         warning('trials have not been assigned to tonic/phasic/artefact') ;
-% %      end
-%      
-%      clear tonic_trial phasic_trial
-        
-    %% calculate mean across all good, phasic and tonic trials
-            
-        if ~isempty(good_trialndx)
-            mifq_trials_alpha = squeeze(nanmean(ifq_trials_alpha(good_trialndx,:,:),1));            
-            ifq.mifq_trials_alpha(s,con,:,:) = mifq_trials_alpha;
-            ifq.off_ifq_alpha(s,con,:) = nanmean(mifq_trials_alpha(:,1:6*fs),2);
-            ifq.on_ifq_alpha(s,con,:) = nanmean(mifq_trials_alpha(:,6*fs+1:12*fs),2);
-        else
-            ifq.mifq_trials_alpha(s,con,:,:) = NaN(128,9000);
-            ifq.off_ifq_alpha(s,con,:) = NaN(128,1);
-            ifq.on_ifq_alpha(s,con,:) = NaN(128,1); 
-        end
-        
-        if ~isempty(phasic_goodtrialndx)         
-            mifq_trials_alpha_phasic = squeeze(nanmean(ifq_trials_alpha(phasic_goodtrialndx,:,:),1));
-            ifq.mifq_trials_alpha_phasic(s,con,:,:) = mifq_trials_alpha_phasic;
-            ifq.off_ifq_alpha_phasic(s,con,:) = nanmean(mifq_trials_alpha_phasic(:,1:6*fs),2);
-            ifq.on_ifq_alpha_phasic(s,con,:) = nanmean(mifq_trials_alpha_phasic(:,6*fs+1:12*fs),2);
-        else
-            ifq.mifq_trials_alpha_phasic(s,con,:,:) = NaN(128,9000);
-            ifq.off_ifq_alpha_phasic(s,con,:) = NaN(128,1);
-            ifq.on_ifq_alpha_phasic(s,con,:) = NaN(128,1); 
-        end
-        
-        if ~isempty(tonic_goodtrialndx)  
-            mifq_trials_alpha_tonic = squeeze(nanmean(ifq_trials_alpha(tonic_goodtrialndx,:,:),1));
-            ifq.mifq_trials_alpha_tonic(s,con,:,:) = mifq_trials_alpha_tonic;
-            ifq.off_ifq_alpha_tonic(s,con,:) = nanmean(mifq_trials_alpha_tonic(:,1:6*fs),2);
-            ifq.on_ifq_alpha_tonic(s,con,:) = nanmean(mifq_trials_alpha_tonic(:,6*fs+1:12*fs),2);
-        else
-            ifq.mifq_trials_alpha_tonic(s,con,:,:) = NaN(128,9000);
-            ifq.off_ifq_alpha_tonic(s,con,:) = NaN(128,1);
-            ifq.on_ifq_alpha_tonic(s,con,:) = NaN(128,1); 
-        end
-        
-        
-        
-         if ~isempty(good_trialndx)
-            mifq_trials_theta = squeeze(nanmean(ifq_trials_theta(good_trialndx,:,:),1));            
-            ifq.mifq_trials_theta(s,con,:,:) = mifq_trials_theta;
-            ifq.off_ifq_theta(s,con,:) = nanmean(mifq_trials_theta(:,1:6*fs),2);
-            ifq.on_ifq_theta(s,con,:) = nanmean(mifq_trials_theta(:,6*fs+1:12*fs),2);
-        else
-            ifq.mifq_trials_theta(s,con,:,:) = NaN(128,9000);
-            ifq.off_ifq_theta(s,con,:) = NaN(128,1);
-            ifq.on_ifq_theta(s,con,:) = NaN(128,1); 
-        end
-        
-        if ~isempty(phasic_goodtrialndx)         
-            mifq_trials_theta_phasic = squeeze(nanmean(ifq_trials_theta(phasic_goodtrialndx,:,:),1));
-            ifq.mifq_trials_theta_phasic(s,con,:,:) = mifq_trials_theta_phasic;
-            ifq.off_ifq_theta_phasic(s,con,:) = nanmean(mifq_trials_theta_phasic(:,1:6*fs),2);
-            ifq.on_ifq_theta_phasic(s,con,:) = nanmean(mifq_trials_theta_phasic(:,6*fs+1:12*fs),2);
-        else
-            ifq.mifq_trials_theta_phasic(s,con,:,:) = NaN(128,9000);
-            ifq.off_ifq_theta_phasic(s,con,:) = NaN(128,1);
-            ifq.on_ifq_theta_phasic(s,con,:) = NaN(128,1); 
-        end
-        
-        if ~isempty(tonic_goodtrialndx)  
-            mifq_trials_theta_tonic = squeeze(nanmean(ifq_trials_theta(tonic_goodtrialndx,:,:),1));
-            ifq.mifq_trials_theta_tonic(s,con,:,:) = mifq_trials_theta_tonic;
-            ifq.off_ifq_theta_tonic(s,con,:) = nanmean(mifq_trials_theta_tonic(:,1:6*fs),2);
-            ifq.on_ifq_theta_tonic(s,con,:) = nanmean(mifq_trials_theta_tonic(:,6*fs+1:12*fs),2);
-        else
-            ifq.mifq_trials_theta_tonic(s,con,:,:) = NaN(128,9000);
-            ifq.off_ifq_theta_tonic(s,con,:) = NaN(128,1);
-            ifq.on_ifq_theta_tonic(s,con,:) = NaN(128,1); 
-        end
-        
-        clear mifq_trials_alpha mifq_trials_alpha_phasic mifq_trials_alpha_tonic mifq_trials_theta mifq_trials_theta_phasic mifq_trials_theta_tonic
-     
-    
-     end
-     
-%      save([Folderpath,sub_Folderpath(s).name,filesep,nm_good_file(1).name],'nm','-append','-v7.3');
-    clear nm
-
- %% calculate mean frequency for wake eve, eyes open, and eyes closed
-  
-     
-%     nm_good_file = dir([Folderpath,sub_Folderpath(s).name,filesep,'*_wake_e*_nm_good.mat']);
-%     load([Folderpath,sub_Folderpath(s).name,filesep,nm_good_file(1).name]); 
-%     
-%     goodwake_file = dir([Folderpath,sub_Folderpath(s).name,filesep,'*_wake_e_fil_czref_goodwake.mat']);
-%     load([Folderpath,sub_Folderpath(s).name,filesep,goodwake_file(1).name]);
-%        
-%     freq.ifq_wake_e_alpha(:,end:end+1) = NaN;
-%     freq.ifq_wake_e_theta(:,end:end+1) = NaN;
-%     
-%     freq_allwake_e_alpha = freq.ifq_wake_e_alpha;
-%     freq_allwake_e_theta = freq.ifq_wake_e_theta;
-%    
-%     eye_closure_samp = nm.startsamp_ERP{2}(1);
-%     EO_samp = 1:eye_closure_samp;
-%     EC_samp = eye_closure_samp+1:size(freq_allwake_e_alpha,2);
-%     EO_e_goodsamp2 = intersect(EO_samp,wake_goodsamp2);
-%     EC_e_goodsamp2 = intersect(EC_samp,wake_goodsamp2);
-% 
-%     ifq.wake_e_alpha_allsub(s,:) = nanmean(freq_allwake_e_alpha(:,wake_goodsamp2),2);
-%     ifq.wake_e_EO_alpha_allsub(s,:) = nanmean(freq_allwake_e_alpha(:,EO_e_goodsamp2),2);
-%     ifq.wake_e_EC_alpha_allsub(s,:) = nanmean(freq_allwake_e_alpha(:,EC_e_goodsamp2),2);
-%     
-%     ifq.wake_e_theta_allsub(s,:) = nanmean(freq_allwake_e_theta(:,wake_goodsamp2),2);
-%     ifq.wake_e_EO_theta_allsub(s,:) = nanmean(freq_allwake_e_theta(:,EO_e_goodsamp2),2);
-%     ifq.wake_e_EC_theta_allsub(s,:) = nanmean(freq_allwake_e_theta(:,EC_e_goodsamp2),2);
-%     
-%    clear nm
-   
- %% calculate mean frequency for wake mor, eyes open, and eyes closed
-  
-     
-%     nm_good_file = dir([Folderpath,sub_Folderpath(s).name,filesep,'*_wake_m*_nm_good.mat']);
-%     load([Folderpath,sub_Folderpath(s).name,filesep,nm_good_file(1).name]); 
-%     
-%     goodwake_file = dir([Folderpath,sub_Folderpath(s).name,filesep,'*_wake_m_fil_czref_goodwake.mat']);
-%     load([Folderpath,sub_Folderpath(s).name,filesep,goodwake_file(1).name]);
-%        
-%     freq.ifq_wake_m_alpha(:,end:end+1) = NaN;
-%     freq.ifq_wake_m_theta(:,end:end+1) = NaN;
-%     
-%     freq_allwake_m_alpha = freq.ifq_wake_m_alpha;
-%     freq_allwake_m_theta = freq.ifq_wake_m_theta;
-%    
-%     eye_closure_samp = nm.startsamp_ERP{2}(1);
-%     EO_samp = 1:eye_closure_samp;
-%     EC_samp = eye_closure_samp+1:size(freq_allwake_m_alpha,2);
-%     EO_m_goodsamp2 = intersect(EO_samp,wake_goodsamp2);
-%     EC_m_goodsamp2 = intersect(EC_samp,wake_goodsamp2);
-% 
-%     ifq.wake_m_alpha_allsub(s,:) = nanmean(freq_allwake_m_alpha(:,wake_goodsamp2),2);
-%     ifq.wake_m_EO_alpha_allsub(s,:) = nanmean(freq_allwake_m_alpha(:,EO_m_goodsamp2),2);
-%     ifq.wake_m_EC_alpha_allsub(s,:) = nanmean(freq_allwake_m_alpha(:,EC_m_goodsamp2),2);
-%     
-%     ifq.wake_m_theta_allsub(s,:) = nanmean(freq_allwake_m_theta(:,wake_goodsamp2),2);
-%     ifq.wake_m_EO_theta_allsub(s,:) = nanmean(freq_allwake_m_theta(:,EO_m_goodsamp2),2);
-%     ifq.wake_m_EC_theta_allsub(s,:) = nanmean(freq_allwake_m_theta(:,EC_m_goodsamp2),2);    
-    
-      
-end
-
-
-save([Savefolder,'freqalphatheta_allsub_',date,'.mat'],'ifq','-v7.3');
-
diff --git a/analyses/resultant/A_Phase_sub_allch.sub b/analyses/resultant/A_Phase_sub_allch.sub
deleted file mode 100644
index fe0e9699b33ed0f593fee27dd95ecff574168cfc..0000000000000000000000000000000000000000
--- a/analyses/resultant/A_Phase_sub_allch.sub
+++ /dev/null
@@ -1,30 +0,0 @@
-#!/bin/bash
-
-#SBATCH --job-name=array
-#SBATCH --array=1-19         #Array job indices/range for $SLURM_ARRAY_TASK_ID (can be incremented if desired)
-#SBATCH --time=72:00:00
-#SBATCH --partition=high_mem
-#SBATCH --ntasks=1
-#SBATCH --mem=4G
-#SBATCH --error=array_%A_%a.err   #Error file label by job number and index.
-
-# Print the task id.
-echo "My SLURM_ARRAY_TASK_ID: " $SLURM_ARRAY_TASK_ID > test_"$SLURM_ARRAY_TASK_ID"
-
-
-cd $SLURM_SUBMIT_DIR
-
-module use /opt/software/user-pkgs/modulefiles/all
-module use /opt/software/user-pkgs/modulefiles/psychology
-module load setcondaenv
-
-conda info -e
-
-# module load amica
-module load matlab/2021a
-
-matlab -nosplash -nodesktop -nodisplay -r "run('/users/nemo/projects/RSN/git/RSN/analyses/resultant/A_Phase_sub_allch($SLURM_ARRAY_TASK_ID)');exit"
-
-
-
-