% ONE3   Exercise 3 from CHAPTER 1 of Dayan and Abbott

%   NEUR 1680
%   Last updated 02/08/2012

% Compute autocorrelation histograms of spike trains generated by:
% A. a Poisson generator with a constant fring rate of 100 Hz, 
% B. a constant firing rate of 100 Hz together with a refractory period modeled as in exercise 2 with tau_ref = 10 ms,
% C. a variable firing rate rate(t) = 100(1+cos(2*pi*t/25)) (t is in ms and rate(t) is in Hz).

% Plot the histograms over a range from 0 to 100 ms.



% A. CONSTANT FIRING RATE

clear
close all
set(0,'DefaultFigureWindowStyle','docked')

%%%%%%  PARAMETERS   

% try different values for duration and lags...
% try zero_mean_flag = 0 or 1

duration = 100;                     % duration of recording [sec]
dt = .001;                         % discrete time step    [sec]
lags = -.1/dt:.1/dt;                % lags, in bin units
zero_mean_flag = 1;                 % set to 1 to get an auto-covariogram
r_0 = 100;                          % firing rate (Hz)

%%%%%%%%%%%%%%%%%%%

N = round(duration/dt);             % number of time steps in recording

p = r_0*dt;                         % probability of observing a spike in any bin
p = min(p,1);

% The Poisson spike train is generated by independent Bernoulli random variables:
spike_train = rand(1,N) < p;        % independent Bernoulli trials,  


autoco = autocorrelation(spike_train, dt, lags, zero_mean_flag);
figure(1)
plot(lags, autoco)
grid
title('AUTOCORRELOGRAM: Poisson spike train')
xlabel(['bin = ' num2str(1000*dt) ' ms'])
ylabel('spikes/sec  *  spikes/bin')


%% B. CONSTANT FIRING RATE WITH REFRACTORY PERIOD

clear

%%%%%%  PARAMETERS   

% try different values for duration and lags...
% try zero_mean_flag = 0 or 1

duration = 100;                     % duration of recording [sec]
dt = .001;                          % discrete time step    [sec]
lags = -.1/dt:.1/dt;                % lags, in bin units
zero_mean_flag = 1;                 % set to 1 to get an auto-covariogram
r_0 = 100;                          % firing rate (Hz)
tau_refractory = .010;              % time constant of recovery of firing rate after a spike (sec)

%%%%%%%%%%%%%%%%%%%

N = round(duration/dt);             % number of time steps in recording

rate = r_0;

for n = 1:N
    p = rate*dt;                    % probability of observing a spike in bin n
    spike = rand < p;
    spike_train(n) = spike;
    rate = (1 - spike) * (rate + (dt/tau_refractory)*(r_0 - rate)); % depending on spike, rate resets or recovers 
end

autoco = autocorrelation(spike_train, dt, lags, zero_mean_flag);
figure(2)
plot(lags, autoco)
grid
title('AUTOCORRELOGRAM: Poisson with refractory period')
xlabel(['bin = ' num2str(1000*dt) ' ms'])
ylabel('spikes/sec  *  spikes/bin')



%% C. SINUSOIDAL FIRING RATE

clear

%%%%%%  PARAMETERS   

% try different values for duration and lags...
% try zero_mean_flag = 0 or 1

duration = 100;                     % duration of recording [sec]
dt = .001;                          % discrete time step    [sec]
lags = -.1/dt:.1/dt;                % lags, in bin units
zero_mean_flag = 1;                 % set to 1 to get an auto-covariogram
r_0 = 100;                          % firing rate (Hz)
tau_refractory = .010;              % time constant of recovery of firing rate after a spike (sec)

%%%%%%%%%%%%%%%%%%%

N = round(duration/dt);             % number of time steps in recording


rate = r_0*(1 + cos(2*pi*(dt:dt:duration)/.025));  % sinusoidal firing rate (Hz)

p = rate*dt;                        % probability of observing a spike in a given bin
p = min(p,1);


% The Poisson spike train is generated by independent Bernoulli random variables:

spike_train = rand(1,N) < p;    % independent Bernoulli trials,  


autoco = autocorrelation(spike_train, dt, lags, zero_mean_flag);
figure(3)
plot(lags, autoco)
grid
title('AUTOCORRELOGRAM: sinusoidal firing rate')
xlabel(['bin = ' num2str(1000*dt) ' ms'])
ylabel('spikes/sec  *  spikes/bin')

