hooke/plugin/peakspot.py

   1 # Copyright (C) 2007-2010 Fabrizio Benedetti
   2 #                         Massimo Sandal <devicerandom@gmail.com>
   3 #                         W. Trevor King <wking@drexel.edu>
   4 #
   5 # This file is part of Hooke.
   6 #
   7 # Hooke is free software: you can redistribute it and/or
   8 # modify it under the terms of the GNU Lesser General Public
   9 # License as published by the Free Software Foundation, either
  10 # version 3 of the License, or (at your option) any later version.
  11 #
  12 # Hooke is distributed in the hope that it will be useful,
  13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 # GNU Lesser General Public License for more details.
  16 #
  17 # You should have received a copy of the GNU Lesser General Public
  18 # License along with Hooke.  If not, see
  19 # <http://www.gnu.org/licenses/>.
  20
  21 """a library of helping functions for spotting force spectroscopy peaks.
  22 """
  23
  24 import numpy as np
  25
  26 def conv_dx(data,vect):
  27     '''
  28     Returns the right-centered convolution of data with vector vect
  29     '''
  30     dim=len(data)
  31     window=len(vect)
  32     temparr=[0.0]*dim
  33
  34     end=dim-window
  35
  36     for j in range(end):
  37      for k in range(window):
  38       temparr[j]+=data[j+k]*vect[k]
  39
  40     return temparr
  41
  42 def absdev(arr):
  43     '''
  44     Calculates the absolute deviation of a vector
  45     '''
  46     med=0.0
  47     absD=0.0
  48
  49     med=np.mean(arr)
  50     for j in arr:
  51         absD+=abs(j-med)
  52     return absD/len(arr)
  53
  54 def noise_absdev(data,positive=False,maxcut=0.2,stable=0.005):
  55     '''
  56     Returns the standard deviation of the noise.
  57     The logic is: we cut the most negative (or positive) data points until the absolute deviation
  58     becomes stable (it doesn't vary more than 0.005) or we have cut more than maxcut*len(data) points.
  59     Then calculate the absolute deviation.
  60
  61     If positive=True we cut the most positive data points, False=we cut the negative ones.
  62     '''
  63     out=[item for item in data] #we copy just to be sure...
  64     out.sort()
  65     if positive:
  66         out.reverse()
  67
  68     temp_absdev=absdev(out)
  69
  70     for index in range(len(out)):
  71         cutindex=(index+1)*5
  72         cut_absdev=absdev(out[cutindex:]) #we jump five points after five points...
  73         if 1-(cut_absdev/temp_absdev) < stable and cutindex<(maxcut*len(out)):
  74             temp_absdev=cut_absdev
  75         else:
  76             break
  77
  78     return cut_absdev
  79
  80 def abovenoise(convoluted,noise_level,cut_index=0,abs_devs=4):
  81     '''
  82     Generates a vector which is 0 where the vector is less than abs_devs*noise_level ; 1 if not (spike).
  83     '''
  84     #calculate absolute noise deviation
  85     #noise_level=noise_absdev(convoluted[cut_index:])
  86
  87     above=[]
  88
  89     for index in range(len(convoluted)):
  90         if index<cut_index:
  91             above.append(0)
  92         else:
  93             #FIXME: should calculate the *average* (here we assume that convolution mean is 0, which is FALSE!)
  94             if convoluted[index] < -1*noise_level*abs_devs:
  95                 above.append(convoluted[index])
  96             else:
  97                 above.append(0)
  98     return above
  99
 100 def find_peaks(above, seedouble=10):
 101     '''
 102     Finds individual peak location.
 103     abovenoise() finds all points above a given threshold in the convolution. This point is often not unique
 104     but there is a "cluster" of points above a threshold.
 105     Here we obtain only the coordinates of the largest *convolution* spike for each cluster.
 106
 107     above=vector obtained from abovenoise()
 108     seedouble=value at which we want to "delete" double peaks. That is, if two peaks have a distance
 109     < than $seedouble points , only the first is kept.
 110     '''
 111     nonzero=[]
 112     peaks_location=[]
 113     peaks_size=[]
 114
 115     for index in range(len(above)):
 116         if above[index] != 0:
 117             nonzero.append(index)
 118         else:
 119             if len(nonzero) != 0:
 120                 if len(nonzero)==1:
 121                     nonzero.append(nonzero[0]+1)
 122                 peaks_size.append(min(above[nonzero[0]:nonzero[-1]]))
 123                 peaks_location.append(above[nonzero[0]:nonzero[-1]].index(peaks_size[-1])+nonzero[0])
 124                 nonzero=[]
 125             else:
 126                 pass
 127
 128     #recursively eliminate double peaks
 129     #(maybe not the smartest of cycles, but i'm asleep...)
 130     temp_location=None
 131     while temp_location != []:
 132         temp_location=[]
 133         if len(peaks_location)>1:
 134             for index in range(len(peaks_location)-1):
 135                 if peaks_location[index+1]-peaks_location[index] < seedouble:
 136                     temp_location=peaks_location[:index]+peaks_location[index+1:]
 137         if temp_location != []:
 138             peaks_location=temp_location
 139
 140     return peaks_location,peaks_size