%
% Solution exercise 2, part d
% Fifth exercise serie, course part 2
% 
% CAAM 415, 11/13/08
%

%parameters of the receptive field
%center surround extent
r_c = 0.24; %in deg
r_s = 0.96;

%relative weight of center vs. surround
k_s = 0.06;
k_c = 1;

%start by plotting in the frequency domain 
%since we can determine the 
%appropriate scaling factor of the response there

df = 0.01; %sampling step in Hz
f_vect = df:df:2.5;

%circular frequency
cf_vect = 2*pi*f_vect;

%Fourier transform according to the lecture notes
vf_vect = k_c*r_c^2*pi*exp(-r_c^2*cf_vect.^2/4) - ...
    k_s*r_s^2*pi*exp(-r_s^2*cf_vect.^2/4);   

%rescale peak value to 50 as in the lecture note figure
scf = 50/max(vf_vect);
vf_vect = scf*vf_vect;

%plot in log log coordinates
figure(3);
loglog(f_vect,vf_vect);
title('Fourier transform of spatial receptive field');
xlabel('Spatial frequency c/deg');
ylabel('Contrast sensitivity');
xlim([0.01 10]);
ylim([1 100]);

%sampling step dx below chosen to give 
%approx +/- 3 deg around center
%of ref with 512 points in each side
dx = 0.006;
x_vect = (-1023:1:1024)*dx;

%receptive field according to lecture notes
vx_vect = k_c*r_c*sqrt(pi)*exp(-(x_vect/r_c).^2) - ...
    k_s*r_s*sqrt(pi)*exp(-(x_vect/r_s).^2);

%scale appropriately to obtain the change in firing rate
%in spk/sec
vx_vect = scf * vx_vect;

%sinusoidal grating spatial and temporal frequencies
f_x = [0.1 1 2]; %in cycles/deg

f_t = 3;%in cyles/sec

%time vector
t_vect = 0:999;
l_t = length(t_vect);

for j= 1:3
  %generate storage space for the grating
  gr = zeros(1000,2048);

  for i = 1:l_t
    %generate grating at each time point
    gr(i,1:2048) = cos(2*pi*f_x(j)*x_vect - 2*pi*f_t*t_vect(i)*1e-3);
  end;

  %compute change in firing rate response
  r = gr*vx_vect'*dx;
  
  %store peak firing rate
  r_x(j) = max(r);
end;

%plot on theoretical transfer function to verify
%result
figure(3);
hold on;
plot(f_x, r_x,'or');