AuditoryHierarchy_DataSetAndCode_NotUploaed/Code/Compute_similarity_matrix.m

function [Corr_WiNs,Corr] = Compute_similarity_matrix(AmpResp,AllSamples)

%% This function calculates the similarity matrix which estimates the correlation between pairs
%% of sounds after correction for trial to trial noise
%% INPUT
% AmpResp = the responses to be correlations : ROI * stim * trials
% AllSamples = 1 (default) calculates all possible pairs of trials, 0
% instead does this for a random 50 combiations which is a lot quicker

%% OUTPUT
% Corr : matrix containing the correlation matrix between trial pairs for all
% combinations : permutations * sounds * sounds
% Corr : similarity matrix : sounds * sounds

%% intialize
if not(exist('AllSamples'))
    AllSamples = 1;
end

% clean up NaN cells
NaN_cells = find(isnan(sum(squeeze(nanmean(AmpResp,3)'))));
AmpResp(NaN_cells,:,:) = [];
disp([num2str(length(NaN_cells)) ' ROIs out of ' num2str(size(AmpResp,1)) 'excluded because they had NaNs'])

% Not all data sets have all of the trials, get the maximum number of
% trials shared across datasets 
maxgoodtrials = find(squeeze(sum(sum(isnan(AmpResp),1),2))==0,1,'last');
disp([num2str(maxgoodtrials) ' trials without NaNs'])

% Get the trial pairs
AllTrialPairs = combnk(1:maxgoodtrials,2);

if AllSamples
    % all the samples
    TrialPairsToUse = AllTrialPairs;
else
    % a subset of trial pairs
    TrialPairsToUse = AllTrialPairs(randperm(size(AllTrialPairs,1),30),:);
end

%% Calculate correlation for all pairs
for perm = 1:size(TrialPairsToUse,1)
    Corr_WiNs(perm,:,:) = corr((squeeze(AmpResp(:,:,TrialPairsToUse(perm,1)))),(squeeze(AmpResp(:,:,TrialPairsToUse(perm,2)))));
end

%% Calculate noise corrected correlation
for k = 1:size(Corr_WiNs,2)
    for kk = 1:size(Corr_WiNs,2)
        if k == kk
            Corr(k,kk) = NaN;
        else
            XCorr = abs(squeeze(nanmean(Corr_WiNs(:,k,k),1)));
            YCorr = abs(squeeze(nanmean(Corr_WiNs(:,kk,kk),1)));
            XYCorr = squeeze(nanmean(Corr_WiNs(:,kk,k),1));
            Corr(kk,k) =  XYCorr./(sqrt(XCorr*YCorr));
        end
    end
end



end