Extraction of non-sinusiodal repetition rates

I have an auto-correlation function that was generated from a signal, and I am trying to extract its 'repetition rate' in order to calculate the dominant frequency of the pulse, but I am not exactly clear how to do this.

Here are two cases, labelled 'good' and 'bad' to mark best/worst case scenarios. What methodologies exist that I can use here? Since the repetitions are not sinusoidal, the spectrum does not seem to yield good information.

Here is the original 'good' case:

followed by its auto-correlation.

Similarly, here is the original 'bad' case:

again followed by its auto-correlation:

EDIT: Raw data files as requested:

Raw1

Raw2

EDIT On Feedback:

Spectral analysis was suggested, however I did not have too much luck with it as the spectrum does not look too sharp around the region where you would expect it to. I think this may to due to the fact that the signal is not a repetitive sinusoid and so projecting it on sinusoidal bases does not do too well.

• Could you share the data used to create these plots?
– Emre
Apr 19 '12 at 23:06
• Isn't the period just the position of the first spike? 650 samples or so? Can you show the original signal? I'm not sure why you're saying one of these is "good" and the other "bad". Apr 20 '12 at 0:33
• @Emre Ok, I have added the original signals. Apr 20 '12 at 3:28
• I would be surprised to find that spectral analysis doesn't provide any useful results for the above datasets, but you might also try calculating the autocorrelation of the autocorrelation signal itself, then searching for a peak there to find the lag that maximizes the autocorrelation function's self-similarity. Apr 20 '12 at 3:29
• @endolith Please see, I have added the original signals. You are right that the repetition is around 650 by how exactly do I compute that automatically? Seems like a peak-picking problem to me? Or is there some other methods that can be used? Apr 20 '12 at 3:29

You are right that the repetition is around 650 by how exactly do I compute that automatically? Seems like a peak-picking problem to me? Or is there some other methods that can be used?

Yes, it's just peak-picking. Your period is the x value of the first strong peak:

Your peaks are all similar in height, probably because you're doing the autocorrelation using the FFT? So it's a circular autocorrelation. You can

1. make it non-circular by zero-padding before doing the FFT, or
2. could skew the plot you currently have by adding a linear function that emphasizes the peaks close to 0 and reduce the height of the farther ones.

It depends on what you're specifically looking for. In either case you then have to pick the highest peak that's not at 0 lag.

I adapted my script from here and it (surprisingly) worked without any tweaking:

raw1 x = 709.23

raw2 x = 710.77

from pylab import *

# Import matlab files
import scipy.io

# Calculate autocorrelation (same thing as convolution, but with
# one input reversed in time)
corr = fftconvolve(sig, sig[::-1], mode='full')

# throw away the negative lags, so x = 0 means 0 lag
corr = corr[len(corr)/2:]

# Find the first low point from derivative
d = diff(corr)
start = find(d > 0)[0]

# Starting at the first low point, find the highest peak
peak = argmax(corr[start:]) + start

# Fit parabola to estimate more precise location of peak
px, py = parabolic(corr, peak)

plot(corr)
axvline(px, color='r')

print px


The variable peak is the integer x value at the peak (709). Then the parabolic function fits a parabola to the peak and tries to get a more precise estimate (709.232...), which may or may not be relevant to what you're doing.

The python fftconvolve() function does zero-padding internally to produce a linear convolution, like convolve() does

• I just did an xcorr honestly. My question to you is how are you exactly picking this (next) biggest peak? (For an automation to do it, it seems to be harder than it looks)... Apr 20 '12 at 4:06
• @Mohammad: Yes, xcorr is circular. When you do zero-padding first, the autocorrelation becomes sloped like this, with closer peaks higher than farther ones. Then I find the first local minima, then search for the highest peak from there. It's not foolproof by any means, but it works for these examples. Apr 20 '12 at 4:22
• Ah! Thats very nice and creative actually... I see what you did. (And instead of diff you can maybe filter with [-1 0 0 0 1] to make it more resilient to fake minima as you descend...). I like it. :-) I think I get everything - then what is this last 'parabolic' command? Last thing, do you think that all the other repetitive peaks give us any extra information? Apr 20 '12 at 4:32
• @Mohammad: Actually, this says that xcorr is not circular. I'm not sure. Sorry, parabolic function is here gist.github.com/255291#file_parabolic.py It just fits a parabola to the peak to try to estimate the true peak more accurately between samples. Apr 20 '12 at 4:36
• I should add, the reason my autocorrelation is not skewed is because I was doing an 'unbiased' xcorr instead of the default 'biased' one. xcorr(signal,'unbiased'); Apr 20 '12 at 4:40

One would expect such a sequence to have a spectrum consisting of lines, as it is almost periodic (if it was periodic, it would have a Fourier series representation, even though it is not sinusoidal). As a quick example:

load raw1.mat
% calculate "unbiased" normalized cross-correlation; adjusted for
% regions where there isn't full overlap
corr = xcorr(raw1,'unbiased');
% calculate oversampled FFT to help in illustration
C = fft(corr, 65536);
% plot FFT to see peaks
figure; plot((-32768:32767)/65536, abs(fftshift(F))); grid on;
% zoom in to interesting area (X-axis is in normalized frequency)
set(gca,'XLim',[0 0.01]);


