%
% Solution to exercise 1, parts 3 and 4, homework 2, part 2
% CAAM 415, 10/23/08
%

rho = 40; %in units of spk/sec

mu = 1000/40; %ISI in msec

%draw 1000 exponentially distributed intervals
r = exprnd(mu,1,1000);

%add a refractory period to each interval
r_ref = r + 2;
mean_isi = mean(r_ref);

v1 = r_ref(2:end) - mean_isi;
v2 = r_ref(1:end-1) - mean_isi;

%compute the serial correlation coefficient from the 
%formulas in the lecture notes
scc1 = mean(v1.*v2)/(std(v1)*std(v2));

info_str = sprintf('serial correlation coefficient first spike train %.2f',scc1);
disp(info_str);

%Generate the modified ISI as given in the exercises 
r_ref_mod = zeros(1,1000);
r_ref_mod(1) = r_ref(1);
for i=2:1000
    if ( r_ref(i-1) < mean_isi)
        r_ref_mod(i) = r_ref(i) + 5;
    else
        r_ref_mod(i) = r_ref(i) - 5;
    end;
end;

%compute serial correlation  coefficient of modified ISI sequence
mean_isi_mod = mean(r_ref_mod);
v1_mod = r_ref_mod(2:end) - mean_isi_mod;
v2_mod = r_ref_mod(1:end-1) - mean_isi_mod;

scc1_mod = mean(v1_mod.*v2_mod)/(std(v1_mod)*std(v2_mod));

info_str = sprintf('serial correlation coefficient second spike train: %.2f',scc1_mod);
disp(info_str);