smoothing - Finding local peaks in-between samples

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

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.