%
% Solution exercise 1, homework 7, part 2
% 
% CAAM 415, 11/25/08
%

%parameters here are k_x, s_x, k_t, s_t, corresponding ...
%to the peak spatial (c/deg) and temporal (c/sec) values of the
%filter and their spread. 

%circular frequency corresponding to 4.2 cycles/deg
k_x = 2*pi*4.2;

%in degrees
s_x = 0.1; 

dx = 0.004;            %units are degrees
nx = 128;

x = ((-nx/2+1)*dx:dx:nx/2*dx);  %units are degrees, 128 points, column vector

%compute the spatial filters
r_es = g_es(k_x,s_x,x);
r_os = g_os(k_x,s_x,x);

%circular frequency corresponding to 8 cycles/sec
k_t = 2*pi*8;

%in sec
s_t = 0.031;

%temporal filter
dt = 0.002;            %units are seconds
nt = 128;
t = ((-nt/2+1)*dt:dt:nt/2*dt)';  %units are seconds

t_et = f_es(k_t,s_t,t);
%t_es = zeros(128,1);
%t_es(64,1) = 1;

t_ot = f_os(k_t,s_t,t);
%t_os = zeros(128,1);
%t_os(64,1) = 1;

%grating time vector
t_gr = (0:dt:2)'; %row vector
n_tgr =length(t_gr);

%grating parameters
eta_x = 2*pi*4; %4 cycles/deg
eta_t = 2*pi*8; %8 cycles/sec

%
%grating moving to the right
%

% first dim is time (rows); second dim is space (columns)
gr = cos(eta_x*repmat(x,[n_tgr 1]) - eta_t*repmat(t_gr,[1 nx]));

%
%grating moving to the left
%
gl = cos(eta_x*repmat(x,[n_tgr 1]) + eta_t*repmat(t_gr,[1 nx]));

%%%%%%%%%%%%%%%%%%%%%%%%%
%
% separable responses
%
%%%%%%%%%%%%%%%%%%%%%%%%%

%
%even filter response right moving 
%

%spatial component of response; even spatial filter
gr_es = gr*r_es';

%temporal convolution
intgr_es_et = conv(gr_es,t_et);

%correct for time shift using conv
gr_es_et(1:n_tgr) = intgr_es_et(64:n_tgr+64-1)*dt;

%
%odd filter response right moving
%

%spatial component; odd spatial filter
gr_os = gr*r_os'; 

%temporal convolution
intgr_os_ot = conv(gr_os,t_ot);

%correct for time shift using conv
gr_os_ot(1:n_tgr) = intgr_os_ot(64:n_tgr+64-1)*dt;

%
%complex cell response; right drifting grating
%
c_sep_r = gr_es_et.^2 + gr_os_ot.^2;


%
%even filter response left moving
%

%spatial component
gl_es = gl*r_es';

%temporal convolution
intgl_es_et = conv(gl_es,t_et);

%correct for time shift using conv
gl_es_et(1:n_tgr) = intgl_es_et(64:n_tgr+64-1)*dt;

%
%odd filter response left moving
%

%spatial component
gl_os = gl*r_os';

%temporal convolution
intgl_os_ot = conv(gl_os,t_ot);

%correct for time shift using conv
gl_os_ot(1:n_tgr) = intgl_os_ot(64:n_tgr+64-1)*dt;

%
% complex cell response; left drifting grating
%
c_sep_l = gl_es_et.^2 + gl_os_ot.^2;

figure;
h1 = plot(t_gr,c_sep_r);
hold on;
h2 = plot(t_gr,c_sep_l,'r.');

xlabel('time (s)');
ylabel('response (arbitrary units)');
title('Response of ME model based on separable RFs to left/right moving gratings'); 

%demonstrate frequency doubling of the response
mean_r = mean(c_sep_r);
max_r = max(c_sep_r);
cos_f2 = mean_r + (max_r-mean_r)*cos(2*eta_t*t_gr);
hc = plot(t_gr,cos_f2,'k.');

legend([h1 h2 hc],{'right moving' 'left' '2F cosine'});
%%%%% NOTE %%%%%%%%
%These responses are oscillatory because the convolution with t_et and t_ot introduces an additional 
%90 degrees phase shift. This results in a total 180 phase shift and frequency doubling of the response. The
%following responses, based on an even temporal filter combined with even and odd spatial filters, are not oscillatory. 
%Similarly, one can generate equivalent non-oscillatory responses using the odd temporal filter combined with odd/even 
%spatial ones or a fixed spatial filter with a combination of odd/even temporal ones. 
%%%%%%%%%%%%%%%%%%%

%temporal convolution
intgr_os_et = conv(gr_os,t_et);

%correct for time shift using conv
gr_os_et(1:n_tgr) = intgr_os_et(64:n_tgr+64-1)*dt;

%
%complex cell response; right drifting grating
%
d_sep_r = gr_es_et.^2 + gr_os_et.^2;

%temporal convolution
intgl_os_et = conv(gl_os,t_et);

%correct for time shift using conv
gl_os_et(1:n_tgr) = intgl_os_et(64:n_tgr+64-1)*dt;

%
%complex cell response; right drifting grating
%
d_sep_l = gl_es_et.^2 + gl_os_et.^2;

figure;
h1a = plot(t_gr,d_sep_r);
hold on;
h2a = plot(t_gr,d_sep_l,'r.');
legend([h1a h2a],{'right moving' 'left'});
xlabel('time (s)');
ylabel('response (arbitrary units)');
title('Response of ME model based on separable RFs to left/right moving gratings'); 

%%%%%%%%%%%%%%%%%%%%%%%%%
%
% non-separable responses
%
%%%%%%%%%%%%%%%%%%%%%%%%%

%response of first oriented (non-separable,even) filter to a 
%right moving grating
gr_ns_e_st = gr_es_et + gr_os_ot;

%second oriented filter

%temporal convolution
intgr_os_et = conv(gr_os,t_et);

%correct for time shift using conv
gr_os_et(1:n_tgr) = intgr_os_et(64:n_tgr+64-1)*dt;


%temporal convolution
intgr_es_ot = conv(gr_es,t_ot);

%correct for time shift using conv
gr_es_ot(1:n_tgr) = intgr_es_ot(64:n_tgr+64-1)*dt;

gr_ns_o_st = gr_os_et - gr_es_ot;

c_nsep_r = gr_ns_e_st.^2 + gr_ns_o_st.^2;


%response of first oriented filter to a 
%left moving grating
gl_ns_e_st = gl_es_et + gl_os_ot;

%second oriented filter

%temporal convolution
intgl_os_et = conv(gl_os,t_et);

%correct for time shift using conv
gl_os_et(1:n_tgr) = intgl_os_et(64:n_tgr+64-1)*dt;


%temporal convolution
intgl_es_ot = conv(gl_es,t_ot);

%correct for time shift using conv
gl_es_ot(1:n_tgr) = intgl_es_ot(64:n_tgr+64-1)*dt;

gl_ns_o_st = gl_os_et - gl_es_ot;

c_nsep_l = gl_ns_e_st.^2 + gl_ns_o_st.^2;

figure;
h3 = plot(t_gr,c_nsep_r);
hold on;
h4 = plot(t_gr,c_nsep_l,'r.');
legend([h3 h4],{'right moving' 'left'});
xlabel('time (s)');
ylabel('response (arbitrary units)');
title('Response of ME model based on non-separable RFs to left/right moving grating');