%% Does the actual reliability calculations. Ruff removed HRTF's have less 
% directions, so downsample normal files. Makes things go faster too.

clear 
%update these paths
cd('F:\MATLAB\')
addpath('F:\MATLAB\Brian_Reliability\Scripts')
% % Load HRIRs

owl_folder{1}='Merry';
owl_folder{2}='Owl19';
owl_folder{3}='Owl39';
owl_folder{4}='Ugly';
owl_folder{5}='Vespucci';

haircut_id{1}='normal';
haircut_id{2}='ohne'; % maximum cut
date_folder = '06-06-22 Data'; %change to current date if rerunning next section
%% Full analysis of reliability. Takes a while... 
% Once this is done once, this and the next section can be skipped. Could
% add another if loop to automatically skip this section.

tic
for o=2:5
    for h=1:2
        clearvars -except h haircut_id o owl_folder date_folder
        load(['HRTF_Wagner_Lab\results\' owl_folder{o} '\' owl_folder{o} '__' haircut_id{h} '__hrir_result.mat'])
        
        newazimuths = -160:20:160;
        newelevations = -60:20:80;
        
        dir=NaN(2,length(azimuths)*length(elevations));
        TF1=NaN(length(azimuths)*length(elevations),size(hrir_l,3));
        TF2=NaN(length(azimuths)*length(elevations),size(hrir_l,3));
        count=0;
        for a=1:length(azimuths)
            for e=1:length(elevations)
                count=count+1;
                dir(:,count)=[elevations(e);azimuths(a)];
                TF1(count,:)=hrir_l(a,e,:);
                TF2(count,:)=hrir_r(a,e,:);
            end
        end        
                
        newdir=NaN(2,length(newazimuths)*length(newelevations));
        count=0;
        for a=1:length(newazimuths)
            for e=1:length(newelevations)
                count=count+1;
                newdir(:,count)=[newelevations(e);newazimuths(a)];
            end
        end                
        
        dir = dir';
        newdir = newdir';
        downsample = ismember(dir, newdir, 'rows');
        TF1 = TF1(downsample, :);
        TF2 = TF2(downsample, :);
        dir = newdir';
        azimuth = newazimuths;
        elevation = newelevations;
        
        clear newdir newazimuth newelevation
        
        % % Analysis parameters
        
        %center frequencies of cochlear filters (Hz)
        cF = 400:200:10000;
        %cF = linspace(1000,9000, 9);
        n_cF=length(cF);
        
        %Sampling frequency for HRTFs (Hz)
        Fs = samplingrate;
        
        %Factor by which you will upsample
        Factor = 5;
        %Number of trials for each condition
        Nt=5;
        % % Get cochlear filters
        %Upsampled frequency
        FS=Fs*Factor;
        %Get filters
        fcoefs=getFilterCoefs(cF,FS);
        
        
        % % Number of directions and initialize        
        nd_tot=size(TF1,1);
        n_dir = size(TF1,1);
        
        IPDm_all = zeros(n_dir, n_cF, Nt);
        IPDs_all = zeros(n_dir, n_cF, Nt);
%         ITDm=zeros(n_dir,n_cF);
%         ITDs=zeros(n_dir,n_cF);
        ILDm_all = zeros(n_dir, n_cF, Nt);
        ILDs_all = zeros(n_dir, n_cF, Nt);
        
        gainr_all = zeros(n_dir, n_cF, Nt);
        gainl_all = zeros(n_dir, n_cF, Nt);
        
        
        % %
        tic
        
        TF1_res = resample(TF1',Factor,1)';
        TF2_res = resample(TF2',Factor,1)';
        
        parfor idir = 1:n_dir % TODO parfor
            
            %Resample the HRIRs
            hL=TF1_res(idir,:);
            hR=TF2_res(idir,:);
            %Run over trials
            
            for tt=1:Nt
                
                %Generate a target sound
                s=genSignal(100,1/30,1,2,2*pi*12);
                s1=resample(s,Factor,1);
                
                %Generate a masker
                s=genSignal(100,1/30,1,2,2*pi*12);
                s2=resample(s,Factor,1);
                
                %Convolve target with HRIR
                sL1=conv(hL,s1);
                sR1=conv(hR,s1);
                
                Z=zeros(n_dir,n_cF);
                GL=zeros(n_dir,n_cF);
                GR=zeros(n_dir,n_cF);
                ITDp=zeros(n_cF,n_dir);
                
                %Run over masking directions
                %disp(num2str([toc o h idir tt]))
                
                for k=1:nd_tot
                    %disp(num2str([o h idir tt k]))
                    
                    %Resample the HRIRs
                    hL2=TF1_res(k,:); %#ok<PFBNS>
                    hR2=TF2_res(k,:); %#ok<PFBNS>
                    
                    %Convolve masker with HRIR and add to target
                    sL=sL1+conv(hL2,s2);
                    sR=sR1+conv(hR2,s2);
                    
                    %cochlear (gammmatone) filterbank
                    vL=ERBFilterBank(sL,fcoefs);
                    vR=ERBFilterBank(sR,fcoefs);
                    
                    %Run over frequency and compute cross correlation
                    for icF = 1:n_cF
                        [x,l]=xcorr(vL(icF,:),vR(icF,:));
                        [~,j]=max(x);
                        ITDp(icF,k)=l(j)*1/Fs*1e6/Factor;
                    end
                    
                    %Compute ILD
                    R=rms(vR,2);
                    L=rms(vL,2);
                    Z(k,:)=20*log10(R./L);
                    GR(k,:)=20*log10(R);
                    GL(k,:)=20*log10(L);
                    
                end
                
                
                %Run over frequency and compute mean and SD
                for n=1:n_cF
                 
                    IPDp=ITDp(n,:)*cF(n)/1e6*2*pi; %#ok<PFBNS>
                    IPDs_all(idir,n,tt) = circstd(IPDp/Nt);
                    IPDm_all(idir,n,tt) = circmean(IPDp/Nt);
                    
                end
                
                %compute mean and SD for ILD. 
                ILDm_all(idir,:,tt) = mean(Z,1);
                ILDs_all(idir,:,tt) = std(Z,[],1);
                
                %compute gain
                gainr_all(idir,:,tt) = mean(GR,1);
                gainl_all(idir,:,tt) = mean(GL,1);

                
            end     
        end
        IPDs = sum(IPDs_all, 3);
        IPDm = sum(IPDm_all, 3);
        ITDs = (1e6 * IPDs) ./ (2 * pi * cF);
        ITDm = (1e6 * IPDm) ./ (2 * pi * cF);
        ILDs = sum(ILDs_all, 3)/Nt;
        ILDm = sum(ILDm_all, 3)/Nt;
        gainl = sum(gainl_all, 3)/Nt;
        gainr = sum(gainr_all, 3)/Nt;

        toc
        save(['Shadron_Pena_2022\', date_folder, '\', owl_folder{o} '__' haircut_id{h} '__hrtf_fig1.mat'])
    end
end
%% Some minor editing

for o = 1:5
    for h = 1:2
        clearvars -except h haircut_id o owl_folder date_folder
        load(['Shadron_Pena_2022\', date_folder, '\', owl_folder{o} '__' haircut_id{h} '__hrtf_fig1.mat'])        
        
        %have to wrap ITD's to make things make sense
        %Earlier had to convert 17x8x563 to a 136x563, so now have to revert
        %back to 17x8x20.
        
        x = size(hrir_l, 3);
        y = size(cF, 2);
            
            IxD.IPDm = zeros(17,8,y);
            IxD.IPDs = zeros(17,8,y);
            IxD.ITDm = zeros(17,8,y);
            IxD.ITDs = zeros(17,8,y);
            IxD.ILDm = zeros(17,8,y);
            IxD.ILDs = zeros(17,8,y);
            IxD.gainl = zeros(17,8,y);
            IxD.gainr = zeros(17,8,y);
            
            for n = 1:y
                for a = 1:17
                    for e = 1:8
                        counter = (a-1)*8 + e;
                        IxD.IPDm(a,e,n) = IPDm(counter,n);
                        IxD.IPDs(a,e,n) = IPDs(counter,n);
                        IxD.ITDm(a,e,n) = ITDm(counter,n);
                        IxD.ITDs(a,e,n) = ITDs(counter,n);
                        IxD.ILDm(a,e,n) = ILDm(counter,n);
                        IxD.ILDs(a,e,n) = ILDs(counter,n);
                        IxD.gainl(a,e,n) = gainl(counter,n);
                        IxD.gainr(a,e,n) = gainr(counter,n);
                    end
                end
            end

        clear n a e count Factor Fs FS idir n n_cF n_dir nd_tot Nt tt x ILDm ILDs IPDm IPDp IPDs ITDp ITDs ITDm
        clear samplingrate s1 s2 sL1 sR1 TF1 TF2 Z ax counter fcoef dir s fcoefs hL hR hrir_l hrir_r 
        
        if size(IxD.IPDm,2) > 8 %Loop occurs if spatial locations are > 17 azimuth and > 9 elevation positions
            %Resamples to have locations that are
            %-160:20:160 in azimuth and -60:20:80 in elevation
            for n = 1:y
                for a = 1:2:34
                    for e = 2:2:16
                        IxD2.IPDm((a+1)/2, e/2, n) = IxD.IPDm(a,e,n);
                        IxD2.IPDs((a+1)/2, e/2, n) = IxD.IPDs(a,e,n);
                        IxD2.ITDm((a+1)/2, e/2, n) = IxD.ITDm(a,e,n);
                        IxD2.ITDs((a+1)/2, e/2, n) = IxD.ITDs(a,e,n);
                        IxD2.ILDm((a+1)/2, e/2, n) = IxD.ILDm(a,e,n);
                        IxD2.ILDs((a+1)/2, e/2, n) = IxD.ILDs(a,e,n);
                    end
                end
            end
            
            
            clear IxD
            IxD = IxD2;
            clear IxD2
            
            azimuths = -160:20:160;
            elevations = -60:20:80;
        end
        
        
        save(['Shadron_Pena_2022\', date_folder, '\', owl_folder{o} '__' haircut_id{h} '__stats_fig1.mat'])
        
    end
end
%% Compile all data into one struct
clearvars -except h haircut_id o owl_folder date_folder

condition = [];
for o = [1 2 3 4 5]
    for h = 1:2
        clearvars -except h haircut_id o owl_folder condition date_folder
        load(['Shadron_Pena_2022\', date_folder, '\', owl_folder{o} '__' haircut_id{h} '__stats_fig1.mat'])
        condition.(owl_folder{o}).(haircut_id{h}) = IxD;
         
    end
end
            
clear IxD name

%average the owls
for h=1:2
for o = [1 2 3 4 5]
        if o == 1
            condition.(haircut_id{h}).avg.IPDm = condition.(owl_folder{o}).(haircut_id{h}).IPDm;
            condition.(haircut_id{h}).avg.IPDs = condition.(owl_folder{o}).(haircut_id{h}).IPDs;
            condition.(haircut_id{h}).avg.ITDm = condition.(owl_folder{o}).(haircut_id{h}).ITDm;
            condition.(haircut_id{h}).avg.ITDs = condition.(owl_folder{o}).(haircut_id{h}).ITDs;
            condition.(haircut_id{h}).avg.ILDm = condition.(owl_folder{o}).(haircut_id{h}).ILDm;
            condition.(haircut_id{h}).avg.ILDs = condition.(owl_folder{o}).(haircut_id{h}).ILDs;
            condition.(haircut_id{h}).avg.gainl = condition.(owl_folder{o}).(haircut_id{h}).gainl;
            condition.(haircut_id{h}).avg.gainr = condition.(owl_folder{o}).(haircut_id{h}).gainr;

        else
            condition.(haircut_id{h}).avg.IPDs = condition.(haircut_id{h}).avg.IPDs + condition.(owl_folder{o}).(haircut_id{h}).IPDs;
            condition.(haircut_id{h}).avg.ITDm = condition.(haircut_id{h}).avg.ITDm + condition.(owl_folder{o}).(haircut_id{h}).ITDm;
            condition.(haircut_id{h}).avg.ITDs = condition.(haircut_id{h}).avg.ITDs + condition.(owl_folder{o}).(haircut_id{h}).ITDs;
            condition.(haircut_id{h}).avg.ILDm = condition.(haircut_id{h}).avg.ILDm + condition.(owl_folder{o}).(haircut_id{h}).ILDm;
            condition.(haircut_id{h}).avg.ILDs = condition.(haircut_id{h}).avg.ILDs + condition.(owl_folder{o}).(haircut_id{h}).ILDs;
            condition.(haircut_id{h}).avg.gainl = condition.(haircut_id{h}).avg.gainl + condition.(owl_folder{o}).(haircut_id{h}).gainl;
            condition.(haircut_id{h}).avg.gainr = condition.(haircut_id{h}).avg.gainr + condition.(owl_folder{o}).(haircut_id{h}).gainr;
        end
end
condition.(haircut_id{h}).avg.IPDm = condition.(haircut_id{h}).avg.IPDm / 5;
condition.(haircut_id{h}).avg.IPDs = condition.(haircut_id{h}).avg.IPDs / 5;
condition.(haircut_id{h}).avg.ITDm = condition.(haircut_id{h}).avg.ITDm / 5;
condition.(haircut_id{h}).avg.ITDs = condition.(haircut_id{h}).avg.ITDs / 5;
condition.(haircut_id{h}).avg.ILDm = condition.(haircut_id{h}).avg.ILDm / 5;
condition.(haircut_id{h}).avg.ILDs = condition.(haircut_id{h}).avg.ILDs / 5;
condition.(haircut_id{h}).avg.gainl = condition.(haircut_id{h}).avg.gainl/5;
condition.(haircut_id{h}).avg.gainr = condition.(haircut_id{h}).avg.gainr/5;
end

%% Makes the reliability maps across azimuth as in Cazettes et al (2014)

lowf = 1200; % The lowest frequency you want to include in the normalization
highf = 8000;
%low = ((lowf-1000)/200); 
low = find(cF == lowf) - 1;
high = find(cF == highf);

for h = 1:2
    condition.(haircut_id{h}).avg.reliability.IPD = (condition.(haircut_id{h}).avg.IPDs.^-1);
    condition.(haircut_id{h}).avg.reliability.ILD = (condition.(haircut_id{h}).avg.ILDs.^-1);
    condition.(haircut_id{h}).avg.reliability.normalizedIPD = zeros(size(azimuths,1),high-low);
    condition.(haircut_id{h}).avg.reliability.normalizedILD = zeros(size(azimuths,1),high-low);
    condition.(haircut_id{h}).avg.reliability.normIPDaz = zeros(size(azimuths,1),2);

    
    for a = 1:17
        [big, pt] = max(condition.(haircut_id{h}).avg.reliability.IPD(a,4,1+low:high));
        for n = 1:high-low
            condition.(haircut_id{h}).avg.reliability.normalizedIPD(a,n) = condition.(haircut_id{h}).avg.reliability.IPD(a,4,n+low)/big;
        end
        condition.(haircut_id{h}).avg.reliability.normIPDaz(a,1) = mean(condition.(haircut_id{h}).avg.reliability.IPD(a,4,1+low:high).^-1);
        condition.(haircut_id{h}).avg.reliability.normIPDaz(a,2) = std(condition.(haircut_id{h}).avg.reliability.IPD(a,4,1+low:high).^-1);
    end
    for a = 1:17
        [big, pt] = max(condition.(haircut_id{h}).avg.reliability.ILD(a,4,1+low:high));
        for n = 1:high-low
            condition.(haircut_id{h}).avg.reliability.normalizedILD(a,n) = condition.(haircut_id{h}).avg.reliability.ILD(a,4,n+low)/big;
        end
    end

end

short = cF(1+low:high);

% Can run to see overall pattern before normalization
% figure
% normalmap = imagesc(azimuths, short, condition.(haircut_id{1}).avg.reliability.normalizedIPD'); axis xy
% title(['IPD Reliability Along Azimuth: ', haircut_id{1}, ' (n = 5)']); colorbar; colormap jet
% 
% figure
% ohnemap = imagesc(azimuths, short, condition.(haircut_id{2}).avg.reliability.normalizedIPD'); axis xy
% title(['IPD Reliability Along Azimuth: ', haircut_id{2}, ' (n = 5)']); colorbar; colormap jet

%% Figure 1a and 1b

posazi = 0:20:160;
shortKHz = short/1000;

colors = f_gsafecmap(256);
colormap(flipud(colors));

d = figure;
set(gcf,'Position',[100 100 600 450]);
posnormalmap = imagesc(posazi, shortKHz, condition.(haircut_id{1}).avg.reliability.normalizedIPD(9:17,:)'); axis xy
title('Normal')
%title(['IPD Reliability Along Azimuth: ', haircut_id{1}, ' (n = 5)']); 
posnormalmapc = colorbar; colormap(flipud(colors)); clim([0.2 1]);
%posnormalmapc.Label.String = 'Reliability'; posnormalmapc.Label.FontSize = 18;
ylabel('Frequency (KHz)'); xlabel('Azimuth (°)');
ylabel('Frequency (KHz)', 'FontSize', 22); xlabel('Azimuth (°)', 'FontSize', 22);
set(gca,'XTick',[0 45 90 135 180], 'XTickLabel', [0 45 90 135 180], 'fontsize', 18);
set(d,'Units','Inches');
set(d,'PaperPositionMode','Auto','PaperUnits','Inches','PaperSize',[2.8, 2.24])
text(-35,8.6, 'a', 'FontSize', 28);
%saveas(posnormalmap, 'F:\Various Files\My Papers\Shadron and Pena 2022\Matlab\normreliab.png')

figure
set(gcf,'Position',[100 100 600 450]);
posohnemap = imagesc(posazi, shortKHz, condition.(haircut_id{2}).avg.reliability.normalizedIPD(9:17,:)'); axis xy
title('Ruffcut')
%title(['IPD Reliability Along Azimuth: ', haircut_id{2}, ' (n = 5)']); 
posohnemapc = colorbar; colormap(flipud(colors)); clim([0.2 1]);
posohnemapc.Label.String = 'Reliability'; posohnemapc.Label.FontSize = 22;
%ylabel('Frequency (KHz)'); xlabel('Azimuth (deg)');
xlabel('Azimuth (°)', 'FontSize', 22);
set(gca,'XTick',[0 45 90 135 180], 'XTickLabel', [0 45 90 135 180], 'fontsize', 18);
text(-35,8.6, 'b', 'FontSize', 28);
%saveas(posohnemap, 'F:\Various Files\My Papers\Shadron and Pena 2022\Matlab\ruffcutreliab.png')


%For non normalized:
%imagesc(posazi, shortKHz, permute(squeeze(condition.(haircut_id{1}).avg.reliability.IPD(9:17,4,:)),[2 1]));colorbar
%% Figure 1c

diffIPD = condition.(haircut_id{2}).avg.reliability.IPD - condition.(haircut_id{1}).avg.reliability.IPD;
diffILD = condition.(haircut_id{2}).avg.reliability.ILD - condition.(haircut_id{1}).avg.reliability.ILD;
dIPD = squeeze(diffIPD(:,4,:));
dILD = squeeze(diffILD(:,4,:));

figure
set(gcf,'Position',[100 100 600 450]);
IPDd = imagesc(azimuths(9:end), cF(low:high), dIPD(9:end,low:high)'); axis xy; 
IPDdc = colorbar; clim([-.3 .3]); colors = f_gsafecmap(256); colormap(flipud(colors));
ylabel('Frequency (KHz)', 'FontSize', 22); 
xlabel('Azimuth (°)', 'FontSize', 22);
IPDdc.Label.String = 'Ruffcut IPD - Normal IPD';
IPDdc.Label.FontSize = 22;
IPDdc.FontSize = 16;
set(gca,'XTick',[0 45 90 135 180], 'XTickLabel', [0 45 90 135 180], 'fontsize', 18);
set(gca,'YTick',2000:2000:8000, 'YTickLabel', [2 4 6 8], 'fontsize', 18);
text(-35,8600, 'c', 'FontSize', 24);
%saveas(IPDd, 'F:\Various Files\My Papers\Shadron and Pena 2022\Matlab\reliabdiff.png')

%set(IPDdc,'Units','Inches');

% colors = f_gsafecmap(256);
% colormap(flipud(colors));
% ILDd = imagesc(azimuths(9:end), cF(4:end), dILD(9:end,4:end)'); axis xy; colorbar; caxis([-.5 .5]);
% ylabel('Frequency (KHz)','FontSize', 20); 
% xlabel('Azimuth (Deg)','FontSize', 20);
% title('Difference in ILD reliability', 'FontSize', 20)

%% ILD, not shown in manuscript

posazi = 160:-20:0;
shortKHz = short/1000;

colors = f_gsafecmap(256);
colormap(flipud(colors));

figure
negnormalmap = imagesc(posazi, shortKHz, condition.(haircut_id{1}).avg.reliability.normalizedILD(1:9,:)'); axis xy
title('Normal', 'FontSize', 20); 
negnormalmapc = colorbar;  clim([0.2 1]);
negnormalmapc.Label.String = 'Reliability'; negnormalmapc.Label.FontSize = 18;
colormap(flipud(colors));
ylabel('Frequency (KHz)', 'FontSize', 20); xlabel('Azimuth (deg)', 'FontSize', 20);
set(gca,'XTick',[0 45 90 135 180], 'XTickLabel', [0 45 90 135 180], 'fontsize', 16);

figure
negohnemap = imagesc(posazi, shortKHz, condition.(haircut_id{2}).avg.reliability.normalizedILD(1:9,:)'); axis xy
title('Ruff Cut', 'FontSize', 20); 
negohnemapc = colorbar; colormap(flipud(colors)); clim([0.2 1]);
negohnemapc.Label.String = 'Reliability'; negohnemapc.Label.FontSize = 18;
%ylabel('Frequency (KHz)','FontSize', 20); 
xlabel('Azimuth (deg)', 'FontSize', 20);
set(gca,'XTick',[0 45 90 135 180], 'XTickLabel', [0 45 90 135 180], 'fontsize', 16);

%For non normalized:
%imagesc(posazi, shortKHz, permute(squeeze(condition.(haircut_id{1}).avg.reliability.IPD(9:17,4,:)),[2 1]));colorbar

%% ILD reliability, not used in manuscript
colors = f_gsafecmap(256);
colormap(flipud(colors));

figure
posnormalmap = imagesc(posazi, shortKHz, condition.(haircut_id{1}).avg.reliability.normalizedILD(9:17,:)'); axis xy
title('Normal: ILD reliability')
%title(['ILD Reliability Along Azimuth: ', haircut_id{1}, ' (n = 5)']); 
posnormalmapc = colorbar; colormap(flipud(colors));
posnormalmapc.Label.String = 'Reliability'; posnormalmapc.Label.FontSize = 12;
ylabel('Frequency (KHz)'); xlabel('ITD (us)');

figure
posohnemap = imagesc(posazi, shortKHz, condition.(haircut_id{2}).avg.reliability.normalizedILD(9:17,:)'); axis xy
title('Ruffcut: ILD reliability')
%title(['IPD Reliability Along Azimuth: ', haircut_id{2}, ' (n = 5)']); 
posohnemapc = colorbar; colormap(flipud(colors));
posohnemapc.Label.String = 'Reliability'; posohnemapc.Label.FontSize = 12;
ylabel('Frequency (KHz)'); xlabel('ITD (us)');

%For non normalized:
%imagesc(posazi, shortKHz, permute(squeeze(condition.(haircut_id{1}).avg.reliability.IPD(9:17,4,:)),[2 1]));colorbar

%% s.d. IPD averaged across freq for each azimuth, not used in manuscript

IPDstd = figure;
hold on
scatter(azimuths, condition.(haircut_id{1}).avg.reliability.normIPDaz(:,1)/180, 40, 'blue', 'filled', 'o')
scatter(azimuths, condition.(haircut_id{2}).avg.reliability.normIPDaz(:,1)/180, 40, 'red', 'filled', 'o')


errorbar(azimuths, condition.(haircut_id{1}).avg.reliability.normIPDaz(:,1)/180, condition.(haircut_id{1}).avg.reliability.normIPDaz(:,2)/180, ...
     'blue', 'LineStyle', 'none')
errorbar(azimuths, condition.(haircut_id{2}).avg.reliability.normIPDaz(:,1)/180, condition.(haircut_id{2}).avg.reliability.normIPDaz(:,2)/180, ...
     'red', 'LineStyle', 'none')
xlabel('Azimuth (deg)')
ylabel('IPD standard deviation')
xlim([-10,90])

legend('Normal', 'Ruff-Removed', 'Location', 'southeast')
fontsize(gca, 16,"points")