function [coeff,lags]=cch(ia,ib,T,resolution)
% Calculating the cross-correlation histogram of two spike interval sequence
% with N bins and a time span of T.
% Inputs:
% ia and ib: the two input sequence of firing intervals, in msec
% T: the maximum time lag, in msec
% resolution: the time resolution of the histogram, in msec
% Outputs:
% corr: cross-correlation
% lags: time lags
% coeff: normalized cross-correlation
%
%Written by Ning Jiang, Institute of Biomedical Engineering, Univesity of New
%Brunswick, NB, Canada, E3B 5A3
%Email: ning.jiang@unb.ca
%
%Date: Mar 5, 2007
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
eta=1/resolution;
corr=0;
reference=0;
referenee=ia(1);
ia=round(ia/resolution);
ib=round(ib/resolution);
%creating referenee train, an increasing sequence representing the timing of the reference unit.
for l=1:length(ia)-1
temp=referenee(l)+ia(l+1);
referenee=[referenee,temp];
end
%finding firings in reference train with that is in (-T,T) range of a firing in referenee train
for n=1:length(ib)
reference=reference+ib(n);
temp=referenee-reference;
corr=[corr,temp(find(abs(temp)<T))];
%corr=[corr,temp(find(abs(temp)<T*eta))];
end
corr=corr(2:length(corr));
%generate the bin edges of the histogram
lags=(-T*eta+1:T*eta)/eta;
%get histogram and get correlation. coounts equals the correlation calculated by Function :xcorr(ia,ib,T), i.e. not normalized correlation.
counts=histc(corr,lags);
lags=lags(1:length(lags)-1);
counts=counts(1:length(counts)-1);
corr=counts*eta;%compensate for eta
% length of recording, in msec
t=referenee(length(referenee));
%getting statistics of a firing train of 1's (a firing) and 0's (time between firings)
%mean of the two firing trains
mean_a=length(ia)/(t*eta);
mean_b=length(ib)/(t*eta);
%variance of the two firing trains
var_a=((1-mean_a)^2*length(ia)+mean_a^2*(t-length(ia)))/(t*eta);
var_b=((1-mean_b)^2*length(ib)+mean_b^2*(reference-length(ib)))/(reference*eta);
crossvar=sqrt(var_a*var_b);
%cova of the two firing trains: normalized correlation minus product of the
%two means
cova=counts/(t*eta)-mean_a*mean_b;
%correlation coefficient
coeff=cova/crossvar;