%
%
%Solution to Assignment 7, part 2, exercise 2e, second 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);

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

%repeat, but with a much lower cut-off frequency stimulus
%effective cut-off frequency of 50 Hz
new_stim_vect = normrnd(0,1,1,l_tvect/20);

new_stim_vect2 = resample(new_stim_vect,20,1);
new_stim_lgn = conv(new_stim_vect2,wt_lgn);
new_stim_lgn = new_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
new_stim_lgn_min = abs(min(new_stim_lgn));

%compute scale factor so that lowest firing rate change is -50 spk/sec 
new_scale_fact = 50/new_stim_lgn_min;
new_stim_lgn2 = new_scale_fact*new_stim_lgn; 

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

%compute the corresponding Poisson spike train
[new_spk_vect,new_v_vect] = inhom_pp(dt,new_stim_lgn2);

%find spike times
new_spk_inds = find(new_spk_vect > 0.5);

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

for i = 1:length(new_spk_inds)
  if ( (new_spk_inds(i) > l_lgn) )
    new_est_lgn(1:l_lgn) = new_est_lgn + new_stim_vect2(new_spk_inds(i):-1:new_spk_inds(i)- l_lgn+1);
  end;
end;

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

%plot with original filter
figure(2);
plot(t_lgn,new_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)');