Finding an approximate local maximas with noisy data in Matlab

Question

The matlab FAQ describes a one-line method for finding the local maximas:

index = find( diff( sign( diff([0; x(:); 0]) ) ) < 0 );

But I believe this only works if the data is more or less smooth. Suppose you have data that jumps up and down in small intervals but still has some approximate local maximas. How would you go about finding these points? You could divide the vector into n pieces and find the largest value not on the edge of each but there should be a more elegant and faster solution.

A one-line solution would be great, too.

Edit: I'm working with noisy biological images which I am attempting to divide into distinct sections.

Answer 1

I'm not sure what type of data you're dealing with, but here's a method I've used for processing speech data that could help you locate local maxima. It uses three functions from the Signal Processing Toolbox: HILBERT, BUTTER, and FILTFILT.

data = (...the waveform of noisy data...);
Fs = (...the sampling rate of the data...);
[b,a] = butter(5,20/(Fs/2),'low');  % Create a low-pass butterworth filter;
                                    %   adjust the values as needed.
smoothData = filtfilt(b,a,abs(hilbert(data)));  % Apply a hilbert transform
                                                %   and filter the data.

You would then perform your maxima finding on smoothData. The use of HILBERT first creates a positive envelope on the data, then FILTFILT uses the filter coefficients from BUTTER to low-pass filter the data envelope.

For an example of how this processing works, here are some images showing the results for a segment of recorded speech. The blue line is the original speech signal, the red line is the envelope (gotten using HILBERT), and the green line is the low-pass filtered result. The bottom figure is a zoomed in version of the first.

SOMETHING RANDOM TO TRY:

This was a random idea I had at first... you could try repeating the process by finding the maximas of the maximas:

index = find(diff(sign(diff([0; x(:); 0]))) < 0);
maxIndex = index(find(diff(sign(diff([0; x(index); 0]))) < 0));

However, depending on the signal-to-noise ratio, it would be unclear how many times this would need to be repeated to get the local maxima you are interested in. It's just a random non-filtering option to try. =)

MAXIMA FINDING:

Just in case you're curious, another one-line maxima-finding algorithm that I've seen (in addition to the one you listed) is:

index = find((x > [x(1) x(1:(end-1))]) & (x >= [x(2:end) x(end)]));

Answer 2

Depending on what you want, it's often helpful to filter the noisy data. Take a look at MEDFILT1, or using CONV along with FSPECIAL. In the latter approach, you'll probably want to use the 'same' argument to CONV and a 'gaussian' filter created by FSPECIAL.

After you've done the filtering, feed it through the maxima finder.

EDIT: runtime complexity

Let's say the input vector has length X and the filter kernel length is K.

The median filter can work by doing a running insertion sort, so it should be O(XK + K log K). I have not looked at the source code and other implementations are possible, but basically it should be O(XK).

When K is small, conv uses a straight-forward O(X*K) algorithm. When X and K are nearly the same, then it's faster to use a fast fourier transform. That implementation is O(X log X + K log K). Matlab is smart enough to automatically choose the right algorithm depending on the input sizes.

Answer 3

If your data jumps up and down a lot, then the function will have many local maxima. So I'm assuming that you don't want to find all the local maxima. But what is your criteria for what a local maximum is? If you have a criteria, then one can design a schem or algorithm for that.

I'd guess right now that maybe you should apply a low-pass filter to your data first, and then find the local maxima. Although the positions of the local maxima after filtering may not exactly correspond to those before.

Answer 4

There are two ways to view such a problem. One can look at this as primarily a smoothing problem, using a filtering tool to smooth the data, only afterwards to interpolate using some variety of interpolant, perhaps an interpolating spline. Finding a local maximum of an interpolating spline is an easy enough thing. (Note that you should generally use a true spline here, not a pchip interpolant. Pchip, the method employed when you specify a "cubic" interpolant in interp1, will not accurately locate a local minimizer that falls between two data points.)

The other approach to such a problem is one that I tend to prefer. Here one uses a least squares spline model to both smooth the data and to produce an approximant instead of an interpolant. Such a least squares spline has the advantage of allowing the user a great deal of control to introduce their knowledge of the problem into the model. For example, often the scientist or engineer has information, such as monotonicity, about the process under study. This can be built into a least squares spline model. Another, related option is to use a smoothing spline. They too can be built with regularizing constraints built into them. If you have the spline toolbox, then spap2 will be of some utility to fit a spline model. Then fnmin will find a minimizer. (A maximizer is easily obtained from a minimization code.)

Smoothing schemes that employ filtering methods are generally simplest when the data points are equally spaced. Unequal spacing might push for the least squares spline model. On the other hand, knot placement can be an issue in least squares splines. My point in all of this is to suggest that either approach has merit, and can be made to produce viable results.