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"))
