%
%
%Solution to Assignment 7, part 2, exercise 2c&d, first reverse correlation
%
%F. Gabbiani 11/25/08
%

dt = 1; %time step in msec
tmax = 25000; %in msec
t_vect = 0:dt:tmax-dt;
l_tvect = length(t_vect);

stim_vect = normrnd(0,1,1,l_tvect);

[wt_lgn, t_lgn] = da(5.5,350,dt);
l_lgn = length(wt_lgn);

stim_lgn = conv(stim_vect,wt_lgn);
stim_lgn = stim_lgn(1:l_tvect);


%rescale the output so as to have a mean of 50 spk/sec 
%with a minimal value of 0 spk/sec

%compute minimum value
stim_lgn_min = abs(min(stim_lgn));

%compute scale factor so that lowest firing rate change is -50 spk/sec 
scale_fact = 50/stim_lgn_min;
stim_lgn2 = scale_fact*stim_lgn; 

%add 50 spk/sec to obtain a mean of 50 spk/sec
stim_lgn2 = stim_lgn2 + 50;

%compute the corresponding Poisson spike train
[spk_vect,v_vect] = inhom_pp(dt,stim_lgn2);

%find spike times
spk_inds = find(spk_vect > 0.5);

%compute reverse correlation over the 
%length of the filter
est_lgn = zeros(1,l_lgn);

for i = 1:length(spk_inds)
  if ( (spk_inds(i) > l_lgn) )
    est_lgn(1:l_lgn) = est_lgn + stim_vect(spk_inds(i):-1:spk_inds(i)- l_lgn+1);
  end;
end;

%normalize peak to 1
est_lgn = est_lgn/max(est_lgn);

%plot with original filter
figure(1);
plot(t_lgn,est_lgn);
hold on;
plot(t_lgn,wt_lgn,'r');
title('Comparison of original and reconstructed LGN filter');
xlabel('time (ms)');
ylabel('firing rate (arbitrary units)');