-- Peakfinder Description & example                                  28-Apr-98
-- U. Wienands, Stanford Linear Accelerator Center
-- uli@slac.stanford.edu
~
Peakfinder provides a prodecure "Peak" that locates peaks above a given threshold in
a dataset defined by two vectors, x and y, (e.g., channel and counts). The routine
is most useful in identifying peaks in fft spectra. 

Peak centroids and  peak height are estimated from a parabolic fit through the
three largest y values for a given peak and returned in a matrix Peaks as {centroid,height} pairs, one for each peak found. A vector Peakindex has the indices
of the largest y value for each peak found. Note that, in general, this is NOT exactly
the peak centroid.

Notes:
------
This is a peak finder, not a hump finder. Presence of a peak is determined only by the
data crossing a threshold. It works well with narrow peaks with a well defined maximum,
not well with broad humps. Plot the data before you trust the peaks found.

"Half peaks" at the beginning and end of the data set are ignored.

I have used & tested it for half a year. But, as always, don't trust the results blindly...

User interface:
---------------
Procedure Peak(x,y,threshold) finds peaks & returns them in Peaks[n,2] as [x,y] pairs
                    Peak centroid & heights are from parabolic fit through peak data.

Parameters:     x            The vector of independent values
                y            The vector of dependent values
                threshold    Threshold for peakfinding
Return values:  Peaks        Matrix of {x,y} pairs, one for each peak centroid & height
                Peakindex    Vector of indices for center channel for each peak
~

---------------------------------------------------------------
--Example (this example needs the fft xfun to generate the data):

-- generate some data with noise so we can see how well Peakfinder does
-- Peaks in fft should be at 10.2 and 300.5

dataf[i] = sin(i/1024*360*10.2+rand(0.1))+
           0.5*sin(i/1024*360*300.5+rand(0.1))+
           rand(0.1) dim[1024]
vecf[i] = i dim[1024]
x:=vecf-1:;y:=dataf:

plotline {x,y}

-- fft this data to make peaks

newaxis
Ylogaxis:=on:
pwrspec(y):-- returns spectrum in "transform"

plotline {x[1:513],transform}

-- Find the peaks

threshold:=0.0001:

plotline threshold

Peak(x,transform,threshold):
table Peaks

table Peakindex
-- note that peak at ~10.1 is off by about 0.1. This is 1/10 of the channel width.

-- now lower threshold to see some of the noise peaks & 
-- verify real peaks don't move

threshold:=5*1e-6:

plotline threshold

Peak(x,transform,threshold):
table Peaks

---------------------------------------------------------
----  Procedure Peak(x,y,threshold) ----

Peak(x,y,threshold) = iii:=1, jjj:=1, Peak_ind:=x-x,
  (iii:=iii+1) while y[iii] > threshold, -- skip initial high channels
  ((iii:=iii+1) while (y[iii] < threshold and iii<=count(y)), -- find 1st channel above threshold
   (Peak_ind[jjj]:=iii,Peak_ind[jjj+1]:=iii,
    (iii:=iii+1,
     Peak_ind[jjj+1]:=iii when y[iii]>y[Peak_ind[jjj+1]]) while y[iii] > threshold, -- find peak
    Peak_ind[jjj+2]:=iii-1,
    jjj:=jjj+3) when iii<count(y)) while iii < count(y), jtest:=jjj, -- done
   (jjj:=jjj-3) when (Peak_ind[jjj-1]=count(y) and Peak_ind[jjj-1]=Peak_ind[jjj-2]),, -- correct for end channel
   jmax:=(jjj-1)/3,
  jjj:=1, Peaks:={0[:2]}[:jmax], Peakindex:=0[:jmax],-- now find centroid using parabola fit through peak
(Peakindex[jjj]:=Peak_ind[3*jjj-1],
                Peaks[jjj,1]:=-Peak_b1(x,y,Peak_ind[3*jjj-1])/2/
                Peak_c1(x,y,Peak_ind[3*jjj-1]) +
                x[Peak_ind[3*jjj-1]], -- position of dy/dx=0
 Peaks[jjj,2]:=Peak_a1(y,Peak_ind[3*jjj-1])+
               (Peaks[jjj,1]-x[Peak_ind[3*jjj-1]])*
               Peak_b1(x,y,Peak_ind[3*jjj-1]) +
               (Peaks[jjj,1]-x[Peak_ind[3*jjj-1]])^2*
               Peak_c1(x,y,Peak_ind[3*jjj-1]), -- Crest of parabola
               jjj:=jjj+1) while jjj<=jmax

-- Coefficients of parabola through peak
Peak_a1(y,idx)=y[idx]
Peak_b1(x,y,idx)=(y[idx-1]-y[idx])*(x[idx+1]-x[idx])/
            ((x[idx+1]-x[idx-1])*(x[idx-1]-x[idx])) -
            (y[idx+1]-y[idx])*(x[idx-1]-x[idx])/
            ((x[idx+1]-x[idx-1])*(x[idx+1]-x[idx]))
Peak_c1(x,y,idx)=(y[idx-1]-y[idx])/((x[idx-1]-x[idx])*(x[idx-1]-x[idx+1])) -
   (y[idx+1]-y[idx])/((x[idx+1]-x[idx])*(x[idx-1]-x[idx+1]))