%This script visualizes the EDA specific information for all 9 subjects to
%replicate Figs 2-4 and figures in the S2 Appendix
%Figures 1-9 are for subjects 1-9 respectively

load('Data/EDA_deepdive.mat');

subjs = 1:9;

for s = 1:length(subjs)
    %Identify indices
    t = (1:size(EDA_temp_amp(s).X,1))*20/60;
    mu_amp = EDA_temp_amp(s).X(:,4);
    sigma_amp = EDA_temp_amp(s).X(:,5);
    mu_BR = EDA_temp_amp(s).X(:,2);
    sigma_BR = EDA_temp_amp(s).X(:,3);
    
    %Smooth signals before computing correlation
    mu_amp_smoothed = smooth(mu_amp,9,'rlowess');
    sigma_amp_smoothed = smooth(sigma_amp,9,'rlowess');
    mu_BR_smoothed = smooth(mu_BR,9,'rlowess');
    sigma_BR_smoothed = smooth(sigma_BR,9,'rlowess');
    
    %Replace NaNs
    mu_amp_smoothed(isnan(mu_amp_smoothed)) = nanmean(mu_amp_smoothed);
    sigma_amp_smoothed(isnan(sigma_amp_smoothed)) = nanmean(sigma_amp_smoothed);
    mu_BR_smoothed(isnan(mu_BR_smoothed)) = nanmean(mu_BR_smoothed);
    sigma_BR_smoothed(isnan(sigma_BR_smoothed)) = nanmean(sigma_BR_smoothed);
    
    mu_amp(isnan(mu_amp)) = nanmean(mu_amp);
    sigma_amp(isnan(sigma_amp)) = nanmean(sigma_amp);
    mu_BR(isnan(mu_BR)) = nanmean(mu_BR);
    sigma_BR(isnan(sigma_BR)) = nanmean(sigma_BR);

    %Compute correlation
    [r_mean,p_mean] = movcorr(mu_amp_smoothed,mu_BR_smoothed,30);
    [r_std,p_std] = movcorr(sigma_amp_smoothed,sigma_BR_smoothed,30);
     
    %Plot
    figure;

    subplot(6,1,1)
    plot(t, mu_BR,'LineWidth',1)
    ylabel('Pulses/min')
    title('Mean pulse rate')
    xlabel('Time (min)')
    ax1 = gca;
    
    subplot(6,1,2)
    plot(t, mu_amp,'LineWidth',1)
    ylabel('Pulse prominence') 
    title('Mean pulse amplitude')
    xlabel('Time (min)')
    ax2 = gca;
    
    subplot(6,1,3)
    plot(t, r_mean,'LineWidth',2)
    ylabel('Correlation') 
    title('Correlation between Mean Pulse Rate and Amplitude')
    xlabel('Time (min)')
    ax3 = gca;

    subplot(6,1,4)
    plot(t, sigma_BR,'LineWidth',1)
    ylabel('Pulses/min')
    title('Standard Deviation of Pulse Rate')
    xlabel('Time (min)')
    ax4 = gca;
    
    subplot(6,1,5)
    plot(t, sigma_amp,'LineWidth',1)
    ylabel('Pulse prominence') 
    title('Standard Deviation of Pulse Amplitude')
    xlabel('Time (min)')
    ax5 = gca;
    
    subplot(6,1,6)
    plot(t, r_std,'LineWidth',2)
    ylabel('Correlation')
    title('Correlation between Standard Deviation of Pulse Rate and Amplitude')
    xlabel('Time (min)')
    ax6 = gca;
    
    linkaxes([ax1,ax2,ax3,ax4,ax5,ax6],'x')
end