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
- finds the derivative of the signal by convolving the signal with the 4th order SGF coefficients for the derivative
- finds pair of samples (i,i+1) where the derivative changes sign
- 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!