hooke/libpeakspot.py

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