lib/peakspot.py

   1 #!/usr/bin/env python
   2
   3 '''
   4 peakspot.py
   5
   6 A library of helping functions for spotting force spectroscopy peaks.
   7
   8 Copyright 2007 by Fabrizio Benedetti and Massimo Sandal
   9
  10 This program is released under the GNU General Public License version 2.
  11 '''
  12
  13 from numpy import mean
  14
  15 def conv_dx(data,vect):
  16     '''
  17     Returns the right-centered convolution of data with vector vect
  18     '''
  19     dim=len(data)
  20     window=len(vect)
  21     temparr=[0.0]*dim
  22
  23     end=dim-window
  24
  25     for j in range(end):
  26         for k in range(window):
  27             temparr[j]+=data[j+k]*vect[k]
  28
  29     return temparr
  30
  31 def absdev(arr):
  32     '''
  33     Calculates the absolute deviation of a vector
  34     '''
  35     med=0.0
  36     absD=0.0
  37
  38     med=mean(arr)
  39     for j in arr:
  40         absD+=abs(j-med)
  41     return absD/len(arr)
  42
  43 def noise_absdev(data,positive=False,maxcut=0.2,stable=0.005):
  44     '''
  45     Returns the standard deviation of the noise.
  46     The logic is: we cut the most negative (or positive) data points until the absolute deviation
  47     becomes stable (it doesn't vary more than 0.005) or we have cut more than maxcut*len(data) points.
  48     Then calculate the absolute deviation.
  49
  50     If positive=True we cut the most positive data points, False=we cut the negative ones.
  51     '''
  52     #TOD: check if this is necessary
  53     out=[item for item in data] #we copy just to be sure...
  54     out.sort()
  55     if positive:
  56         out.reverse()
  57
  58     temp_absdev=absdev(out)
  59
  60     for index in range(len(out)):
  61         cutindex=(index+1)*5
  62         cut_absdev=absdev(out[cutindex:]) #we jump five points after five points...
  63         if 1-(cut_absdev/temp_absdev) < stable and cutindex<(maxcut*len(out)):
  64             temp_absdev=cut_absdev
  65         else:
  66             break
  67
  68     return cut_absdev
  69
  70 def abovenoise(convoluted,noise_level,cut_index=0,abs_devs=4):
  71     '''
  72     Generates a vector which is 0 where the vector is less than abs_devs*noise_level ; 1 if not (spike).
  73     '''
  74     #calculate absolute noise deviation
  75     #noise_level=noise_absdev(convoluted[cut_index:])
  76
  77     above=[]
  78
  79     for index in range(len(convoluted)):
  80         if index<cut_index:
  81             above.append(0)
  82         else:
  83             #FIXME: should calculate the *average* (here we assume that convolution mean is 0, which is FALSE!)
  84             if convoluted[index] < -1*noise_level*abs_devs:
  85                 above.append(convoluted[index])
  86             else:
  87                 above.append(0)
  88     return above
  89
  90 def find_peaks(above, seedouble=10):
  91     '''
  92     Finds individual peak location.
  93     abovenoise() finds all points above a given threshold in the convolution. This point is often not unique
  94     but there is a "cluster" of points above a threshold.
  95     Here we obtain only the coordinates of the largest *convolution* spike for each cluster.
  96
  97     above=vector obtained from abovenoise()
  98     seedouble=value at which we want to "delete" double peaks. That is, if two peaks have a distance
  99     < than $seedouble points , only the first is kept.
 100     '''
 101     nonzero=[]
 102     peaks_location=[]
 103     peaks_size=[]
 104
 105     for index in range(len(above)):
 106         if above[index] != 0:
 107             nonzero.append(index)
 108         else:
 109             if len(nonzero) != 0:
 110                 if len(nonzero)==1:
 111                     nonzero.append(nonzero[0]+1)
 112                 peaks_size.append(min(above[nonzero[0]:nonzero[-1]]))
 113                 peaks_location.append(above[nonzero[0]:nonzero[-1]].index(peaks_size[-1])+nonzero[0])
 114                 nonzero=[]
 115             else:
 116                 pass
 117
 118     #recursively eliminate double peaks
 119     #(maybe not the smartest of cycles, but i'm asleep...)
 120     temp_location=None
 121     while temp_location != []:
 122         temp_location=[]
 123         if len(peaks_location)>1:
 124             for index in range(len(peaks_location)-1):
 125                 if peaks_location[index+1]-peaks_location[index] < seedouble:
 126                     temp_location=peaks_location[:index]+peaks_location[index+1:]
 127         if temp_location != []:
 128             peaks_location=temp_location
 129
 130     return peaks_location,peaks_size