Finding local peaks in-between samples

Question

I have $n$ discrete samples of a seismic signal $y[n]$:

I want to find local maxima in the signal.

A naive test for if $y[n]$ is a maximum would be: $$y[n]: maxima \textbf{ if } y[n] > y[n-1] \textbf{ and } y[n] > y[n+1]$$

However the maxima are probably located in between samples, e.g. there may be a maximum at $i=4.25$.

In order to find maxima in between samples I believe that I need to interpolate $y[n]$.

• How do I find maxima using interpolation ?
• What form of interpolation should I use ?

As you can see my signal is not very noisy, however it would be good if the method also did a bit of filtering so that the maxima exceed a treshold and have a certain width (no spikes).

My biggest problem however is just to find peaks in between samples. Any suggestions for a good way to do this ?

Thanks in advance for any answers!

Answer 1

Getting a sub-sample resolution

A very cheap (in terms of code size) solution is just to upsample your signal. In matlab, this can be done with interp(y ,ratio). A slightly more complicated solution consists in naively detecting peaks ; and for each peak, fitting a parabola through y[peak - 1], y[peak], y[peak + 1] ; then using the point at which this parabola is maximal as the true peak position.

Regarding peak detection

A bunch of techniques which help:

• As suggested by Hilmar, convolving the signal by a Gaussian or Hann window, the width of which is roughly equal to half the minimum interval you want to see between detected peaks. Since temporal accuracy seems essential to your application, make sure that you take into account the time delay introduced by the filtering, though!
• Subtract to your signal a median filtered version of itself (with a fairly large observation window) ; and divide the result by a standard-deviation filtered version of itself. This gets rid of trends and allows the thresholds to be expressed in units of standard deviations.
• For peak-picking, I formulate that using a "top-hat" filter. Define the top-hat filtered version of your signal as yt[n] = max(y[n - W], y[n - W + 1], ..., y[n + W - 1], y[n + W]) ; and use as peaks the points where y[n] == yt[n] and y[n] > threshold.

All this can be very efficiently implemented in Matlab with a few passes of nlfilter.

Answer 2

Try a lossy peak detector:

y[n] = max(abs(x[n]),a*y[n-1]);

where "a" is a number smaller than 1 that controls how fast the detector decays. It determines how close neighboring peaks can be without smooshing into a single one. Then do a threshold detection.

Answer 3

You can also use a mean shift algorithm to refine your peaks using a kernel function.

This can especially useful and efficient if you have very few peaks that you care about compared to the size of the signal, so you won't have to upsample the entire signal (which can be very expensive).

For example here is a simplified MATLAB or GNU Octave implementation:

function peak = meanshift(y, sigma, peak_guess)
    peak = peak_guess;
    x = 1:size(y, 2);

    for iter = 1:10
        %TODO: ignore areas more than 3 sigma away
        %      since those don't affect the output
        weights = exp(-(x - peak).^2 / sigma^2) .* y;
        off = sum(weights .* x) / sum(weights);
        peak = off;
    end
endfunction

although I must emphasize that you have to implement the thing described in the "TODO" in order for it to be efficient.