Extraction of non-sinusiodal repetition rates

Question

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.

Answer 1

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
#sig = scipy.io.loadmat('raw1.mat')['raw1'][0]
sig = scipy.io.loadmat('raw2.mat')['raw2'][0]

# 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

Answer 2

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.

Answer 3

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.