%
% Solution exercise 2, part a, b, c 
% Fifth exercise serie of part 2
% 
% CAAM 415, 11/13/08
%

%parameters of the RGC 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(1);
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;

%plot receptive field
figure(2);
plot(x_vect,vx_vect);
title('Spatial receptive field');
xlabel('Spatial angle (deg)');
ylabel('Firing frequency (rel. spont. spk/s)');

%store receptive field in wrap-around order
vx2_vect(1:1025) = vx_vect(1024:2048);
vx2_vect(2048:-1:1026) = vx_vect(1023:-1:1);

%fast fourier transform and scale appropriately
vf2_vect = fft(vx2_vect)*dx;

%compute sampling step in the frequency domain
df = 1/(2048*dx);
f2_vect = (0:1:1025)*df;

%compare with theoretical result
figure(1);
hold on;
loglog(f2_vect,real(vf2_vect(1:1026)),'r');