%
%
% Comparison of theoretical and numerical computation
% of the LIF transfer function
%
% Fabrizio Gabbiani 11/03/08
%

%sampling step in the time domain
dt_ms = 1; %in msec
dt_s = dt_ms*1e-3; %in sec

%time constant of the LIF in msec
tau_ms = 30; %in msec
tau_s = tau_ms*1e-3; %in sec;

%Nyquist frequency 
f_nyquist = 1/(2*dt_s); %in Hz

%sampling vector in the time domain in msec
p_2 = 9; %pick the closest power of 2 
t_vect = 0:dt_ms:2^p_2-1;
n_samples = length(t_vect);

%compute the theoretical transfer function in the time domain
g_t = exp(-t_vect/tau_ms);

figure(3); 
plot(t_vect,g_t);
xlabel('time (ms)');
ylabel('transfer function');

%sampling step in the frequency domain
f_samp = 1/(n_samples*dt_s); %in Hz

%corresponding frequency vector in Hz
f_vect = 0:f_samp:f_nyquist;
n_f = length(f_vect);

%circular frequency vector
omega_vect = 2*pi*f_vect;

%theoretical real and imaginary part of the transfer function in the
%frequency domain
g_hat_r =  tau_s./(1+ (tau_s*omega_vect).^2);
g_hat_i = -(tau_s^2*omega_vect)./(1 + (tau_s*omega_vect).^2); 

%compute the real and imaginary parts by discrete FFT
g_fft = fft(g_t);
g_fft_r = real(g_fft(1:n_f));
g_fft_i = imag(g_fft(1:n_f));

%convert back discrete FFT to continuous approximation
g_fft_rs = g_fft_r*dt_s;
g_fft_is = g_fft_i*dt_s;

figure(1);
semilogx(f_vect,g_hat_r);
hold on;
semilogx(f_vect,g_fft_rs,'r');
xlabel('frequency [Hz]');
ylabel('real part');

figure(2);
semilogx(f_vect,g_hat_i);
hold on;
semilogx(f_vect,g_fft_is,'r');
xlabel('frequency [Hz]');
ylabel('imaginary part');