Finding periodicity of spikes, matlab

Question

In my time series data, I have often spikes with regular periods. Sometimes, I get two or three different periodical spike sequences (two electrical noise sources), some with T1 seconds period, some with T2 second period and some irregular and at random times.

Lets say I have two groups of overlapping periodical spikes in my data, one with a period of 8 seconds and another with 6. I have zeroed the data except for the times where there is spikes. So, my data now looks like a train of spikes. How can I find the periodicity of the two and separate the two groups? Cross-correlation didn't help as many are just single sample spikes. I was thinking of FFT. But I can't see anything in the spectra at those low frequencies, nor when I use FFT on original data. My data has 500 Hz. sampling frequency. I do not want to remove them, just to find the periodicity. Thanks.

Answer 1

My Scipy code that generates data, plots it and performs the fft is listed below. The code generates two "spike trains" at intervals of 6 and 8 seconds. The following plot shows the signal in time, and its fft (sorry, the fft title got mangled).

The interesting thing to notice is that the frequencies you are looking for are there. Said frequencies are

1/8 Hz or omega = 1/8/500 = 0.00025 omega and

1/6 Hz or omega = 1/6/500 = 0.00033

Which is where the first two peaks are in the frequency domain. Further out are many more repeating peaks/aliases because these are essentially really extreme, intermittent square waves. You could try low-pass-filtering to isolate the fundamental frequencies and eliminate other peaks, but I think that's probably a lost cause at a sample rate of 500 Hz. Trying to get a 0.2 Hz cutoff filter @ 500 Hz seems a bit much. Just zooming in on the frequency region of interest might be good enough. The next step could be to understand exactly how these spike trains are going to cause repetition/aliasing in the frequency domain and look for a certain pattern.

Or try to either LPF and downsample a lot and retry the FFT route, calculate individual DFT bins using straight-up DFT math or use the Goertzel algorithm to calculate DFT bins for suspected frequency locations. Maybe try a Goertzel algorithm from 0.0 Hz to 0.2 Hz with 0.01 spacing - just as an example.
Good luck - it doesn't seem like a trivial problem.

from pylab import *
close('all') #close previous plots

fs=500 # generate signal
s1 = hstack(( 1,zeros(fs*6-1) )) # 1/6 Hz
s2 = hstack(( 1,zeros(fs*8-1) )) # 1/8 Hz
s1r = tile(s1,(12,))
s2r = tile(s2,(9,))
s = s1r + s2r
print "shape of s is %s"%(str(s.shape))

sfft=fft(s) # execute fft
sfreqs=fftshift(fftfreq(len(sfft)))
print "shape of sfft is %s"%(str(sfft.shape))

figure(); # generate plots
subplot(211); plot( 1./fs * arange(len(s)) , s ) #, 'o' )
title('signal, s');xlabel('sec');ylabel('mag')
subplot(212); plot( sfreqs , (abs(sfft)) ) #, 'o' )
title('fft(s)');xlabel('$\omega$');ylabel('linear mag')
axis([0.0,0.0005,0.0,25.])
show()

UPDATE: Thinking about it some more, the best solution probably doesn't involve signal processing. You may want to threshold the samples in to zeros and spikes. Then take the interval between the spike locations 0 and 1 and look for that interval through the rest of the spikes. If that interval is consistent with all of your spikes, then you've found that series of spikes and should "subtract out" that series from the spikes. Then keep iteratively looking at intervals of spikes. In a noise, imperfect system, a good implementation of this will probably be tricky, but it's probably a better solution to this specific problem.