% Computational Neuroscience II: Foundations of Neural Coding, SS 2003
% Solutions for Assignment 2 (22.04.03) 
% (c) written by Richard Kempter 


% The name of this script is 'blatt2.m'.It can be executed by typing
% 'blatt2' (without the suffix .m) at a matlab prompt.

clear all; % clear all variables from memory
           % type "help clear" for more info
% -------------- Plotting with Matlab
t = 0:0.1:6.2; % generate a vector
y = sin(t);    % generate a second vector

figure(1)      % create figure 1 window (next 'plot' will be in this window)
plot(t,y);     % plot vector t versus vector y
xlabel('t')    % always label your axes 
ylabel('sin(t)')       
title('Sine function') % never forget to add a title

% -------------- Logarithm function 
t = 0.01:0.001:5;
f = log(t);
figure(2);   % we open a second figure window
plot(t,f);
xlabel('t');
ylabel('log(t)');
title('Logarithm function');

% ------------- Exponential function
tau=2;           % time constant in units of milliseconds
t=0:0.001:10;    % time vector  
f1=exp(-t/tau); % exponential function
f2=t/(tau^2).*f1;  % alpha function (note the use of '.*') 
figure(3);
plot(t,f1,'r-.');  % plot t versus f1 using a red, dot-dashed line
hold on;           % the next plot command will not erase the previous
                   % graph  
plot(t,f2,'b-','LineWidth',3); % thick blue full line 
xlabel('time (ms)')
legend('Exponential function','Alpha function')
hold off           

% -------------------------
% Simple statistics of neural spike trains

s1 = load('spiketrain1.dat'); % load spike trains, spike times are in units
s2 = load('spiketrain2.dat'); % of milliseconds

n1 = length(s1);  % get the number of spikes  
n2 = length(s2);

T1 = (max(s1) - min(s1))/1000;  % time interval in units of seconds
T2 = (max(s2) - min(s1))/1000;

f1 = (n1-1)/T1;   % firing rates in units of Hz
f2 = (n2-1)/T2;   % ('-1' because the last spike has no associated interval)

%print the results
fprintf('Spike train 1: %d spikes / %f sec = %f Hz (firing rate)\n',n1,T1,f1);
fprintf('Spike train 2: %d spikes / %f sec = %f Hz (firing rate)\n',n2,T2,f2);

% ----------------------
% now we measure the firing rate of the neuron in 
% time bins of size T

T = 500; % in units of milliseconds 
	 % (try different numbers here, e.g. 10, 100, 1000, 10000 ...)

% first we create the EDGES vector for histc
edges1 = min(s1):T:max(s1);   
edges2 =   s2(1):T:s2(n2); % two different methods doing the same thing  
% try also: 'edges2 = s2(1):T:s2(end)' What does 'end' mean here? 

N1 = histc(s1,edges1);   % generate histogram counts     
N2 = histc(s2,edges2);   % if 'histc' is not available, use
                         % 'hist(s2,length(edges2))' instead

% convert numbers of spikes into firing rates in units of Hz  
F1=N1/T * 1000; % factor 1000 for units of Hz
F2=N2/T * 1000;

% ------------ plot the time-binned histograms
figure(4)
bar(edges1/1000,F1,'histc')   % 
xlabel('Time (sec)');  ylabel('Firing rate (Hz)')
title(strcat('Spike train 1, bin size T=',num2str(T),' ms'))

figure(5)
bar(edges2/1000,F2,'histc')
xlabel('Time (sec)'); ylabel('Firing rate (Hz)')
title(strcat('Spike train 2, bin size T=',num2str(T),' ms'))

% --------------------
% let's have a look at the interspike intervals (ISIs)
 
for i =2:length(s1),  % generate a vector of ISIs
  isi1(i) = s1(i) - s1(i-1);
end
for i =2:length(s2),
  isi2(i) = s2(i) - s2(i-1);
end

bins = 100; % the number of bins (try different numbers)

figure(6)
hist(isi1,bins)
xlabel('Time (ms)')
ylabel('Number of spikes')
title('Interspike interval histogram of spike train 1')

figure(7)
hist(isi2,bins)
xlabel('Time (ms)')
ylabel('Number of spikes')
title('Interspike interval histogram of spike train 2')

% save the last two plots
print('-dpsc2','-f6','figure6.ps');
print('-dpsc2','-f7','figure7.ps');