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

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 '12 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 '12 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 '12 at 0:26

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.