# Matlab - How to create an autocorrelogram using a spike-train

I have already posted this question on MATLAB's question forum. However, I haven't received any responses, and am hoping to find some folks who may be a bit more signal processing savvy on here. Here goes... I'm wondering if someone could help clear up some confusion I am having surrounding the xcorr() function since I am struggling to parse how it works exactly.

r = xcorr(x) returns the autocorrelation sequence of x.

I am trying to make an autocorrelogram out of spike train data. Thus, I am inputting timestamps into xcorr() like so... xcorr(spiketimes). I am operating under the assumption that autocorrelating time-stamps for events (spikes) should lead to an autocorrelogram representative of the the autocorrelation between interspike intervals. This autocorrelogram, I believe, should provide the same info as an autocorrelogram of 0's and 1's (0 = no spike, 1 = spike) inputted into a vector who's length is the size of the recording (in milliseconds). When I run that line of code, however, I find that there's a peak towards the length of the timestamps vector. This doesn't make sense. The peak should be at lag zero (where the vector is 100% correlated with itself). What am I doing wrong here? Or is my understanding of an autocorrelogram incorrect? Does anyone here have any experience using xcorr(), and have any advice?

Also, this is less important and a bit technical, but for fellow matlab users, I'm not totally sure how to adjust bin width. If anyone has some pointers, I'd appreciate it. Finally, the fact that I can't center the plot that's created from xcorr() at zero is driving me crazy. How can I fix this?

