Single dimension peak fitting

Question

I have a single dimensional array of floating point values (c# doubles FYI) and I need to find the "peak" of the values ... as if graphed.

I can't just take the highest value, as the peak is actually a plateau that has small fluctuations. This plateau is in the middle of a bunch of noise. I'm looking find a solution that would give me the center of this plateau.

An example array might look like this:

1,2,1,1,2,1,3,2,4,4,4,5,6,8,8,8,8,7,8,7,9,7,5,4,4,3,3,2,2,1,1,1,1,1,2,1,1,1,1

where the peak is somewhere in the bolded section.

Any ideas?

Answer 1

You can apply a low-pass filter to your input array, to smooth out the small fluctuations, then find the peak in the filtered data. The simplest example is probably a "boxcar" filter, where the output value is the sum of the input values within a certain distance from the current array position. In pseudocode, it would look something like this:

for i = 0, samplecount-1
  if (i < boxcar_radius) or (i >= (samplecount - boxcar_radius))  then
       filtered_data[i] = 0 // boxcar runs off edge of input array, don't use
  else
    filtered_data[i] = 0
    for j = i-boxcar_radius, i+boxcar_radius
       filtered_data[i] = filtered_data[i] + input_data[j]
    endfor
  endif
endfor

If you have some idea how wide the "plateau" will be, you can choose the boxcar radius (approximately half the expected plateau width) to detect features at the appropriate scale.

Answer 2

Peak detection is one of the stages in Phase Correlation and other motion estimation algorithms used in places like video compression. One approach is this: consider a candidate for a peak and a window of a certain number of neighbours. Now fit a quadratic function using standard regression. The peak, with subpixel accuracy, is at the maximum of the fitted quadratic.

Answer 3

Obviously exact solution depends on details. If your distribution is always nice as in your example you could have:

def GetPeak(l):
  large = max(l) * 0.8
  above_large = [i for i in xrange(len(l)) if l[i] > large]
  left_peak = min(above_large)
  right_peak = max(above_large)
  return (left_peak, right_peak)