Detecting periodicity in point processes

Question

I have data from a periodic one dimensional point process representing the times of events which is also mixed in with lots of data that I regard as noise. The total number of points is of order one million. I would like to detect (approximate) periodicity in the data to reconstruct the periodic point precess. This question is a little vague but what is a good method for doing this?

As an example, if my points are $1,10,25,30,31,50$ I would like to detect $10, 30, 50$ as equally spaced points. As my values are reals I need this to tolerate small errors in the data.

Answer 1

What you are looking for are pulse trains in noise. The pulse train is its own Fourier transform (sort of) and noise (usually) is spread out evenly in the frequency domain.

1. Make vector of all zeros. Call it x.
2. For each time value, set that index of x to 1. For instance, if your points are 1,10,25,30,31,50, then set x(1)=1, x(10)=1, x(25)=1, and so on.
3. Take the DFT of this sequence and plot it's amplitude.
4. There will be a series of regularly-spaced spikes corresponding to the periodic components of the input. Their spacing will be the reciprocal of the period in the time domain.

Here's is a quick and dirty example I put together in Matlab:

N = 1e6;
x = zeros(N,1);
x(1:33:end) = 1;             % Pulse train 1.
x(1:137:end) = 1;            % Pulse train 2.
randinds = randperm(N);      % Random permutation for noise.
x(randinds(1:10e4)) = 1;     % Add some noise points.
X = fftshift( fft( x ) );
f = (-0.5:1/N:0.5-1/N);
plot( f, abs( X ) );

There are two series of pulse trains visible in the plot of the magnitude of the FT of x. By visual inspection, one is spaced at 0.0303 and the other at 7.299e-3, which gives time domain periods of 33.0033 and 137.005, respectively.

You may have to tailor the above to your exact situation, but I think something along these lines is doable and easy to implement.