The resulting plot:

The first peak is at zero frequency, due to the nonzero mean of the autocorrelation signal. The next, fundamental frequency peak is at approximately $0.00141$, which compares favorably to the previously-determined estimate of $\frac{1}{650} = 0.00153$. This implies that the overall period of the autocorrelation sequence is approximately $\frac{1}{0.00141} = 709$ samples.

• Thanks! I demeaned it and got the large peak at 0.00141 as well. Some follow ups? 1) Why are you oversampling the FFT here? I understand the visual reason but whats the DSP reason? 2) Do the other peaks give any additional information thats relevant? Thanks again Apr 20 '12 at 4:07
• 1: It's for visual effect only. You can absolutely use a smaller transform length and interpolate as endolith showed in his answer. 2: In this case, not really. The other peaks will be located at multiples of the fundamental frequency. If you're just concerned with estimating that frequency, then in a noisy environment you might be able to use all the peak locations to get a better estimate of the fundamental spacing, but I doubt this would be too valuable. Apr 20 '12 at 13:04

If you can model your noise and signal then a good technique might well be bayesian inference for a general linear model. This book has an excellent introduction to the techniques and demonstrates them applied to inferring the frequency of a single sinusoid in (lots of) noise. Although this is not really my area, I think it will apply to a more complex model such as yours (as long as you can formulate it as a GLM). It might even work on the single sinusoid model. The link is to Amazon's page where they kindly provide a partial preview up to the chapter of interest (chapter 2). Another link by the same author which is more complete compared to the amazon page is on a change point detector using a GLM

I should prefix the following by reiterating that this is not my area, so I hope there are no errors, but I don't guarantee it. For a proper treatment, see the references.

You can represent your signal as a general linear model if you can represent it as follows: $$\mathbf{d} = \text{G}\mathbf{b} + \mathbf{e}$$ where $\mathbf{d}$ is the Nx1 vector of measured noisy data points, $\mathbf{e}$ is the Nx1 vector of noise samples. $\text{G}$ is the NxQ matrix of basis vectors describing the signal. This is parameterised as you see fit (ideally with as few parameters as possible). Think of this as the function that can you can use to create your measured signal in the absence of noise. You get one $\text{G}$ for each parameter you choose, and you need to compute the probability at that point, so too many parameters makes the problem intractable. $\mathbf{b}$ is the Qx1 weighting vector of the basis functions given by the columns of $\text{G}$. Imagine you signal is some sinusoid, then the basis functions are a cosine and a sine wave (of the same frequency). If it was a pure cosine, then the weighting for the sine basis would be zero. If it was a pure sine wave, then the weighting for the cosine basis would be zero. Anything else and you'd need to sum a bit of each.

Under some basic assumptions about the noise being Gaussian and the priors on the other parameters (discussed in the reference book), one can come up with a marginal posterior distribution on the basis parameters:

$$p\left(\mathbf{w}_m|\mathbf{d},M\right) \propto \frac{\left[\mathbf{d}^{\text{T}}\mathbf{d} - \mathbf{d}^{\text{T}}\text{G}\left(\text{G}^{\text{T}}\text{G}\right)^{-1}\text{G}^{\text{T}}\mathbf{d}\right]^{\frac{-\left(N-Q\right)}{2}} }{\sqrt{\text{det}\left(\text{G}^{\mathbf{T}}\text{G}\right)}}$$

I'm a bit flakey as to what $M$ represents, but I think it's a placeholder for "what we know and assumed". $\mathbf{w}_m$ denotes the parameters of the basis functions (which essentially define $\text{G}$). By calculating the above for lots of $\mathbf{w}_m$ vectors, you get something akin to a probability distribution over the space of little $m$ (albeit, unnormalised).

Essentially what you yield with these techniques is a probability distribution over the parameter that you care about (in this case, frequency), from which you can pick the "most likely" value. The point is that you factor knowledge about what to expect into the inference, this means you're massively restricting your search space.

• Thanks Henry. Can you please expand on what you mean exactly by "(as long as you can formulate it as a GLM)". When you say, 'model it as a generalized linear model', what does that mean exactly? Model the spikes? Linear in what sense? Apr 20 '12 at 17:19
• Right, I've found another paper that describes it more fully (it's the same author, just without the bits I noticed were missing from the Amazon preview!). The paper describes a change point detection algorithm. I'll try to add a bit more information to the answer. Apr 20 '12 at 18:00
• Ah! Thanks a lot for that. Seems like a great read... I will have to study it further. I have traditionally shied away from Bayesian methods because I never understood how one goes, explicitly, from data points in a vector to a probability distributions described. Perhaps this paper will help in that regard. Apr 20 '12 at 18:10
• I given my synopsis of the techniques. I'm open to being corrected on any aspect of it. Apr 20 '12 at 18:38
• Is little 'm' here the index of the family of basis functions G you are cycling over to create the probability density function? Apr 20 '12 at 19:54