How do I use a Savitzky Golay filter to find local maxima (in between samples) in a discretely sampled 1D signal?

Question

I have a seismic signal y(i):

Here I have found one maximum: i=152.54, y=222.29 manually and plotted it in red.

I want to find all maxima automatically.

I read that the Savitzky Golay Filter (SGF) can be used to find smoothed estimates of both a signal and its derivatives, and that one of the benefits of the SGF is that it preserves minima and maxima much better than other filters. This sounds great for my use.

I found a Matlab script that generates SGF coefficients. And used this to find that the 4th order SGF coefficients for the derivative. I coded a small Matlab script that

1. finds the derivative of the signal by convolving the signal with the 4th order SGF coefficients for the derivative
2. finds pair of samples (i,i+1) where the derivative changes sign
3. finds zero crossing of derivative by linear interpolation between i and i+1

Script:

function [maxX,maxY] = findLocalMax(y)
% Kernel for 4th order Savitzky-Golay filter for finding derivative:
d4 = [0.0724 -0.1195 -0.1625 -0.1061 0 0.1061 0.1625 0.1195 -0.0724];

dy = conv(y,d4,'same'); % derivative

[m n] = size(dy);
maxX = [];
maxY = [];
for i = 1 : n - 1
  if dy(i) < 0 && dy(i+1) > 0 % max somewhere between i and i+1
    a = dy(i)/(dy(i) - dy(i+1)); % linear interpolation
    mx = i + a;
    maxX = [maxX mx];
    my = y(i)*(1-a) + y(i+1)*a; % linear interpolation
    maxY = [maxY my];
  end
end

Im my script I had to test for if the derivative changes from negative to positive to get the function to give the desired result, however this confuses me. Should not the derivative for a maximum go from positive to negative? Is there some better way to distinguish between maxima and minima?

Below is the result of using this function to find the maxima on my signal:

The results look good, but I notice that some maxima are not found: i= 143.13, 190.88, 256.97.

Is this because they are to close to other maxima?

How can I control the nearest two maxima may be?

Thanks in advance for any answers!

Answer 1

While I'm not familiar with this specific type of filter, based on the plot you've shown, I would guess that the maxima that aren't found by your process are just butting up against the time resolution inherent in the process. Any kind of "smoothing" implies that there is some time-local smearing of the signal of interest, such that if there are two nearby peaks, it's possible that they will be coalesced into one. It's possible that a lower-order filter might present less of this behavior, likely at the expense of the amount of smoothing that you get.