% Example 1: Generate a pair of Poisson spike trains with
%            a parabolic cross-correlation
% To call this function from Matlab, just type Example1
% into the command line and press enter.
% This program might take a while to run because
% computing cross-correlation functions takes a 
% while.


% First define some parameters:
rate=.3;    % Rate of the processes
T=250000;   % Total length of each process
dt=.1;      % bin size
winrad=20;  % window radius over which to compute and 
            % plot the cross-correlation

% First define a probability mass function (w)
% over a domain (x)
% We normalize w to make sure that it sums to 1

x=0:.1:10;
w=-(x-5).^2+25;
w=w/sum(w);

% Plot this function so that we can compare it to 
% the cross-correlation function
subplot(2,1,1);
plot(x,rate*w/dt);
xlabel('x');
ylabel('w');
axis([-winrad winrad 0 max(w*rate/dt)]);

% Now create the first Poisson spike train
% It is a million time-units long with a bin
% size of .1 time units
% It has a rate of .2 spikes per time unit
p1=Poisson(rate,T,dt);

% Now create the second spike train by jittering
% the first with random numbers pulled from the
% distribution w
% Pay attention to the syntax for passing the function
% handle to jitter.  It's a bit tricky.
p2=jitter(p1,@()randnum(w,x),dt);

% Calculate the cross-correlation function
% cc12 is the cc function and tau is its domain
[cc12,tau]=cross_corr(p1,p2,winrad,dt);

% Now plot the cross-correlation function
% Bar graphs work well for this data (since 
% it is basically a normalized histogram)
subplot(2,1,2);
bar(tau,cc12);
xlabel('time');
ylabel('cross-correlation');
axis([-winrad winrad 0 max(cc12)]);