• why is your input spiketimes and not the sequence of spikes itself? Oct 20 '15 at 20:44
• @panthyon Is my error as simple as that? If so, awesome. I just thought that if it was running an autocorrelation, then whether I used time-stamps or actual spikes did not matter. Autocorrelating time-stamps for events (spikes) should lead to an autocorrelogram representative of the the autocorrelation between interspike intervals. This would provide the same info as an autocorrelogram of 0's and 1's (0 = no spike, 1 = spike) inputted into a vector who's length is the size of the recording (in milliseconds), would it not? Oct 20 '15 at 21:20
• @ChrisEndemann Maybe you could edit the question and include in your assumption that Autocorrelating time-stamps for events (spikes) should lead to an autocorrelogram representative of the the autocorrelation between interspike intervals. This would provide the same info as an autocorrelogram of 0's and 1's (0 = no spike, 1 = spike) ? Oct 20 '15 at 21:40
• @ChrisEndemann "is my error as simple as that?" i don't know, why don't you try it? (: also the number of 'bins' are related to the length of the sequence (in fact $2N-1$ bins). if you modify the x-axis to indicate the bins in seconds, you might need to code a grouping mechanism that takes into account the sampling rate. i would not use the ISI as the input but indicate the lag on the x-axis of your autocorrelation. Oct 21 '15 at 0:20
• can you post plots and maybe some snippets of representative data? have you tried substituting timestamps for 1's? Oct 21 '15 at 15:38

When I run that line of code, however, I find that there's a peak towards the length of the timestamps vector.

So I'm assuming you're running:

r = xcorr(x);


and you're perplexed by the peak at r(length(x)) ? That's because xcorr returns a vector that is the autocorrelation from a lag of -length(x) to a lag of +length(x) (unless told otherwise) so that the middle value is the zero lag.

However you seem interested in finding correlations between the pulse times. That's a whole other kettle of fish... especially when there are pulse trains of different periods or there are missing data. To see the complexity that can happen, check out this paper.

What are you trying to achieve?

To see that there's a big difference between the autocorrelation function of the pulse arrival times and the actual pulses themselves, I've generated an example in R. The resulting autocorrelation functions are plotted below. --

Ok, so you have a periodic(ish) pulse train and you want to figure out what the pulse period is.

• The autocorrelation of the pulse times will not help. This is because the times are just a sorted array, so you'll just get a straight line as the autocorrelation.

• The autocorrelation of the pulses themselves should help. However, you'll get a peak at 0 lag, and at lags of $nP$ where $P$ is your pulse period and $n$ is an integer.

• The post you reference is looking at the histogram of the inter-pulse times (time differences between successive pulses). Depending on how varied your pulse train is, this might be the best approach.

If I change my simulation to use a pulse that is periodic (with zero noise on the pulse times!), then I get the second lot of code below, and the results are graphed here: For this particular run, the Tperiod = 9.350267. As you can see, both the ACF of the pulses and the histogram method get a peak close to where we want it (though quantization effects mean we don't get exactly at 9.350267; better selection of "bins" might improve that).

R Code Below

#26547

Npulses <- 100
Tmax <- 1000

pulseTimes <- sort(runif(Npulses,0,Tmax))

# install.packages("Hmisc")
library(Hmisc)

xc_times = acf(pulseTimes, plot = FALSE)

times <- 0:.1:Tmax

pulses <- times*0

for (pulseNo in 1:length(pulseTimes))
{
idx <- which.min(abs(pulseTimes[pulseNo] - times))
pulses[idx] = 1
}

xc_pulses = acf(pulses, plot = FALSE)

plot(xc_times$lag,xc_times$acf, type="l", ylim=c(-0.1,1))
lines(xc_pulses$lag, xc_pulses$acf, col="red")
legend(15, 0.9, c("ACF of times", "ACF of pulses"), lwd=c(2.5,2.5),col=c("black","red"))


Second version of R code below

#26547

Npulses <- 100
Tmax <- 1000

Tperiod <- runif(1,9,10)
Tstart <- runif(1,0,1)

#pulseTimes <- sort(runif(Npulses,0,Tmax))
pulseTimes <- Tstart + (0:Npulses)*Tperiod

# install.packages("Hmisc")
library(Hmisc)

xc_times = acf(pulseTimes, plot = FALSE)

times <- 0:.1:Tmax

pulses <- times*0

for (pulseNo in 1:length(pulseTimes))
{
idx <- which.min(abs(pulseTimes[pulseNo] - times))
pulses[idx] = 1
}

xc_pulses <- acf(pulses, plot = FALSE)

histogram <-  hist(diff(pulseTimes), breaks=xc_pulses$lag, plot=FALSE) plot(xc_times$lag,xc_times$acf, type="l", ylim=c(-0.1,1)) lines(xc_pulses$lag, xc_pulses$acf, col="red") lines(histogram$mids, histogram$counts / max(histogram$counts), col="green")

legend(12, 1.0, c("ACF of times", "ACF of pulses", "Histogram of inter-pulse times"), lwd=c(2.5,2.5, 2.5),col=c("black","red", "green"))

• Thanks so much for clearing that up. Pardon the confusion, but I do not think that I am specifically interested in the the correlations between pulse times. I simply want to know whether or not the neurons I have data on show predictability in their firing patterns. I think that I could do this by running an autocorrelation of the spike-train vector, or by running an autocorrelation between the time-stamps of events (spikes). Honestly, I don't see the functional difference between the two given the output that I wish to obtain. Care to offer me some insight? Oct 21 '15 at 21:49
• Thank you very much for the visual comparison between the two functions. Could you elaborate a bit more, in words, on how they differ? I only ask because the following post, stats.stackexchange.com/questions/148865/…, uses an ACF of times, but still gets a graph that looks very similar to an ACF of pulses. Oct 21 '15 at 22:17
• That diagram is a histogram of time differences, which is different again. I'll try to reply in more detail, but it'll be an hour or two. Have to put the babies to bed. :-)
– Peter K.
Oct 21 '15 at 22:20
• @ Peter K. You sir, have been a great deal of help. Thank you for your time and effort. I have often heard people refer to an autocorrelation of spike times as an autocorrelation histogram (web.mit.edu/hst.723/www/Labs/LabANF.htm). I suppose the histogram of time differences may be the most common option. Anyways, I have written some code to produce an "autocorrelation histogram" with help from this post: neuro.bio.lmu.de/teaching/statistical-models_ss2013/solutions/… . If you find the time, could you verify that my code is correct? Posted in answers. Oct 24 '15 at 22:36

Chris, since you say you are not from signal processing the definition of a spike for this link might not make sense, but the spike train by definition is a series of impulses time shifted by each spike time so its like $[1 \, 0\, 0\, 1\, 0\, 1\, 0\, 0]$ and not $[t_{1}\, 0 \, 0 \, t_{2} \, 0 t_{3} \, 0]$ etc. see http://icwww.epfl.ch/~gerstner/SPNM/node34.html#SECTION02424200000000000000

You say you are interested in finding the degree of randomness of your spike train, a time-series process. Autocorrelation is part of that but will not alone tell you with statistical exactness the odds that a give correlation would be spuriously generated (i.e. greater than an expected value of the ACF).

However, as a quick check, you can look at the first lag of your spike train and see how close to 1 it is in the normalized autocorrelation. Values closer to 1 will tell you there is a stronger (or weaker) predictability to the time-series process used as your input. At lag zero in normalized autocorrelation the value will always be 1.

function [ac,x] = Auto_corr(ts, binsize, n_lags)
% INPUT: ts = timestamps (sorted)
%        binsize = binsize for binning timestamps.
%        n_lags = n_lags in the xcorr.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Cowen 2021
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
plot_it = false;
scaleopt = 'none';

ac = nan(n_lags,1);
x = linspace(binsize/2,(n_lags*binsize) - binsize/2,n_lags)' ; % bin centers.

if isempty(ts)
return
end

edges = ts(1):binsize:ts(end);
c = histcounts(ts,edges);
[ac] = xcorr(c,n_lags,scaleopt);
ac = ac((n_lags+2):end);