function p=scrambledpvalue(A,B)
% scrambledpvalue calculates the negative log of the p value when 
% correlating two input signals A and B. A is a two-dimensional signal 
% (e.g. gamma band power over all channels) and B is a one-dimensional
% signal (e.g. speech power). The p-values are computed using a
% randomization approach, where A is flipped, circularly shifted by a 
% random amount and correlated with B for a large number of simulations
% (default being 500). Then the distribution of r-values found from
% these simulations are fit to a beta distribution to find the
% parameters of this random distribution. The p-value is then found by 
% finding where the true r-value from correlating A with B lies on this
% random distribution.

% CALLING SEQUENCE:
%       scambledpvalue(A,B)
% INPUT:
%       A: Two-dimensional signal (e.g. gamma band across all ECoG channels) 
%       that we want to correlate with B
%       B: One-dimensional signal (e.g. speech) or two-dimensional signal
%       (e.g. covert gamma band power) that we want to correlate with A.

% Specify the number of simulations. Default is 500. 
numsim=100;

% Flip A and circularly shift it by a random number of samples. Repeat this
% for the number of simulations defined above. Populate a matrix 'data'
% with the scaled r-values for each iteration.
d=flipud(A);
data=zeros(numsim,size(A,2));
for i=1:numsim
   r=randperm(size(A,1));
   g=circshift(d,r(1));
   parfor j=1:size(A,2)
       if size(B,2)==1
         data(i,j)=(corr(g(:,j),B(:,1))/2)+0.5;
       else
         data(i,j)=(corr(g(:,j),B(:,j))/2)+0.5;
       end
   end
end

%Set the nan values to zero.
ind=isnan(data);
index= ind(1,:)==1;
data(:,index)=0;

% Fit a beta distribution to the random distribution of scaled r-values
% obtained above. Find p-value using the true scaled r-value from 
% correlating A and B and fitting it to this random beta distribution.
p=zeros(1,size(A,2));
for i=1:size(data,2)
   if length(unique(data(:,i)))~=1
     param=mle(data(:,i),'distribution','beta');
     if size(B,2)==1
        p(1,i)=1-betacdf((corr(A(:,i),B(:,1))/2)+0.5,param(1),param(2));
     else
         p(1,i)=1-betacdf((corr(A(:,i),B(:,i))/2)+0.5,param(1),param(2));
     end
    else
      p(1,i)=1;
    end
end

% Add an epsilon value (10^-20) to the p-values to get rid of zeros prior
% to finding the -log of the p-values. These -log p values are then 
% returned to the calling function/script.
p=p';
p=p+(10^-20);
p=-log(p);
end