I have $n$ discrete samples of a seismic signal $y[n]$: enter image description here

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!

  • 1
    $\begingroup$ Maybe look at question 1 and question 2. $\endgroup$
    – Geerten
    Mar 21, 2012 at 11:00
  • $\begingroup$ Several methods for frequency spectra: dspguru.com/dsp/howtos/how-to-interpolate-fft-peak $\endgroup$
    – endolith
    Mar 26, 2012 at 20:51
  • $\begingroup$ That second one doesn't have an answer @Geerten ;-) $\endgroup$
    – Ivo Flipse
    Mar 28, 2012 at 9:32
  • 1
    $\begingroup$ Oh..haha, good point. Well I'll reference back to this question on that question ;) $\endgroup$
    – Geerten
    Mar 28, 2012 at 9:37

3 Answers 3


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.

  • $\begingroup$ The combination of upsampling plus parabolic interpolation may work better than either alone. $\endgroup$
    – hotpaw2
    Mar 22, 2012 at 19:32

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.

  • $\begingroup$ You have both a x[n] and y[n] in your equation. Is this correct or should it be just y[n]? $\endgroup$
    – Andy
    Mar 21, 2012 at 17:11
  • $\begingroup$ x[n] is the input, y[n] is the output. Bad answer overall, there is a typo (fixed now) and I misunderstood the question. Apologies $\endgroup$
    – Hilmar
    Mar 22, 2012 at 0:26

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;

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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.