Finding an approximate local maximas with noisy data in Matlab

Question

Hi,

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

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 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(X*K + K log K). I have not looked at the source code and other implementations are possible, but basically it should be O(X*K).

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

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.