[Tutor] Peaks detection

Mon May 31 04:21:34 EDT 2004

```Marco wrote:

>Hey all,
>i am tryin to write in python a peak detection algorithm to use on spectra taken from a data analyser in a lab.
>The spectrum comes in a text file with 2 columns, one with the channel number and the other one with number of counts for that channel.
>Making the calibration to the energy is not a problem, except for the fact that the algorithm _should_ recognize the peaks of the spectrum (those that have a fair gaussian form) and notate the centroid of the gaussian.
>
>The problem is: i can't find an algorithm except for a labview "procedure" on the net.
>Once found the "how" implementing it in python shouldn't be that hard :)
>
>
>so the question is:
>Anyone has a source for an algorithm of that kind? Any suggestion?
>
>
Hi Marco,

I'm only going to give you an outline of the solution, because this
problem is a lot easier to solve in theory than in practice (code is
untested).

Lets't define two functions:

def is_local_maximum(counts, current_channel):
"""Return True if counts in current_channel are larger than
in neighbouring channels"""

if (counts[current_channel] > counts[current_channel-1] and
counts[current_channel] > counts[current_channel+1]):
return True
else:
return False

def measure_peak(counts, channels, max_channel, width)
"""Return the mean (centroid) and stddev of a peak at max_channel"""

mean = 0.0
stdev = 0.0
total_counts = 0.0
for i in range(max_channel-width, max_channel+width):
total_counts += counts[i]
mean += counts[i]*channels[i]
mean = mean/total_counts
for i in range(max_channel-width, max_channel+width):
stdev += counts[i]*(channels[i]-mean)**2
stdev = math.sqrt(stdev/total_counts)
return mean, stdev

To analyze the dataset we are going to assume that you know the
background signal and the rms of the background signal (measuring them
is left as an exercise)

def analyze(counts, channels, background, background_rms, min_sigma,
window_width):
minimum_signal = background+background_rms*min_sigma
for i in range(window_width, len(channels)-window_width):
if counts[i] < minimum_signal:
# Ignore noise
continue
if is_local_maximum(counts, i):
mean, stdev = measure_peak(counts, channels, i, window_width)
print "Detected peak", counts[i], mean, stdev

Now for the caveats:

- I assume that there is going to be only one peak within the window
around a given maximum
- I assume that the peaks are properly sampled (channel_width <
width_of_peak/2)
- I assume that the window width >> width of the peak
- I assume that the total signal within the window is dominated by the peak
- I assume that the background is flat and the rms of the background is
constant

These assumptions are generally not (all) true in practice

Hope this help,

Roeland

```